[Eeglablist] Units in output of timefreq - wavelet normalization

[Andreas Widmann]

Dear Niko,

I’m puzzled by this difference since a long time too (and as you have written a book chapter on WT actually I would have hoped you could help resolving this issue ;).

(Morlet) wavelet normalization always appeared somewhat arbitrary to me (as signal amplitude will never be directly reflected across the whole TF plane for peaky spectra/time courses). To my understanding the most common normalization for wavelets is unit energy (and Gabor). The help text for timefreq states that dftfilt3 is "exact Tallon Baudry“. TB (1998, JNeurosci) states that "Wavelets were normalized so that their total energy was 1,…“.

The wavelets produced by dftfilt3 appear to always have an energy of srate (thus they are *not* unit energy normalized?!):
[wavelet,cycles,freqresol,timeresol] = dftfilt3(F, ncycles, srate);
E = sum( abs( wavelet{ 1 } ) .^ 2 )
Consequently, to my understanding the correct „normfactor“ should be sqrt( E ) or better sqrt( srate ).

You might want to confirm by looking at the TF transform of the (real part of the) wavelet itself
wavelet{ 1 } = wavelet{ 1 } / sqrt( srate ); % Unit energy normalize
mysig = zeros(1,D*srate);
delay = ceil( ( length( mysig ) - length( wavelet{ 1 } ) ) / 2 );
mysig( delay:delay + length(  wavelet{ 1 } ) - 1 )  = real( wavelet{ 1 } ) * 2; % Discard imag part
normfactor = sqrt( srate );
which should now have a peak amplitude of 1 (independent of sampling rate and signal duration).

As the issue appears to be not only in EEGLAB but also other implementations, I always assumed my reasoning to be incorrect. Is it? What am I missing?

Best,
Andreas

> Am 10.08.2016 um 11:58 schrieb Niko Busch <niko.busch at gmail.com>:
> 
> Dear Makoto (and everyone who replied to me personally regarding this post),
> 
> thank you for your reply! I see that the result of the wavelet transform inside the timefreq function is dependent on the length of the signal, which in turn is dependent on the number of cycles and sampling rate. However, simply dividing by the length of the wavelet does not seem to be the solution either. I modified the code below by including a "normalization factor", which currently is simply the length of the wavelet. Dividing the wavelet transformed amplitudes by this factor gives the right order of magnitude, but the results are still quite off. By increasing the sampling rate or number of cycles, the results are even more off. I believe we are on the right track, but something is still missing. Do you have any ideas?
> 
> Cheers,
> Niko
> 
> %% Create sine wave
> clear all
> D = 4;       % total signal duration in seconds.
> sigD = 1;    % duration of the test oscillation within the signal.
> F = 10;      % frequency of the test oscillationin Hz.
> P = .25;     % Phase of the test oscillation. 2 pi radians = 360 degrees
> srate = 256; % sampling rate, i.e. N points per sec used to represent sine wave.
> T = 1/srate; % sampling period, i.e. for this e.g. points at 1 ms intervals
> t = [T:T:D]; % time vector.
> 
> sigpoints = length(t)/2 - (sigD*srate)/2:(length(t)/2 + (sigD*srate)/2)-1;
> mysig = zeros(1,D*srate);
> mysig(sigpoints) = sin(2*F*t(sigpoints)*pi+ 2*pi*P);
> 
> %% TF computation
> ncycles = 4;
> 
> [wavelet,cycles,freqresol,timeresol] = dftfilt3(F, ncycles, srate);
> normfactor = length(wavelet{1});
> 
> [tf, outfreqs, outtimes] = timefreq(mysig', srate, ...
>     'cycles', ncycles, 'wletmethod', 'dftfilt3', 'freqscale', 'linear', ...
>     'freqs', F);
> 
> %% Plot
> figure; hold all
> plot(t,mysig); 
> plot(outtimes./1000,abs(tf)./normfactor)
> xlabel('Time (seconds)');
> ylabel('Amplitude');
> legend('input signal', 'wavelet result')
> 
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to eeglablist-request at sccn.ucsd.edu