Question on Freqanalysis_tfr

[jan-mathijs schoffelen]

Dear Paco,

This looks like a bug to me.
One way of normalizing fft'ed data, is to obtain a so-called spectral
density, i.e. the amount of power per frequency bin. Because you have
data.fsample/2 (=Nyquist frequency) frequency bins, this
normalization would be:

2*(abs(cTmp).^2)/data.fsample.

In the code it seems abs(cTmp)^2 is normalized with the square root
of 2./data.fsample.

Thanks for finding this; I will fix it and the fixed version will be
available in the downloadable toolbox as of tomorrow. However,
probably a good reason why this has gone unnoticed that long, is that
freqanalysis_tfr is a very old piece of code (and very slow), and the
exact same functionality should be covered by freqanalysis_wltconvol
(but this guy is much faster). Even better, freqanalysis_mtmconvol
operates in a very similar way (and you can even get it to mimick a
classical waveletanalysis, given the proper settings), but here you
have much more flexibility in defining your time-frequency resolution
by playing with 'multitapers'.
To give you a look in the kitchen: freqanalysis_tfr operates by
applying a convolution in the time domain, between the time domain
data and the wavelet. Convolution in the time domain is equivalent to
multiplication in the frequency domain. Freqanalysis_mtmconvol/
wltconvol operate by performing this multiplication in the frequency
domain, rather than the slow convolution.

Yours,

Jan-Mathijs

On Feb 16, 2009, at 2:58 PM, Paco Diaz wrote:

> Dear fieldtrip users,
>
>   While reading carefully what the function freqanalysis_tfr do, I
> have
> found something that seems very extrange to me. I think that the
> general
> procedure for the wavelet transformation is clear but I can't get
> why do you
> multiply by 2, and divide by the sampling rate, the absolute value
> of the
> convolution.
>   I have still used the function with very good results, but I
> would like to
> understand that step in order to sleep well.
>
> Here is a piece of the code:
>
>   for k=1:size(dat,1)
>     for j=1:length(cfg.foi)
>       cTmp = conv(dat(k,:),M{j});
>       cTmp = (2*abs(cTmp)/data.fsample).^2;
>       cTmp = cTmp(ceil(length(M{j})/2):length(cTmp)-floor(length(M
> {j})/2));
>       cTmp = cTmp(:,1:cfg.downsample:end);
>       if strcmp(cfg.keeptrials, 'yes')
>         freq.powspctrm(i,k,j,:) = cTmp';
>       else
>         freq.powspctrm(k,j,:) = squeeze(freq.powspctrm(k,j,:)) +
> cTmp';  %
> compute the running sum
>       end
>     end
>   end
>
>
> Thank you very much in advance, Paco.
>
> ----------------------------------
> The aim of this list is to facilitate the discussion between users
> of the FieldTrip  toolbox, to share experiences and to discuss new
> ideas for MEG and EEG analysis. See also http://listserv.surfnet.nl/
> archives/fieldtrip.html and http://www.ru.nl/fcdonders/fieldtrip.

----------------------------------
The aim of this list is to facilitate the discussion between users of the FieldTrip  toolbox, to share experiences and to discuss new ideas for MEG and EEG analysis. See also http://listserv.surfnet.nl/archives/fieldtrip.html and http://www.ru.nl/fcdonders/fieldtrip.