Back to the main page.
Bug 399 - Bug with estimating power of Nyquist freq with ft_specest_mtmfft.m
| Status | CLOSED FIXED |
| Reported | 2011-01-14 00:29:00 +0100 |
| Modified | 2011-03-18 00:37:06 +0100 |
| Product: | FieldTrip |
| Component: | specest |
| Version: | unspecified |
| Hardware: | Macintosh |
| Operating System: | Mac OS |
| Importance: | P1 minor |
| Assigned to: | Roemer van der Meij |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
David Groppe - 2011-01-14 00:29:26 +0100
Created attachment 16 small data set and alt function for computing power spectra ft_specest_mtmfft.m overestimates the power at the Nyquist frequency by a factor of two when the number of time points used for estimation is a factor of two. I believe the problem is that it doesn't take into account the fact that the Nyquist frequency is unique (like the DC component) when the number of time points is a factor of two and should not have its power doubled like the other non-DC frequencies. To illustrate, I've attached a small data set and a function from our lab for estimating power spectra. If you run the following code on the attached data set: load dataPP2 %Estimate power across entire time window cfg = []; cfg.method = 'mtmfft'; cfg.output = 'pow'; cfg.foilim=[0 125]; % from 0 Hz to Nyquist cfg.taper='boxcar'; freq=ft_freqanalysis(cfg,dataPP2); %estimate power with our code [pwr f]=pwr_spec_temp(dataPP2.trial{1}',250,'boxcar',-1); n_time_pts=size(dataPP2.trial{1},2); rect_wind=ones(1,n_time_pts); rect_wind=rect_wind./norm(rect_wind); %normalize taper as is done by ft_specest_mtmfft.m tapered_data=dataPP2.trial{1}.*rect_wind; %apply taper to data fprintf('Sum Sqrs Time Domain: %g, Fieldtrip Sum Pwr: %g, Our Sum Pwr: %g\n', ... sum(tapered_data.^2),sum(freq.powspctrm),sum(pwr)); You should get the following command line output: "Sum Sqrs Time Domain: 139.156, Fieldtrip Sum Pwr: 139.183, Our Sum Pwr: 139.156" Thus fieldtrip is over-estimating the power in the data. The difference is in the Nyquist frequency which is twice as big as what it should be: figure;plot(freq.freq,freq.powspctrm./pwr');
David Groppe - 2011-01-14 00:59:38 +0100
Created attachment 17 non-ftrip function for estimating power spectra I think I forgot to include this in the original bug report.