[Eeglablist] Units in output of timefreq - wavelet normalization

[Makoto Miyakoshi]

Dear Arno,

No, it's not easy for me. Primarily because I'm not a math person. But
another reason is that there seems at least a few good candidates to be a
default. After all, which is the best reasonable normalization do you think?

Makoto

On Thu, Sep 1, 2016 at 3:23 PM, Arnaud Delorme <arno at ucsd.edu> wrote:

> I am following on this discussion here and agree with Andreas
> recommendation.
> Andreas/Makoto, is it easy for you to send us a modified version of
> dftfilt3 with the proper scaling. We will name it dftfilt4 and, after
> testing, make it the default.
> For others, again, this should not modify any time-frequency plot where
> power is shown relative to baseline (which is the default).
>
> Thank you,
>
> Arno
>
> > On Aug 26, 2016, at 9:06 AM, Andreas Widmann <widmann at uni-leipzig.de>
> wrote:
> >
> > Dear all,
> >
> > for me the discussion has not yet come to a satisfying conclusion. I
> would like to argue that it is essential that we agree on a common
> terminology. I expect that some (normalization) issues will automagically
> disappear (but as always some new issues come up ;)
> >
> > (1) There are two families of transformations involved and intermixed in
> the discussion: Wavelet vs. STFT
> >
> > In wavelet transforms time-frequency decomposition is performed with
> template waveforms (the mother wavelet) which are scaled (dilated or
> compressed) and translated. Due to scaling the time and frequency
> resolution changes across the time-frequency-plane (TF-plane): high
> frequency and low time resolution at low frequencies and low frequency and
> high time resolution at high frequencies. The wavelet transform is a
> multiresolution analysis (MRA) and this is not a side effect but the core
> feature. There exist various wavelets with different properties; Morlet
> wavelets are sine based and thus closely related to STFT and FT. However,
> by scaling the number of cycles in a Morlet wavelet is never changed. I
> have never read a definition of (Morlet) wavelet transform in a textbook
> changing cycles instead of scaling a wavelet. If you find one, please, send
> me the reference.
> >
> > In contrast to the wavelet transform the STFT family has a constant
> fixed time and frequency resolution across the whole TF plane. There are
> various versions, in particular application of different window types but
> also variations of other parameters. As in Morlet wavelets a sinusoid is
> multiplied by a window function. To achieve constant frequency resolution
> of the atoms the number of cycles has to change with frequency (not
> arbitrarily but in a fixed relation as already annotated by Norman) while
> window width is kept constant.
> >
> > There are time-frequency transformations not clearly belonging to any of
> these two categories as "multiresolution STFT" and others. Neither wavelet
> transform nor STFT is to be preferred in general. They are equally valid
> but targeting different applications. MRA is preferred if frequency
> resolution is considered more important at low frequencies and time
> resolution at high frequencies as e.g. in audition. Googling „wavelet vs.
> STFT“ gives dozens of presentations and papers comparing both methods
> including their (dis-)advantages and applications.
> >
> > (2) Normalization has different requirements for wavelet transforms vs.
> STFT
> >
> > Scaling wavelets changes their energy. This is disadvantageous as now
> the energy of the output of the wavelet transform would change with scale.
> Therefore, most commonly sets of (scaled) wavelets are normalized to have
> all the same energy. It not necessarily matters whether this is unit energy
> or not. But personally I would prefer unit energy for two reasons. First,
> because it gives a straight-forward and intuitively accessible
> interpretation of the output amplitudes: the amplitude of a sample in the
> output reflects the integrated energy of the input in the surface area
> covered by the wavelet in the TF plane at the given scale and delay
> weighted by the window function (also see my previous message how to
> compute the wavelet transform amplitude of a sinusoid). Second, because
> data and amplitudes can be more easily compared across studies and methods
> (also see example code below).
> >
> > The atoms of the STFT, however, are already energy normalized "by
> design" as their window width is constant. Thus, energy normalization is
> not a core issue in STFT. Unit energy normalization might be advantageous
> for the reasons stated above but it will work also without. If one agrees
> to the above differentiation I would consider the fourth version in the
> code contributed by Mike ("Define empirical FWHM of wavelet in the
> frequency domain“) rather as a member of the STFT family (time and
> frequency resolution is constant). I have added a single set of these (once
> with and once without unit energy normalization) to my previous code
> hopefully demonstrating the differences between transforms and
> normalizations.
> >
> > (3) Suggestions for EEGLAB
> >
> > I would suggest to modify dftfilt3 so that by defaults wavelets are
> scaled to unit energy (as indirectly suggested by the documentation) and
> adding an option for backward compatible amplitude scaling. I would not add
> Gabor transform equivalent normalization to the multiresolution Morlet
> wavelet code (in case this is really needed this is a one-liner on the
> command line). Rather in my opinion a nice option would be to separately
> implement a textbook Gabor transform STFT in a separate function (actually
> reflecting Mike’s suggestion with minor changes; see Norman’s remarks).
> >
> > (4) General remark
> >
> > Sorry, I cannot withhold a more general „philosophical“ remark. Imho,
> the highest priority is replicability. Whatever analysis one performs
> should be described most thoroughly and detailed allowing for perfect
> replication and comparison of results across datasets, studies,
> researchers, labs, etc. This is considerably easier using textbook methods
> and terminology. For new methods and own approaches it is essential to
> describe them in a sufficiently detailed and accessible way. In my view
> this is the biggest problem with past and present literature applying TF
> transforms to electrophysiological data.
> >
> > Best,
> > Andreas
> >
> > fs = 256;
> > T = 4;
> > pnts = T * fs;
> > t = ( 0:pnts - 1 ) / fs;
> > cycles = 7;
> > F = 4:2:70;
> >
> > % White noise
> > rng( 0 )
> > signal = randn( fs * T, 1 );
> > % f = [ 10 50 ];
> > % signal = sin( 2 * pi * f( 1 ) * t )' + sin( 2 * pi * f( 2 ) * t )';
> >
> > %% 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_gabor = A .* cmorlf;
> >
> > tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* cmorlf_gabor );
> >
> > h = subplot( 4, 3, 3, 'YDir', 'normal' );
> > imagesc( t, F, abs( tf' ) .^ 2 )
> > colorbar
> > set( h, 'YDir', 'normal' );
> > subplot( 4, 3, 2 );
> > plot( F, mean( abs( tf ) .^ 2 ) )
> > subplot( 4, 3, 1 )
> > plot( (0:pnts / 2) * fs / pnts, cmorlf_gabor( 1:pnts / 2 + 1, : ) )
> >
> > % Unit energy
> > A = sqrt( pnts ./ ( sqrt( pi ) * sigma_f ) );
> > cmorlf_energy = A .* cmorlf;
> >
> > tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* cmorlf_energy );
> >
> > h = subplot( 4, 3, 6, 'YDir', 'normal' );
> > imagesc( t, F, abs( tf' ) .^ 2 )
> > colorbar
> > set( h, 'YDir', 'normal' );
> > subplot( 4, 3, 5 );
> > plot( F, mean( abs( tf ) .^ 2 ) )
> > subplot( 4, 3, 4 )
> > plot( (0:pnts / 2) * fs / pnts, cmorlf_energy( 1:pnts / 2 + 1, : ) )
> >
> > %% STFT
> > % adapted from https://sccn.ucsd.edu/pipermail/eeglablist/2016/
> 011536.html
> > % by mikexcohen at gmail.com
> > hz = (0:pnts - 1) * fs / pnts;
> > fwhm = 10;
> > for iF = 1:length( F )
> >    s  = fwhm*(2*pi-1)/(4*pi);       % normalized width
> >    x  = hz-F( iF );                 % shifted frequencies
> >    fx = exp(-.5*(x/s).^2)';         % gaussian
> >    fxArray( :, iF )  = fx./max(fx); % gain-normalize
> > end
> >
> > tf = sqrt( 2 ) * ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .*
> fxArray );
> >
> > h = subplot( 4, 3, 9, 'YDir', 'normal' );
> > imagesc( t, F, abs( tf' ) .^ 2 )
> > colorbar
> > set( h, 'YDir', 'normal' );
> > subplot( 4, 3, 8 );
> > plot( F, mean( abs( tf ) .^ 2 ) )
> > subplot( 4, 3, 7 )
> > plot( (0:pnts / 2) * fs / pnts, fxArray( 1:pnts / 2 + 1, : ) )
> >
> > % Unit energy
> > fxArray_unit = fxArray * sqrt( pnts / max( sum( fxArray .^ 2 ) ) );
> >
> > tf = ifft( repmat( fft( signal ), [ 1 length( F ) ] ) .* fxArray_unit );
> >
> > h = subplot( 4, 3, 12, 'YDir', 'normal' );
> > imagesc( t, F, abs( tf' ) .^ 2 )
> > colorbar
> > set( h, 'YDir', 'normal' );
> > subplot( 4, 3, 11 );
> > plot( F, mean( abs( tf ) .^ 2 ) )
> > subplot( 4, 3, 10 )
> > plot( (0:pnts / 2) * fs / pnts, fxArray_unit( 1:pnts / 2 + 1, : ) )
> >
> >> Am 26.08.2016 um 03:11 schrieb Makoto Miyakoshi <mmiyakoshi at ucsd.edu>:
> >>
> >> Dear Norman,
> >>
> >>> 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?
> >>
> >> MJWAWVCN method? That's certainly too long and complicated. It
> describes what it does, but as a definition it's self-contradictory and
> bad. I wonder what Andreas would say on this.
> >>
> >>> But I'm not sure if it is applicable for a wavelet analysis with
> variable cycles.
> >>
> >> Normalization is a normalization, it does not matter whether the cycles
> are fixed or variable. Am I wrong?
> >>
> >> Makoto
> >>
> >>
> >> On Wed, Aug 24, 2016 at 12:15 AM, Norman Forschack <
> forschack at cbs.mpg.de>wrote:
> >> 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.dewrote:
> >>>>
> >>>>> 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
> >>
> >>
> >>
> >> --
> >> 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.
> 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/20160909/fcbb7fd1/attachment.html>