Back to the main page.
Bug 905 - indexing out of range depening on trial-specific number of tapers
| Status | CLOSED DUPLICATE |
| Reported | 2011-08-30 10:32:00 +0200 |
| Modified | 2012-01-11 14:37:05 +0100 |
| Product: | FieldTrip |
| Component: | specest |
| Version: | unspecified |
| Hardware: | Other |
| Operating System: | Linux |
| Importance: | P1 normal |
| Assigned to: | Roemer van der Meij |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Eric Maris - 2011-08-30 10:32:39 +0200
Created attachment 127 test script and small data set that will allow you to recreate the error (see bug description) Dear developers, Over the passed days I stumbled over a number of bugs (or errors from my part) that I at first tried to circumvent by changing some variables' values in my script. However, it didn't become much better. So, here they come. 1. Running ft_freqdescriptives with cfg.trials=[1:ntrials] (i.e., selecting all trials in the freqanalysis output that I pass as an argument to freqdescriptives) gives an error in seloverdim (called by ft_selectdata) at line 82 because indexing operation tmp{j}(sel,:,:,:,:,:,:,:) tries to select elements from tmp{1} that go past the size of the first dimension in tmp{j}. 2. ft_freqanalysis produces an unequal number of tapers for different trials that all have the same number of samples. The difference in number of tapers is always 1. 3. ft_freqdescriptives produces an error that is related to the conversion of my fourier input (obtained from ft_freqanalysis) to power values. The error occurs in the subfunction fixcsd of ft_checkdata (see error message below). Again, this has to do with selecting elements (here, using the variable indx that contains integers that go past the size of the first dimension of data.fourierspctrm. This piece of code is only executed if the fourier data have an unequal number of tapers in the different trials, which they do in my case (see point 2). Error in ==> ft_checkdata>fixcsd at 707 tmpdat = data.fourierspctrm(indx,:,:,:); Error in ==> ft_checkdata at 604 data = fixcsd(data, cmbrepresentation, channelcmb); Error in ==> ft_freqdescriptives at 142 freq = ft_checkdata(freq, 'cmbrepresentation', 'sparsewithpow', 'channelcmb', {}); The occurrence of these problems depends on the settings of the ft_freqanalysis configuration, in particular, the value of tapsmofreq (which I set on a value that is computed from the Fourier frequencies for my data length)and maybe also of foilim. best, Eric Maris
Roemer van der Meij - 2011-09-21 17:01:43 +0200
Hey Eric, It turns out that the different number of tapers per trial (which is not supported), are being caused by round-off errors in the time-axis. If you for example do: data.time(1) - data.time(2) you will see slight variations in the order of 10e-11. In the current implementation, the sampling rate is determined within ft_specest_xxx by using the time-axis. I.e., fsample = 1 / (time(2) - time(1)). Because of the above round-off errors, the computed sampling rate slightly differs per trail, and the number of tapers computed is off because of this. If you do e.g. data.time(i) = round(data.time(i) .* 10000) ./ 10000, it works out. The reason we use the time-axis for computing the sampling rate is.... something I forgot :o. JM (cc), do you remember? Issues with the sampling rate computation have been reported before, I think... The other points you mention are related to the different number of tapers per trial, which is not supported at the moment (i.e. the code in ft_freqdescriptives fails if there are a variable number of tapers per trial). When the number of tapers are constant, the script runs without errors. Cheers, Roemer
Jan-Mathijs Schoffelen - 2011-09-22 10:17:34 +0200
Hi guys, I concur with a part of the diagnosis: the problems seem to be caused by the small round-off issues in the individual trials' time axes, causing different numbers of tapers (e.g. comparing trials 3 and 4 in the example data). The number of tapers is a consequence of the input arguments into the double_dpss function, which occurs in line 87 of ft_specest_mtmfft. Double_dpss is a wrapper around MATLAB's dpss function, and only ensures that the input arguments are double precision, because we had some issues with single precision Slepian sequences in the past. To make a long story short, comparing the second input argument into this function on line 87 for trial 3 and 4 (for one of the iterations over different foilims), I get the following values: 9.749999998476257 versus 9.750000003094783, leading to either 19 or 20 tapers. However, I didn't understand why this should pose any problems. At least, if cfg.output='fourierspctrm' the data matrix contains in the first dimension just the concatenated single taper estimate. It shouldn't care about the fact that trial 3 had 19 tapers, and trial 4 had 20 tapers. In the past this just worked as far as I remember for variable length trials. Looking a bit further I noticed, taking the freq-structure from the second freqbinindx, that the dimensionality of the fourierspctrm is 30x28x7, and this is inconsistent with the expected dimensionality of 32x28x7 (trial 4 and 5 having 4 tapers, leading to 8 trials x 3 tapers + 2 trials x 4 tapers = 32 tapers). So, something unexpected (at least to me) seems to happen in ft_freqanalysis, where the single trial spectra (as computed by ft_specest_mtmfft) are collected. Ft_specest_mtmfft seems to return the variable number of tapers spectra allright, i.e. 3x28x7 for trial 3 and 4x28x7 for trial 4. However, when the fourierspctrm is created (line 674 in ft_freqanalysis) this is done incorrectly because the indexing goes wrong. So, having said all this, here's the verdict: 1 I think that ft_freqanalysis should support variable number of tapers, when cfg.method = 'mtmfft' (because it did so in the past and I don't see an objection against it). It puzzles me why this has gone unnoticed so far. 2 Yet in the present case I don't know whether this is what Eric wants, so we also have the round-off issue with respect to the number of tapers which we need to think about. Regarding 1: I have made a fix in ft_freqanalysis which seems to work fine (no subsequent problems in ft_freqdescriptives and ft_selectdata and ft_checkdata). Yet it may be suboptimal because of memory allocation issues: the initial matrix is allocated assuming a fixed number of tapers per trial. Regarding 2: for the time being, (@Eric) if you ensure each trial to have the same time-axis, i.e. datapart.time(1:end) = datapart.time(1); you should at least get the same number of tapers per trial.