Back to the main page.
Bug 1228 - improve the fsample computation
| Status | CLOSED FIXED |
| Reported | 2011-12-12 18:12:00 +0100 |
| Modified | 2012-02-01 16:53:22 +0100 |
| Product: | FieldTrip |
| Component: | specest |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P3 normal |
| Assigned to: | Roemer van der Meij |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Robert Oostenveld - 2011-12-12 18:12:55 +0100
As reported by Craig (CC), at this time FT does 1./(time(2)-time(1) to compute the sampling rate. This is sensitive to numerical errors, and 1./mean(diff(time)) would be a better convention. This needs to be fixed in the following functions, where the specest functions are the most important ones manzana> grep 'time(2).*time(1)' *.m */*.m */*/*.m besa2fieldtrip.m: data.fsample = 1/(time(2)-time(1)); % time is already in seconds besa2fieldtrip.m: data.fsample = 1/(data.time(2)-data.time(1)); ft_spike_plot_jpsth.m:cfg.winlen = ft_getopt(cfg,'winlen', 5*(jpsth.time(2)-jpsth.time(1))); ft_spike_plot_jpsth.m:sampleTime = jpsth.time(2) - jpsth.time(1); % get the binwidth private/avw_hdr_make.m: hist.exp_date(1:10) = sprintf('%02d-%02d-%04d',datime(3),datime(2),datime(1)); private/specest_nanfft.m:fsample = 1/(time(2)-time(1)); specest/ft_specest_hilbert.m:fsample = 1/(time(2)-time(1)); specest/ft_specest_mtmconvol.m:fsample = 1/(time(2)-time(1)); specest/ft_specest_mtmfft.m:fsample = 1/(time(2)-time(1)); specest/ft_specest_tfr.m:fsample = 1/(time(2)-time(1)); specest/ft_specest_wavelet.m:fsample = 1/(time(2)-time(1));
Roemer van der Meij - 2011-12-12 18:18:23 +0100
So, thus this mean that we will not un-deprecate fsample? Using 1/mean(diff(time)) will of course be more accurate, but it is still not as an accurate or consistent representation as the true fsample (given by the recording system). (I will go ahead and implement the changes)
Roemer van der Meij - 2011-12-12 18:30:27 +0100
(changes committed)
Jan-Mathijs Schoffelen - 2011-12-13 08:43:01 +0100
*** Bug 976 has been marked as a duplicate of this bug. ***
Robert Oostenveld - 2011-12-14 17:36:57 +0100
On 14 Dec 2011, at 14:38, Richter, Craig wrote: The mean diff approach to the fsample determination has fixed my issue. My spectra are silky smooth now. Thanks!!
Craig Richter - 2011-12-19 13:47:30 +0100
Created attachment 202 No Baseline corrected spectrum
Craig Richter - 2011-12-19 13:55:49 +0100
Hey Guys, Sorry about the out of order attachment. I'm having trouble with the fsample fix again. In the attachment you can see a tfr generated using mtmfft. I am using mtmfft rather than mtmconvol for the tfr since mtmconvol was suspected to be creating problems with Granger causality analysis due to its handling of DC, i.e. DC is not zero for each window. I can ensure this with mtmfft, but this required slicing the data manually. I do this with ft_redefinetrial. This is what led to the problems with numerical accuracy, since time(2) - time(1) had a variable result based on the specific part of the original data that the window was defined using ft_redefinetrial. So my tfrs were not smooth. The current fix reduces the size of the errors so they are not perceptible, as can be seen from the tfr, except that these errors do indeed become perceptible when a percent change baseline is used and there is little power in a given band. Here differences in the baseline value, which is an average over the mtmfft spectra that make up the baseline time interval, and the poststimulus window, that are due to numerical error are very apparent. See the second attached tfr. You can see the problem in the DC bin. Since there was no power in the DC bin due to mean subtraction before freqanalysis, the percent change represents only the change in the numerical error, so the effect is greatly amplified. What do you guys suggest for getting around this? I still like using a singular defined fsample. I highly doubt that this will be the last place we see this sort of phenomenon. Best, Craig
Craig Richter - 2011-12-19 13:56:13 +0100
Created attachment 203 TFR with DC numerical error problem
Roemer van der Meij - 2011-12-19 14:24:48 +0100
Hi Craig, I do not entirely understand. What is it specifically in the second attachment that would be due to numerical issues? Is it the background fluctuations you can see in all frequencies? If so, wouldn't those look less smooth in time and frequency if solely due to numerical issues (or at the least with a single gradient on the frequency axis)? The 0 Hz bin looks very noisy indeed. However, it's non-zero-ness could also partly be due to spectra leakage from the bin above it. Also, what taper did you use? Demeaning is done on the raw data, and tapering (in case of mtmfft) is done on the data after demeaning (mtmconvol tapers the wavelets). This also leads to non-zero dc I would say (although the remaining DC would be small).
Craig Richter - 2011-12-19 15:10:35 +0100
It is just the 0 Hz bin. I actually apply a hanning taper and then demean to tackle the issue with the slightly non-zero mean after tapering. I've been doing this on Pascal's advice do deal with some other issues caused by my very short window analysis. In the previous implementation of mtmfft I had the same 'bumpiness' problem across all bins. This occurred at the same time point for all bins based on the calculation of fsample. There was no issue particular to the zero Hz bin. Now you can see this choppy pattern at nearly all time points. This is because the numerical error is present at all time bins due to the calculation of fsample from all the data.time values, which differ in the pattern of accuracy for each data window used for my tfr. If I look at the DC bin without applying a baseline, then it is very smooth, as in the fist image, but after the baseline the choppiness is present. If you look at any of the bins close up, you can see a bumpiness caused by the lack of accuracy in the fsample calculation, but this is such a small percentage of the actual signal that it is not noticeable. In the DC bin, this error is the only component of the signal, so when compared with the baseline, which may also have a slight difference from zero, this will give rise to zero, when the error matches, and positive and negative signals depending on the ratio of the error. This is why you can see that the colors in the baseline occupy the full +/- 1 range.
Jan-Mathijs Schoffelen - 2011-12-19 15:17:47 +0100
Hi Craig, Wouldn't the problem go away if you fool ft_freqanalysis within the for-loop when doing data.time = repmat({[0:nsmp]/fsample},[1 ntrl]) after ft_redefinetrial, but before ft_freqanalysis? Are you sure that the problem is the fsample numerical issue issue? If your windows are becoming very short with respect to the lowest true frequency content in your data I would expect leakage problems in the lowest frequency bin, particularly if you are using cfg.pad > datalength (are you?).
Craig Richter - 2011-12-20 11:47:57 +0100
I've tested further and the fsample fix is good. My issue is with tapering before demeaning. I was trying to address an issue in the low end of the spectrum by reversing the order. It appears to have worked in that now the DC is so flat that I'm seeing some sort of noise after baselining, but this noise is 10e-34 in magnitude. It is equally present in mtmfft and my test fft, based on matlab commands. So the code appears to be operating properly. Sorry for the false alarm!!
Robert Oostenveld - 2011-12-22 09:24:36 +0100
I noticed in the help/description of ft_datatype_raw that % Deprecated fields: % - fsample and in the code for 2011/latest if ~isfield(data, 'fsample') data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); end Given the recent discussions (online and offline, a.o. with Craig and Roemer), I suggest to be lenient for the moment and keep the code as it is. I.e. fsample will be added on the fly if not present. Still, remember that other fieldtrip code should not rely on it. I have changed the specific code and committed: manzana> svn diff utilities/ Index: utilities/ft_datatype_raw.m =================================================================== --- utilities/ft_datatype_raw.m (revision 5068) +++ utilities/ft_datatype_raw.m (working copy) @@ -100,7 +100,7 @@ end if ~isfield(data, 'fsample') - data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); + data.fsample = 1/mean(diff(data.time{1})); end manzana> svn commit utilities/ Sending utilities/ft_datatype_raw.m Transmitting file data . Committed revision 5070.
Roemer van der Meij - 2011-12-22 10:43:01 +0100
Hmmm, I guess that one slipped through the previous grep grep 'time{1}(2).*time{1}(1)' *.m */*.m */*/*.m ft_freqsimulation.m: cfg.fsample = 1/(cfg.time{1}(2) - cfg.time{1}(1)); % determine from time-axis ft_scalpcurrentdensity.m:scd.fsample = 1/(data.time{1}(2) - data.time{1}(1)); ft_spike_triggeredspectrum.m:fsample = 1./ (data.time{1}(2)-data.time{1}(1)); % assumes constancy compat/ft_databrowser_old.m: fsample = 1/(data.time{1}(2)-data.time{1}(1)); private/data2raw.m: data.fsample = 1/(data.time{1}(2)-data.time{1}(1)); private/rejectvisual_channel.m:fsample = 1/(data.time{1}(2) - data.time{1}(1)); private/rejectvisual_summary.m:fsample = 1/(data.time{1}(2) - data.time{1}(1)); private/rejectvisual_trial.m:fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/data2raw.m: data.fsample = 1/(data.time{1}(2)-data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); Shall I change all of these as well to 1./mean(diff(data{1}.time))? fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/private/data2raw.m: data.fsample = 1/(data.time{1}(2)-data.time{1}(1));
Roemer van der Meij - 2011-12-22 10:44:21 +0100
Wooops, something went wrong with the previous copy/paste. The proper list: grep 'time{1}(2).*time{1}(1)' *.m */*.m */*/*.m: ft_freqsimulation.m: cfg.fsample = 1/(cfg.time{1}(2) - cfg.time{1}(1)); % determine from time-axis ft_scalpcurrentdensity.m:scd.fsample = 1/(data.time{1}(2) - data.time{1}(1)); ft_spike_triggeredspectrum.m:fsample = 1./ (data.time{1}(2)-data.time{1}(1)); % assumes constancy compat/ft_databrowser_old.m: fsample = 1/(data.time{1}(2)-data.time{1}(1)); private/data2raw.m: data.fsample = 1/(data.time{1}(2)-data.time{1}(1)); private/rejectvisual_channel.m:fsample = 1/(data.time{1}(2) - data.time{1}(1)); private/rejectvisual_summary.m:fsample = 1/(data.time{1}(2) - data.time{1}(1)); private/rejectvisual_trial.m:fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/data2raw.m: data.fsample = 1/(data.time{1}(2)-data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); fileio/private/ft_datatype_raw.m: data.fsample = 1/(data.time{1}(2) - data.time{1}(1)); utilities/private/data2raw.m: data.fsample = 1/(data.time{1}(2)-data.time{1}(1));
Roemer van der Meij - 2012-01-11 14:36:51 +0100
*** Bug 905 has been marked as a duplicate of this bug. ***