[Eeglablist] Units in output of timefreq - wavelet normalization

[Makoto Miyakoshi]

Dear Jason,

> related to this discussion

Andreas recommended the following book in the discussion above. I read the
Chapter 2, and it is indeed very much related, I thought.

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

I've also read a chapter of Mike X Cohen's book, but not about wavelet.
It's a very good book, many of my colleagues including myself were
impressed!

Makoto

On Wed, Aug 24, 2016 at 12:01 AM, jason roger <jasonroger8 at gmail.com> wrote:

> Dear Andreas, Norman, Makoto, and EEGLAG list,
>
> This is very interesting discussion, can of you highlight strong books
> that are related to this discussion and major EEG information with their implementation in
> Matlab?
>
> Highly appreciate any help.
>
> Jason
>
>
>
> On Wed, Aug 24, 2016 at 8:47 AM, Makoto Miyakoshi <mmiyakoshi at ucsd.edu>
> wrote:
>
>> Dear Norman,
>>
>> Thanks for your comment. I would like to confirm the following with you.
>>
>> About (6),
>>
>> > But in fact, there are a lot of possibilities how steps of these cycle
>> numbers, and therefore the window widths, are defined across frequencies.
>>
>> But technically speaking, you vary the cycle numbers and it's not wavelet
>> anymore, correct? I understand there are several reasonable ways to adjust
>> the cycles numbers, but that's not the main point here.
>>
>> It does not matter even if newtimef() is actually not a wavelet
>> transform, if we use is with varying number of cycles. I just want to give
>> it a correct label and description that points us to exact methods with
>> which engineers and specialists in other fields can understand.
>>
>> About (3), yes, if we see flat frequency responses it means
>> underweighting the high-frequency power. That seems correct intuitively.
>> However, raw output normalization is mainly for convenience anyway, and
>> there could be more than one user need in doing this. It does not hurt to
>> support multiple ways to normalize it, unless we all see that one method is
>> apparently wrong. The energy normalization, as I understand, is not a
>> 'wrong' method (or is it?) but one of reasonable and potentially useful
>> ways to normalize the output--across frequency too.
>>
>> Thank you again for your continued interest and contribution.
>> If we can all agree with the summary, I'll feed it to EEGLAB developers.
>> It was indeed one of the deepest discussion on the list in the past years!
>>
>> Makoto
>>
>>
>>
>> On Tue, Aug 23, 2016 at 3:19 PM, Norman Forschack <forschack at cbs.mpg.de>
>> wrote:
>>
>>> Dear Makoto and list,
>>>
>>> I agree to almost all points with which Makoto summarized our discussion
>>> so well. Yet, after reading this excellent reference provided by Andreas, I
>>> like to further distinguish things within points 3 and 6. Please also
>>> correct me if I'm wrong.
>>>
>>> Starting with 6, I wouldn't say that if the number of cycles in wavelet
>>> increases over frequency, this actually is a Gabor transform/STFT per se.
>>> According to Addison, there are at least two differences. For Gabor 'the
>>> internal frequency, f, is allowed to vary within the gaussian window of
>>> fixed width' (qoting Addison). In other words, time and frequency
>>> resolution stays the same across analysed frequencies once the window
>>> length is given. For wavelets it's possible to 'immitate' this by applying
>>> variable cycle numbers across frequencies (with which it's not a wavelet
>>> transform anymore as Andreas noted). But in fact, there are a lot of
>>> possibilities how steps of these cycle numbers, and therefore the window
>>> widths, are defined across frequencies. I think the most convenient control
>>> for resolution parameters can be gained by ''Define empirical FWHM of
>>> wavelet in the frequency domain' as Mike proposed' (point 4). This would,
>>> however, still require that wavelet amplitude is normalized to unit gain
>>> (i.e. maximum value normalization and not unit energy or, alternatively,
>>> unit energy plus 'normalized to the RMS amplitude of the sinusoid').
>>> Another difference is, 'that the wavelet complex sinusoid is centered at
>>> [the given time point] on the time axis whereas the complex sinusoid
>>> contained within the gabor atom is 'centered' at the origin (t=0). [...So
>>> for] STFT sinusoid remains fixed in relation to the origin [...] whereas
>>> for the wavelet [...] the sinusoid and window move together. This should
>>> get important if you want to extract phase related measures which, without
>>> a correction term, wouldn't be comparable between the two transforms.
>>>
>>> For point 3, I think we haven't reached a consensus yet. I think Mike
>>> was referring to unit energy normalization when he stated: 'Trying to
>>> correct for this [maintain constant power across frequencies] while keeping
>>> the number of cycles constant is actually underestimating power at higher
>>> frequencies, because it involves trying to make the red bars more blue,
>>> even though the redness is accurate.' (if you're still listening Mike, is
>>> that right?) So it might be mentioned that applying unit energy
>>> normalization ensures that for an analysis with fixed cycle number across
>>> frequencies, the area covered by wavelet spectra (energy) is the same for
>>> all analysed center frequencies (but variable cycles is a different thing).
>>> (or maybe there's a more elegant way to put this).
>>>
>>> I also want to thank all the folks contributing to this discussion, I
>>> learned a lot from you!
>>> Best
>>> Norman
>>>
>>>
>>> ----- Original Message -----
>>> From: "Makoto Miyakoshi" <mmiyakoshi at ucsd.edu>
>>> To: "Andreas Widmann" <widmann at uni-leipzig.de>, "Norman Forschack" <
>>> forschack at cbs.mpg.de>, "Niko Busch" <niko.busch at gmail.com>, "Mike X
>>> Cohen" <mikexcohen at gmail.com>
>>> Cc: "eeglablist" <eeglablist at sccn.ucsd.edu>, ramon at sccn.ucsd.edu
>>> Sent: Samstag, 20. August 2016 02:21:12
>>> Subject: Re: [Eeglablist] Units in output of timefreq - wavelet
>>> normalization
>>>
>>>
>>> Update--
>>>
>>>
>>>
>>> Ok, let me confirm conclusions for the EEGLAB developers.
>>>
>>>
>>> 1) This is the bottom line: Let's not use the raw uV^2 outputs since
>>> these changes due to non-interesting factors as Mike mentioned. Instead,
>>> let's use dB changes against baseline etc.
>>>
>>>
>>> 2) If the raw uV^2 is absolutely necessary, there are probably two
>>> reasonable ways to normalize the raw values. One is, in Andreas's sample
>>> code, 'Analog to Gabor transform' in which A = 2. You'll find increasing
>>> power as frequency increases (which is the nature of the wavelet and not a
>>> problem.)
>>>
>>>
>>> 3) The other good one is 'Unit energy' in which A =
>>> sqrt(pnts./(sqrt(pi)*sigma_f). Probably, EEGLAB developers can implement
>>> both as options and leave the choice to the users. You'll find flatter
>>> response over the frequency range.
>>>
>>> 4) Yet another good one is 'Define empirical FWHM of wavelet in the
>>> frequency domain' as Mike proposed. Personally, I found it most reasonable
>>> for intuitiveness. It's also normalized across frequencies.
>>>
>>>
>>> 5) EEGLAB's explanation about dftfilt3() needs to be at least more
>>> detailed, and potentially the default/optional normalization can be changes
>>> to one of the above two, as long as it does not conflict with backword
>>> compatibility.
>>>
>>>
>>> 6) Also, it's probably good to add the following description to EEGLAB
>>> functions that if the number of cycles in wavelet increases in the time
>>> domain as the frequency increases, technically it's not a wavelet transform
>>> but classified as a Gabor transform (special case of short-time FT).
>>>
>>>
>>> If I wrote something wrong, please correct me.
>>>
>>>
>>> Makoto
>>>
>>>
>>>
>>>
>>>
>>>
>>> On Fri, Aug 19, 2016 at 3:00 PM, Makoto Miyakoshi < mmiyakoshi at ucsd.edu
>>> > wrote:
>>>
>>>
>>>
>>> Dear Andreas, Niko, Norman, Mike, and the list,
>>>
>>> What a detailed discussion!
>>> Let me quickly confirm & summarize.
>>>
>>> > 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 ).
>>>
>>> I confirmed with your code that it's E ~= srate which is not normalized.
>>> Yes, if it's all the same, it's better to be normalized using one of these
>>> reasonable ways...
>>>
>>>
>>> I discussed this normalization issue with Rey in Feb 2007 in SCCN as a
>>> visitor (also asked the same question to Arno in August 2006 in email). It
>>> was one of the most puzzling problems that I found that Herrmann's Morlet
>>> equation was different from that of Tallon-Baudry's and Gruber's in the
>>> normalization factor A; to make things worse, EEGLAB's timef()'s dftfilt()
>>> was further different. Rey agreed with me, saying that he was puzzled too.
>>> I was glad to know that I was not the one who was confused! Several months
>>> before I met him, he had already written dftfilt3() to support the standard
>>> Morlet wavelet in EEGLAB, for which I've been grateful.
>>>
>>>
>>> Ok, let me confirm conclusions for the EEGLAB developers.
>>>
>>>
>>> 1) This is the bottom line: Let's not use the raw uV^2 outputs since
>>> these changes due to non-interesting factors as Mike mentioned. Instead,
>>> let's use dB changes against baseline etc.
>>>
>>>
>>> 2) If the raw uV^2 is absolutely necessary, there are probably two
>>> reasonable ways to normalize the raw values. One is, in Andreas's sample
>>> code, 'Analog to Gabor transform' in which A = 2.
>>>
>>>
>>> 3) The other good one is 'Unit energy' in which A =
>>> sqrt(pnts./(sqrt(pi)*sigma_f). Probably, EEGLAB developers can implement
>>> both as options and leave the choice to the users.
>>>
>>>
>>> 4) EEGLAB's explanation about dftfilt3() needs to be at least more
>>> detailed, and potentially the default/optional normalization can be changes
>>> to one of the above two, as long as it does not conflict with backword
>>> compatibility.
>>>
>>>
>>> Let me know if I'm wrong and/or missing something important.
>>>
>>>
>>> As a list member, I'd like to express my deep gratitude for your depth
>>> of knowledge and kindness to put so much effort on this discussion. I
>>> learned a lot!
>>>
>>>
>>> Makoto
>>>
>>>
>>>
>>>
>>> On Thu, Aug 18, 2016 at 10:49 AM, Andreas Widmann <
>>> widmann at uni-leipzig.de > wrote:
>>>
>>>
>>> P.S.: I forgot that I wanted to add to the final paragraph how to
>>> estimate the resulting amplitude of a sinusoid in the unit energy
>>> normalized wavelet transform: the energy of the wavelet with the
>>> corresponding frequency peak amplitude normalized to the RMS amplitude of
>>> the sinusoid, e.g. (for Niko’s example):
>>>
>>> mysig = 5 * sin(2*F*t*pi);
>>> [wavelet,cycles,freqresol,timeresol] = dftfilt3(F, ncycles, srate);
>>> temp = wavelet{1} / max( abs( wavelet{1} ) ) * sqrt( mean( abs( mysig )
>>> .^ 2 ) );
>>> sum( abs( temp ) .^ 2 )
>>> sqrt( sum( abs( temp ) .^ 2 ) )
>>>
>>> Best,
>>> Andreas
>>>
>>>
>>>
>>> > Am 17.08.2016 um 21:07 schrieb Norman Forschack < forschack at cbs.mpg.de
>>> >:
>>> >
>>> > Dear Andreas,
>>> >
>>> > thanks a lot for your elaborate input! Below some further
>>> clarifications.
>>> >
>>> > First of all, I did not intend to further discuss question (1) as I
>>> found it quite convincing from your previous post that dftfilt3 is not
>>> computing results as we expect it to do.
>>> > But probably I should have put my initial question into other words
>>> (e.g. 'One question was' instead of 'The initial question was'), in order
>>> not to confuse things. So I am sorry for the mess!
>>> >
>>> > Coming to question (2) I do not fully agree here. Literally, Nico's
>>> questions were: 'what do values in tf represent and how should I
>>> "normalize" these values to match the amplitude of the input signal?' So, I
>>> was aiming to give some input on the second part. Below, using your code,
>>> Nico's signal is plotted together with both types of amplitude envelopes,
>>> i.e. gabor normalization and unit energy (see also picture attached). It's
>>> clear that gabor transform (i.e. maximal wavelet spectrum power)
>>> normalization closely matches the amplitude of the original signal, while
>>> the unit energy envelope shows much higher values. Unless there's a major
>>> bug in the code below (please check), I'd say gabor is more appropriate in
>>> this scenario and the normalization of a complex wavelet spectrum to
>>> maximal power is much more straight forward (see below).
>>> >
>>> > However, I totally agree that when it comes to EEG, things are
>>> different. As Makoto and Mike already pointed out before, we usually do not
>>> interpret raw EEG power, but rather refer to a baseline or another
>>> condition to compare with. In this case it shouldn't matter which wavelet
>>> normalization has been applied.
>>> >
>>> > 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.
>>> >
>>> > 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?
>>> >
>>> > 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?)
>>> >
>>> > 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?
>>> > 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 am very much looking forward to your comments, thanks a lot Andreas!
>>> > Best
>>> > Norman
>>> >
>>> > figure links:
>>> > .fig: https://www.dropbox.com/s/ktoffojfcxxu6vd/wavelet_normalizat
>>> ion.fig?dl=0
>>> > .png: https://www.dropbox.com/s/v5i4d95y53t2woi/wavelet_normalizat
>>> ion.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.uc
>>> sd.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.uc
>>> sd.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
>>>
>>>
>>>
>>>
>>> --
>>>
>>>
>>> Makoto Miyakoshi
>>> Swartz Center for Computational Neuroscience
>>> Institute for Neural Computation, University of California San Diego
>>>
>>
>>
>>
>> --
>> Makoto Miyakoshi
>> Swartz Center for Computational Neuroscience
>> Institute for Neural Computation, University of California San Diego
>>
>> _______________________________________________
>> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
>> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
>> sd.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/20160825/df249c3a/attachment.html>