[FieldTrip] Padding with mtmfft and mtmconvol

Thu Sep 21 12:36:30 CEST 2017

Don’t worry Vladimir,

Robert and I have discussed these issues at length on several occasions in the past ;)
Best wishes,
JM

On 21 Sep 2017, at 12:29, Vladimir Litvak <litvak.vladimir at gmail.com<mailto:litvak.vladimir at gmail.com>> wrote:

Hi Jan-Mathijs,

Yes, you are right about Robert's example. But if you do:

pad = 10;

fsample = 1000;
time = (1:1000)/fsample;
dat = randn(size(time));

[spectrum1,ntaper1,freqoi1] = ft_specest_mtmfft(dat, time, 'taper', 'hanning', 'pad', pad);

power1 = abs(spectrum1).^2;
power1 = squeeze(power1);

[spectrum2,ntaper2,freqoi2,timeoi2] = ft_specest_mtmconvol(dat, time, 'taper', 'hanning', 'timeoi', mean(time), 'timwin', 0.99, 'pad', pad);

power2 = abs(spectrum2).^2;
power2 = squeeze(power2);

figure
plot(freqoi1, power1);
hold on
plot(freqoi2, power2, 'r');

You will see the problem that I'm talking about. We discussed with Robert yesterday and this is indeed 'a feature' which has to do with the fact that the outputs of mtmfft and mtmconvol have different units. The former is spectral density whereas the latter is spectral power.

Here is what Robert wrote me:

the units of computations (also here) are a known and long-standing issue. I know for a long time that the two have different scaling, but did not think about it for a long time. I recall something like this: To compare TFRs over frequencies, you don't want the bandwidth to affect the estimate. Shorter wavelets have a larger bandwidth, hence the 1/Hz would affect those. E.g. imagine a 10Hz and a 20Hz sine wave, and do a TFR with conventional wavelets: at 20Hz the wavelet is 2x shorter, so the spectral resolution over which the signal(and noise) spreads is different. If you were to compute the TFR in V^2/Hz, the same V at 20Hz would have a different value, because the length of the wavelet affects the 1/Hz.
something related (but nevertheless different) applies to the mtmfft: if you want to estimate broadband activity in a window of 1 second or a window of 2 seconds, you would get different spectral resolutions. The nyquist is the same, but the power gets distributed over more bins between 0 and Fnyquist/2. That would cause the values to appear smaller in the 2-s case. Hence we compute spectral density, which somehow normalizes for this. I never found a really clear explanation, but google got me this https://dsp.stackexchange.com/questions/33957/what-is-the-difference-between-the-psd-and-the-power-spectrum
what confuses me is that power (or variance) is already normalized, i.e. sum of squared values divided by N. So we have energy (which increases with length), power (which does not increase with length), and power density

So one issue is that most people don't know about this including me and possibly you. I think a good solution would be to add an option to specify the output units for all the methods as there might be quite subtle considerations for choosing one over the other as Robert suggests.

Vladimir

On Thu, Sep 21, 2017 at 10:53 AM, Schoffelen, J.M. (Jan Mathijs) <jan.schoffelen at donders.ru.nl<mailto:jan.schoffelen at donders.ru.nl>> wrote:
Hi to all who’s reading along,

Perhaps the two cases will become more similar once the ‘timwin’ is increased in length for the mtmconvol case….

Best wishes,

JM

On 20 Sep 2017, at 16:26, Robert Oostenveld <r.oostenveld at donders.ru.nl<mailto:r.oostenveld at donders.ru.nl>> wrote:

Hi Vladimir,

I suggest that you first start with a simpler case, like this

fsample = 1000;
time = (1:1000)/fsample;
dat = randn(size(time));

[spectrum1,ntaper1,freqoi1] = ft_specest_mtmfft(dat, time, 'taper', 'hanning');

power1 = abs(spectrum1).^2;
power1 = squeeze(power1);

[spectrum2,ntaper2,freqoi2,timeoi2] = ft_specest_mtmconvol(dat, time, 'taper', 'hanning', 'timeoi', mean(time), 'timwin', 0.5);

power2 = abs(spectrum2).^2;
power2 = squeeze(power2);

figure
plot(freqoi1, power1);
hold on
plot(freqoi2, power2, 'r');

Note that these are not the same (albeit similar), which I had expected…

best
Robert

On 20 Sep 2017, at 12:56, Vladimir Litvak <litvak.vladimir at gmail.com<mailto:litvak.vladimir at gmail.com>> wrote:

Dear Fieldtrippers,

I'm looking into an issue of one of SPM users who gets different results when doing TF decomposition compared to computing a spectrum for the same time window. I'm not sure I got to the bottom of it yet but one thing I found is that ft_specest_mtmfft and ft_specest_mtmconvol are affected differently by increasing padding. For short padding the results are similar but with increasing padding there are differences both in offset of the spectrum and its overall shape. See attached images where the top one shows original spectra and the bottom one aligns the lowermost bin to zero.

Is this a bug or a feature?

Below is the script that produces these plots. I could provide the data as well but this could probably be reproduced with any data.

Thanks,

Vladimir

-------------------------------------

pad = 0.5;%1%10

freqoi = 5:45;
timwin = 0.4+0*freqoi;

[spectrum,ntaper,freqoi,timeoi] = ft_specest_mtmconvol(data, time, 'taper', 'hanning', 'timeoi', 1.2, 'freqoi', freqoi,...
    'timwin', timwin, 'tapsmofrq', 1, 'verbose', 0, 'pad', pad);

s1 = squeeze(mean(mean(abs(spectrum), 4), 2));

figure;
subplot(2,1,1)
plot(freqoi, s1);
subplot(2,1,2);
plot(freqoi, s1-s1(1));
%%
[spectrum,ntaper,freqoi] = ft_specest_mtmfft(data, time, 'taper', 'hanning', 'freqoi', freqoi,...
    'tapsmofrq', 1, 'verbose', 0, 'pad', pad);

s2 = squeeze(mean(mean(abs(spectrum), 2), 1));

subplot(2,1,1)
hold on
plot(freqoi, s2, 'r');
subplot(2,1,2)
hold on
plot(freqoi, s2-s2(1), 'r');
<pad0_5.png><pad1.png><pad10.png>_______________________________________________
fieldtrip mailing list
fieldtrip at donders.ru.nl<mailto:fieldtrip at donders.ru.nl>
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip

_______________________________________________
fieldtrip mailing list
fieldtrip at donders.ru.nl<mailto:fieldtrip at donders.ru.nl>
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip

_______________________________________________
fieldtrip mailing list
fieldtrip at donders.ru.nl<mailto:fieldtrip at donders.ru.nl>
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip

_______________________________________________
fieldtrip mailing list
fieldtrip at donders.ru.nl<mailto:fieldtrip at donders.ru.nl>
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20170921/6415344d/attachment-0002.html>