[Eeglablist] Units in output of timefreq - wavelet normalization
Norman Forschack
forschack at cbs.mpg.de
Wed Aug 24 00:15:57 PDT 2016
Dear Makoto,
> But technically speaking, you vary the cycle > numbers and it's not wavelet anymore,
> correct?
Yes, I would say so. I just tried to outline that it's neither a Gabor transform/ STFT. For that, wavelets would need to have unit gain, fixed time and frequency resolution and require a different centering in time. So what is it then? We haven't talked about a label for this type. Is there any? Maybe just wavelet analysis with variable cycles numbers?
About 3, I agree it would be nice to have multiple options available once raw output normalization is needed and I also cannot see why unit energy normalization should be wrong. It's been out for a long while. But I'm not sure if it is applicable for a wavelet analysis with variable cycles. So far, Mike and Andreas were refering to scaled wavelets, i.e. with fixed number of cycles. Maybe they could comment on that.
Norman
---- Makoto Miyakoshi schrieb ----
>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_
>> 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
>>
>>
>>
>>
>> --
>>
>>
>> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20160824/5383601a/attachment.html>
More information about the eeglablist
mailing list