[FieldTrip] Padding with mtmfft and mtmconvol

Thu Sep 21 12:29:12 CEST 2017

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> 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>
> 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>
> 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
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>
>
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>
>
>
> _______________________________________________
> fieldtrip mailing list
> 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/aa6c8cff/attachment.html>