[Eeglablist] Units in output of timefreq - wavelet normalization

[Makoto Miyakoshi]

Dear Andreas,

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

Is that Chapter 2, Fig 2.25? I found the electric one in our library but
has not page numbers.

I'll comment after reading this Chapter 2.

Makoto



On Thu, Aug 18, 2016 at 8:08 AM, Andreas Widmann <widmann at uni-leipzig.de>
wrote:

> 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
>
> _______________________________________________
> 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
>



-- 
Makoto Miyakoshi
Swartz Center for Computational Neuroscience
Institute for Neural Computation, University of California San Diego
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20160818/5567871b/attachment.html>