[Eeglablist] Units in output of timefreq - wavelet normalization

Andreas Widmann widmann at uni-leipzig.de
Thu Aug 18 08:08:16 PDT 2016 

Dear Norman,

> The other points you mentioned, I found super super interesting and helpful and it would be great if you could comment on some emerging questions. Maybe it is also interesting for others to approach this intricate topic on a more detailed and basic level.
> 
> You mentioned the gabor transform is rather suited for narrowband signals and EEG, however, is closer to white noise. But to my eyes EEG signals seem to be a lot more 'structured' e.g. when one thinks of the narrow spectral peaks of the rolandic rhythm in the alpha and beta band, or when a steady state evoked potential paradigm was applied.
I would like to split the question into two parts: (1) Narrowband signals and wavelet transform and (2) EEG.

(1) Problematic for narrowband signals is actually the variable frequency resolution of the Morlet wavelet transform. If you consider a signal composed of two sines, e.g.

f=[ 10 50 ];
signal = sin( 2 * pi * f( 1 ) * t )' + sin( 2 * pi * f( 2 ) * t )’;

with my script from yesterday you will notice that normalization analog to Gabor transform preserves peak amplitude while energy normalization preserves energy. As energy is spread wider for higher frequencies there is either a drop of peak amplitude or overestimation of total energy, both suboptimal solutions.

You can keep frequency (and time) resolution constant by increasing the number of cycles with frequency. This is, however, in the strict sense no longer a wavelet transform (mother wavelet is scaled, that is, number of cycles is fixed) but the Gabor transform (Gaussian window; a special case of STFT). There is an excellent and accessible comparison of Morlet wavelet vs. Gabor transform in [1, pp. 45-51]. Personally, I would likely use Gabor transform for the analysis of this kind of narrowband test signals. It is an interesting question whether Morlet wavelet vs. Gabor transform would be better suited for EEG data. I would like to avoid commenting here as for my own data I would actually neither use Morlet or Gabor transform but the multitaper analysis as implemented in Fieldtrip, which gives excellent control over frequency resolution (in particular for high frequencies) and other aspects.

(2) I did not mean (or write) white noise but broadband and noisy ;) Besides the 1/f spectrum the EEG certainly has regular structures, yes. I fully agree. To my experience they are, however, not as sharp/regular as a sinusoid but more broadband due to fluctuations of amplitude, frequency, and/or phase (we want to detect in tf analysis).

> The gabor wavelet spectrum indeed seems to produce higher average amplitudes for higher frequencies, but the actual average values seem to vary much less than for the unit energy normalization, which then certainly does not produce 'flat' average amplitudes across frequencies. Do you think this is related to scaling issues?
Yes, i think this is a scaling issue.

