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.