Back to the main page.
Bug 2365 - Fourier output of time-frequency methods (mtmconvol, wavelet) does not work with variable-length trials
| Status | CLOSED FIXED |
| Reported | 2013-11-05 22:01:00 +0100 |
| Modified | 2015-07-15 13:31:23 +0200 |
| Product: | FieldTrip |
| Component: | core |
| Version: | unspecified |
| Hardware: | All |
| Operating System: | All |
| Importance: | P1 normal |
| Assigned to: | Jan-Mathijs Schoffelen |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
- 2013-11-05 22:01:57 +0100
Created attachment 555 bug data Run attached code on attached data (5 trials worth, variable length) and note that cfg.output = 'fourier' does not work due to assignment mismatch issues. Data is not getting padded correctly (or something) if you inspect the assignment 'fourierspctrm(currrptind,:,:,:) = spectrum;' in the debugger. Same result for both mtmconvol and wavelet methods. The code works fine with cfg.output = 'pow', so it is probably something wrong with the fftflag or keeprpt = 4 code branches in ft_freqanalysis.
- 2013-11-05 22:02:18 +0100
Created attachment 556 bug script
- 2013-11-05 22:05:05 +0100
I should also note that cfg.output = 'fourier' works fine with method mtmfft (since the trials are averaged over), so it is specifically the time-frequency methods that return fourier spectra are broken for variable-length trials.
Robert Oostenveld - 2013-11-06 11:50:16 +0100
Hi Dean, thanks for the detailed bug report, we will follow this up asap.
Haiteng Jiang - 2013-11-11 13:42:33 +0100
(In reply to comment #0) (In reply to comment #3) Hi Robert & Dean , I just had a look . Indeed , the current ft_freqanalysis doesn't deal with different trial length very well. This problem can be solved if we initialize the fourierspctrm with max trial length . Accordingly , we need to change ft_freqanalysis(cfg, data) at line 533; ntoi = numel(toi) --> ntoi = max(trllength); then update the spectrm matrix with different trial length line 702 spectrum = permute(reshape(spectrum_mtmconvol,[nchan ntoi ntaper(1) nfoi]),[3 1 4 2]); --> spectrum = permute(reshape(spectrum_mtmconvol,[nchan trllength(itrial) ntaper(1) nfoi]),[3 1 4 2]); line 713 fourierspctrm(currrptind,:,:,:) = spectrum; --> fourierspctrm(currrptind,:,:,1:trllength(itrial)) = spectrum; This will influence output 'pow ' as well, so we have to adapt : line 606 spectrum = reshape(permute(spectrum_mtmconvol(:,:,freqtapind{ifoi}),[3 1 2]),[ntaper(ifoi) nchan 1 ntoi]); --> spectrum = reshape(permute(spectrum_mtmconvol(:,:,freqtapind{ifoi}),[3 1 2]),[ntaper(ifoi) nchan 1 trllength(itrial)]); @ Robert: I don't know debug the bug on the bug on our meta since I don't have edit right. At this moment , I do it on my desktop, which is kind of inconvenient . Coud you please point me the common way to do it in DCCN ? Thanks in advance.
Robert Oostenveld - 2013-11-11 14:43:23 +0100
I am surprised to see that the option cfg.pad is not being used in the test script. But the default is maxperlen, which should be ok. When I run the script, I see >> convol = ft_freqanalysis(cfg, data); the input is raw data with 128 channels and 5 trials processing trials trial 2, frequency 37 (37.06 Hz), 1 tapers Subscripted assignment dimension mismatch. Error in ft_freqanalysis (line 713) fourierspctrm(currrptind,:,:,:) = spectrum; 713 fourierspctrm(currrptind,:,:,:) = spectrum; where K>> size(spectrum) ans = 1 128 40 1150 K>> size(fourierspctrm) ans = 5 128 40 1151 there is a one-sample difference. Looking at data.time{1} and data.time{2} I do see that they differ by one sample. But then I also see that cfg.toi = 'all' which is not a supported option according to the documentation of ft_freqanalysis itself. The option is supported by the low-level code, but not the high level code. If I change the test script into cfg.toi = 0:0.1:5 cfg.pad = 'maxperlen' then it computes just fine. So it is not a code error, but an (uncaught) user error. Passing 'all' to the low-level code (ft_specest_mtmconvol) results in that function indeed returning "all" timebins, which are different in trial 1 and 2. Should we explicitly check for this user error?
- 2013-11-11 21:45:56 +0100
I had forgotten that cfg.toi = 'all' was unsupported -- I found it by looking through the source. Is the preferred calling convention for these methods then to find the longest trial and set cfg.toi accordingly to do the analysis on all timepoints? e.g., maxlen = max(cellfun(@length, data.time)); cfg.toi = data.time{cellfun(@length, data.time) == maxlen}; I can confirm that this works on my actual data, but it does seem like it would be easier to make a cfg.toi = 'all' option, to save user the trouble of having to do this (and maybe make it default, to match the default maxperlen padding). Thanks Robert and Haiteng for working on this!
Jan-Mathijs Schoffelen - 2015-02-11 11:08:48 +0100
Since Haiteng won't be following up on this one, I will take over this one in order to get it resolved. @Haiteng, if you still feel like contributing to this relatively straightforward one, it's most welcome.
Jan-Mathijs Schoffelen - 2015-02-11 14:51:02 +0100
looking into this a bit more I noticed that: -ft_specest_mtmconvol supports 'all' as a timeoi specification -the issue occurs when requesting 'fourier' in the output because it runs a different chunk of code in ft_freqanalysis. When asking for cfg.output = 'pow' with cfg.keeptrials results in some extra checking in ft_freqanalysis.
Jan-Mathijs Schoffelen - 2015-02-11 15:15:26 +0100
[jansch@mentat003 fieldtrip]$ svn diff ft_freqanalysis.m Index: ft_freqanalysis.m =================================================================== --- ft_freqanalysis.m (revision 10200) +++ ft_freqanalysis.m (working copy) @@ -253,6 +253,16 @@ error('you must specify a smoothing parameter with taper = dpss'); end end + cfg = ft_checkconfig(cfg, 'required', 'toi'); + if ischar(cfg.toi) && strcmp(cfg.toi, 'all') + % do an educated guess with respect to the requested time bins + begtim = min(cellfun(@min,data.time)); + endtim = max(cellfun(@max,data.time)); + cfg.toi = linspace(begtim,endtim,round((endtim-begtim)./mean(diff(data.time{1})))+1); + + elseif ischar(cfg.toi) + error('cfg.toi should be either a numeric vector, or can be ''all'''); + end case 'mtmfft' cfg.taper = ft_getopt(cfg, 'taper', 'dpss'); [jansch@mentat003 fieldtrip]$ svn add test/test_bug2365.m A test/test_bug2365.m [jansch@mentat003 fieldtrip]$ svn commit -m "bugfix - added check for cfg.toi = 'all'" ft_freqanalysis.m test/test_bug2365.m Sending ft_freqanalysis.m Adding test/test_bug2365.m Transmitting file data .. Committed revision 10202. [jansch@mentat003 fieldtrip]$
Jan-Mathijs Schoffelen - 2015-05-07 13:20:39 +0200
regression function is failing due to a file permissions issue.