> Unit energy normalization seems to increase frequency smearing for higher frequencies. Although one may apply a baseline normalization later, do you think this makes unit energy normalization less comparable to gabor normalization (i.e. contrary to what I've written above?)
No, the normalization does not matter for smearing (frequency resolution).

> You said it doesn’t matter whether you normalize in the time or frequency domain. I'm confused now about this line of unit energy code:
> A = sqrt( pnts ./ ( sqrt( pi ) * sigma_f ) );
> Could you explain how sqrt(srate) normalization term transformed to the line above?
Frequencies are normalized in the first step for convenience. Thus, the srate term is indeed missing here, sorry. Otherwise the terms should be identical for time and frequency domain (replacing sigma_f by sigma_t). I always manually check for unit energy by sum( abs( wavelet ) .^ 2 ) which should result in 1 in the time and pnts in the frequency domain.

> Also, when plotting the envelope based on this normalization, it yields an amplitude which is 25 times higher than in the original signal. When I normalized my complex wavelet with sqrt(srate)
>> m =   exp(-t.^2/(2*s^2)).*exp(1i*2*pi*F.*t) ./ sqrt(srate);
> 
> as you did before for the dft3filt wavelet
>> wavelet_norm{ 1 } = wavelet{ 1 } / sqrt( srate );
> 
> I got an envelope which was just 2.5 times higher than the original signal. So I guess there's an error in my code. could you elaborate on that please?
I cannot replicate „2.5 times higher“, sorry. I did notice two things: you changed cycles from 4 to 7. And we used amplitude in Niko’s and your first but power in my and your second example. Could this possibly explain the differences? With 4 cycles you get a peak amplitude of 3.8 (squared 14.4), or with 7 cycles of 5 (squared 25.3) with both, Niko’s (with normfactor = sqrt(srate)) and your new code.

I would like to note that the amplitudes you get with unit energy normalization are not as arbitrary as they appear. They do reflect the energy in the signal covered by the Heisenberg box [1]. I selected Gaussian white noise for my example to demonstrate this (but I messed it up by overwriting a variable in the first attempt, sorry again). With a higher number of cycles the Heisenberg box is "longer" (higher frequency, lower time resolution) and covers a larger part of the sinusoidal test signal.

Best,
Andreas

[1] Addison, P.S. (2002). The illustrated wavelet transform handbook. IOP Publishing, Bristol.

> I am very much looking forward to your comments, thanks a lot Andreas!
> Best
> Norman
> 
> figure links:
> .fig: https://www.dropbox.com/s/ktoffojfcxxu6vd/wavelet_normalization.fig?dl=0
> .png: https://www.dropbox.com/s/v5i4d95y53t2woi/wavelet_normalization.png?dl=0
> 
> clear
> 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.
> dt = 1/srate; % sampling period, i.e. for this e.g. points at 1 ms intervals
> time = dt:dt:D; % time vector.
> 
> sigpoints = length(time)/2 - (sigD*srate)/2:(length(time)/2 + (sigD*srate)/2)-1;
> mysig = zeros(1,D*srate);
> mysig(sigpoints) = sin(2*F*time(sigpoints)*pi+ 2*pi*P);
> 
> figure; hold all
> plot(time,mysig);
> 
> cycles = 7;
> pnts = D * srate;
> 
> freqs = F * ( pnts / srate );
> sigma_f = repmat( freqs / cycles, [ pnts 1 ] );
> f = repmat( 0:pnts - 1, [ size( freqs, 2 ) 1 ] )' - repmat( freqs, [ pnts 1 ] );
> cmorlf = exp( -f .^ 2 ./ ( 2 * sigma_f .^ 2 ) );
> 
> % Analog to Gabor transform
> A = 2;
> cmorlf_A = A .* cmorlf;
> 
> tf = ifft( fft( mysig' ) .* cmorlf_A );
> 
> plot( time, abs( tf ) .^ 2 , 'r-')
> 
> % Unit energy
> 
> A_unit = sqrt( pnts ./ ( sqrt( pi ) * sigma_f ) );
> cmorlf_Au = A_unit .* cmorlf;
> 
> tf_unit = ifft( fft( mysig' ) .* cmorlf_Au );
> 
> plot( time, abs( tf_unit ) .^ 2, 'm:' )
> legend({'Nico''s signal';'gabor';'unit energy'})
> 
> 
>> ----- On Aug 17, 2016, at 3:05 PM, Andreas Widmann widmann at uni-leipzig.de wrote:
>> 
>>> Dear Norman,
>>> 
>>> sorry, I do not fully agree. First, I think it is important not to mix up the
>>> two questions:
>>> 
>>> (1) Does dftfilt3 provide the expected/documented results?
>>> To my understanding this is not the case. I hope I do not miss anything obvious.
>>> 
>>> vs. (2) What is the recommended/best/optimal wavelet normalization for EEG data
>>> analysis?
>>> I think there cannot be a generally valid recommendation and different solutions
>>> are necessary for different applications. The normalization you suggested is
>>> equivalent to the normalization of the Gabor transform (incorrectly abbreviated
>>> to Gabor normalization in my previous post, sorry). This normalization is well
>>> suited for narrowband signals and peaky spectra and the frequently used for
>>> example in audio signal analysis. For broadband/noise signals, however, the
>>> Gabor transform normalization overestimates high frequencies relative to low
>>> frequencies.
>>> 
>>> Please find below some code directly comparing both normalizations using a white
>>> noise signal (btw. it doesn’t matter whether you normalize in the time or
>>> frequency domain). In the global wavelet spectrum (middle column) power is
>>> increasing with frequency for Gabor transform normalization (top row) but
>>> „flat“ for unit energy normalization (bottom row).
>>> 
>>> As EEG reflects a broadband and noisy signal I personally usually prefer to
>>> apply energy normalization. As always, there are exceptions and other opinions
>>> and special applications and ...
>>> 
>>> Best,
>>> Andreas
>>> 
>>> fs = 256;
>>> T = 4;
>>> pnts = T * fs;
>>> t = ( 0:pnts - 1 ) / fs;
>>> cycles = 7;
>>> F = 2:2:70;
>>> 
>>> % White noise
>>> rng( 0 )
>>> signal = randn( fs * T, 1 );
>>> 
>>> % Morlet wavelet
>>> freqs = F * ( pnts / fs );
>>> sigma_f = repmat( freqs / cycles, [ pnts 1 ] );
>>> f = repmat( 0:pnts - 1, [ size( freqs, 2 ) 1 ] )' - repmat( freqs, [ pnts 1 ] );
>>> cmorlf = exp( -f .^ 2 ./ ( 2 * sigma_f .^ 2 ) );
>>> 
>>> figure
>>> 
>>> % Analog to Gabor transform
>>> A = 2;
>>> cmorlf = A .* cmorlf;
>>> 
>>> tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* cmorlf );
>>> 
>>> h = subplot( 2, 3, 3 );
>>> imagesc( t, F, abs( tf' ) .^ 2 )
>>> set( h, 'YDir', 'normal' );
>>> subplot( 2, 3, 2 );
>>> plot( F, mean( abs( tf ) .^ 2 ) )
>>> subplot( 2, 3, 1 )
>>> plot( (0:pnts / 2) * fs / pnts, cmorlf( 1:pnts / 2 + 1, : ) )
>>> 
>>> % Unit energy
>>> A = sqrt( pnts ./ ( sqrt( pi ) * sigma_f ) );
>>> cmorlf = A .* cmorlf;
>>> 
>>> tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* cmorlf );
>>> 
>>> h = subplot( 2, 3, 6 );
>>> imagesc( t, F, abs( tf' ) .^ 2 )
>>> set( h, 'YDir', 'normal' );
>>> subplot( 2, 3, 5 );
>>> plot( F, mean( abs( tf ) .^ 2 ) )
>>> subplot( 2, 3, 4 )
>>> plot( (0:pnts / 2) * fs / pnts, cmorlf( 1:pnts / 2 + 1, : ) )
>>> 
>>>> Am 15.08.2016 um 18:19 schrieb Norman Forschack <forschack at cbs.mpg.de>:
>>>> 
>>>> Dear all,
>>>> 
>>>> I'd like to contribute from the perspective of a discussion on Mike Cohen's
>>>> Blog.
>>>> 
>>>> The initial question was, how to obtain an amplitude envelope of a given signal
>>>> which has in fact the same amplitude as the given signal, right?
>>>> So coming from Nicos signal:
>>>> 
>>>> 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
>>>> time = T:T:D; % time vector.
>>>> 
>>>> sigpoints = length(time)/2 - (sigD*srate)/2:(length(time)/2 + (sigD*srate)/2)-1;
>>>> mysig = zeros(1,D*srate);
>>>> mysig(sigpoints) = sin(2*F*time(sigpoints)*pi+ 2*pi*P);
>>>> 
>>>> one way to obtain equal amplitudes is to normalize the wavelet by its maximal
>>>> value within the frequency domain:
>>>> 
>>>> % some preparations
>>>> mysig = mysig';
>>>> ss = size(mysig);
>>>> cycles = 4;
>>>> dt = 1/srate;
>>>> sf = F/cycles;
>>>> s = 1./(2*pi*sf);
>>>> t = (-4*s:dt:4*s)';
>>>> nM = length(t);
>>>> halfMsiz = (nM-1)/2;
>>>> hz = linspace(0,srate/2,floor(nM/2)+1);
>>>> Ly = ss(1)*ss(2)+nM-1;
>>>> Ly2=pow2(nextpow2(ss(1)*ss(2)+nM-1));
>>>> 
>>>> % fft of signal
>>>> X=fft(reshape(mysig,ss(1)*ss(2),1), Ly2);
>>>> 
>>>> % building morlet wavelet (without a normalization factor)
>>>> m =   exp(-t.^2/(2*s^2)).*exp(1i*2*pi*F.*t);
>>>> H = fft(m,Ly2); % fft of wavelet
>>>> 
>>>> % normalize wavelet spectrum
>>>> H = H./max(H);
>>>> 
>>>> y =  ifft(X.*H,Ly2);
>>>> y = y(floor(halfMsiz+1):Ly-ceil(halfMsiz));
>>>> y_amp = 2* abs(y);
>>>> 
>>>> figure; plot(time,mysig,'b',time,y_amp,'r')
>>>> 
>>>> This seems to work for any combination of srate and cycles (except when number
>>>> of cycles become large) because the signal spectrum is convolved by spectral
>>>> wavelet values being maximally one.
>>>> 
>>>> I have not fully worked my way through Andreas' example. It normalizes the
>>>> wavelet in time, not in frequency domain as here. So it is probably not
>>>> comparable.
>>>> But doing time domain normalization within the lines above by just replacing the
>>>> kernel formula:
>>>> 
>>>> m =   exp(-t.^2/(2*s^2)).*exp(1i*2*pi*F.*t) ./ sqrt(srate); % unit energy
>>>> 
>>>> and commenting out the max(H) normalization, however, yields an amplitude
>>>> envelope which is 2.5 times larger than the original signal amplitude and
>>>> increases when the number of wavelet cycles is increased. But as Andreas
>>>> suggestion referred to the dftfilt3 output, the matter becomes more complicated
>>>> as this function uses it's own normalization factor:
>>>> A = 1./sqrt(s*sqrt(pi));
>>>> and there are problably some more relevant differences (not even going into the
>>>> timefreq function).
>>>> 
>>>> In sum, this post may have fostered the general confusion (or at least mine) but
>>>> for a more puristic approach to the matter of wavelet normalization, the lines
>>>> above might be of some value (kudos to mike x cohen, of course).
>>>> 
>>>> All the best
>>>> Norman
>>>> 
>>>> 
>>>> ----- On Aug 12, 2016, at 6:06 PM, Andreas Widmann widmann at uni-leipzig.de wrote:
>>>> 
>>>>> 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
>>>>> 
>>>>> _______________________________________________
>>>>> 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
>>>> _______________________________________________
>>>> 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
> _______________________________________________
> 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

More information about the eeglablist mailing list