[Eeglablist] Units in output of timefreq - wavelet normalization

[Mike X Cohen]

Dear Abdul,

Thanks for that link. It's a nice method. However, I'm not sure that an
approach like this would be optimal for most applications in EEG. You are
correct that EEG contains nonstationarities. However, there is also
variability across trials within an individual, and variability across
individuals. So in fact, highly precise results can be suboptimal. For
example, if trial 1 contains a burst at 5 Hz 250 ms while trial 2 contains
a burst at 6 Hz at 380 ms, then such a precise TF localization would
produce attenuated or no results in trial-averaging. Therefore, some
smoothing is actually beneficial and sometimes necessary. In fact, this
smoothing is the main motivation for applying multitaper analyses to data
above, e.g., 30 Hz.

Nonetheless, there are situations, like SSVEP, where very narrow-band and
precise TF localization is important, and I think a method like what you
propose could be useful for these specific situations.

Mike




On Fri, Aug 26, 2016 at 7:05 AM, Md. Abdul Awal <awalece04 at yahoo.com> wrote:

> Related to the discussion, I would like to add:
>
> As you know, EEG is a non-stationary signal i.e, its spectral contents are
> timevarying. So, joint time-frequency analyis is very useful.
> Joint time-frequency  anlaysis is done by STFT and Wavelet and using
> time-frequency distribution.
>
> STFT and wavelet provide low resolution. One the hand, time-frequency
> distribution offers high resolution.
>
> Another serious problem as discussed in the email bellow that, one should
> never blindly trust a single set of parameters….
>
> This is an age-old problem.  In our recent paper, we solved this problem.
> Users only need the signal itself. The algorithm will automatically provide
> optimal parameters and provide high resolution time-frequency resolution.
>
> I think this could be very useful in this context.
>
> Here is the paper: M. A. Awal and B. Boashash, "An automatic fast
> optimization of Quadratic Time-frequency Distribution using the hybrid
> genetic algorithm," Signal Processing, vol. 131, pp. 134-142, February,
> 2017.
> Link:  http://authors.elsevier.com/a/1TbLWbZX4o2FF
> One can now freely download the full paper: *valid for 50 days, until
> October 12, 2016.*
>
> with regards,
> ------------------------
> Md. Abdul Awal
> The University of Queensland, Australia
>
>
> ------------------------------
> *From:* Makoto Miyakoshi <mmiyakoshi at ucsd.edu>
> *To:* Mike X Cohen <mikexcohen at gmail.com>
> *Cc:* eeglablist <eeglablist at sccn.ucsd.edu>; Andreas Widmann <
> widmann at uni-leipzig.de>; "ramon at sccn.ucsd.edu" <ramon at sccn.ucsd.edu>
> *Sent:* Friday, August 26, 2016 10:57 AM
>
> *Subject:* Re: [Eeglablist] Units in output of timefreq - wavelet
> normalization
>
> Dear Mike,
>
> > Perhaps the more accurate description would be a "set of Morlet
> wavelets" or a "bank of Morlet wavelets," which simply states that there
> are sine waves being tapered by Gaussians.
>
> Yes I agree.
>
> > At the risk of belaboring this point, one should never blindly trust a
> single set of parameters just because it's the default or because it's what
> someone wrote in a paper. This is particularly true for filtering, which is
> really what convolution and wavelets are all about. Always try a range of
> parameters and different analysis approaches to make sure you have a real
> finding, not some bizarre quirk of one specific parameter set.
>
> That's good general advice. As long as they can afford it, yes. On the
> other hand, there are psychologists (me) and medical doctors in EEG users
> who are not fluent in math or coding. For those users, I hope EEGLAB, for
> example, provides the most established and reasonable default choices and
> parameters.
>
> > I suppose the question of what to do in EEGLAB is about finding the
> right balance between what most users would want, while providing a few
> advanced options without overwhelming users.
>
> For that, user feedback is indispensable... hence I promote your post to
> EEGLAB Bugzilla!
>
> Thank you Mike.
>
> Makoto
>
>
>
> On Wed, Aug 24, 2016 at 12:14 AM, Mike X Cohen <mikexcohen at gmail.com>
> wrote:
>
> Hi all. "Morlet wavelet" or "Gabor" is mostly just nomenclature. Perhaps
> the more accurate description would be a "set of Morlet wavelets" or a
> "bank of Morlet wavelets," which simply states that there are sine waves
> being tapered by Gaussians. To speak of a mother wavelet and daughter
> wavelets would imply that they have identical parameters but different sine
> waves.
>
> More generally, there are many right ways to analyze data. There are, of
> course, many wrong ways as well, but if you have a critical eye, the wrong
> ways are usually easy to spot. Different normalizations or parameter
> settings can be more or less appropriate in certain situations. And
> normalization is always a tricky issue, not only in time-frequency analyses
> (and not only in neuroscience). At the risk of belaboring this point, one
> should never blindly trust a single set of parameters just because it's the
> default or because it's what someone wrote in a paper. This is particularly
> true for filtering, which is really what convolution and wavelets are all
> about. Always try a range of parameters and different analysis approaches
> to make sure you have a real finding, not some bizarre quirk of one
> specific parameter set.
>
> I suppose the question of what to do in EEGLAB is about finding the right
> balance between what most users would want, while providing a few advanced
> options without overwhelming users. And if someone really needs a
> custom-designed analysis, they can write their own code. The math and
> programming behind these analyses is generally not difficult (<apologies
> for shameless promotion> Cohen, 2014, MIT Press).
>
> Mike
>
>
> On Wed, Aug 24, 2016 at 2: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,time resol] = 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,time resol] = 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*p i*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/ktof fojfcxxu6vd/wavelet_normalizat
> ion.fig?dl=0
> <https://www.dropbox.com/s/ktoffojfcxxu6vd/wavelet_normalization.fig?dl=0>
> > .png: https://www.dropbox.com/s/v5i4 d95y53t2woi/wavelet_normalizat
> ion.png?dl=0
> <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*p i*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(ha lfMsiz));
> >>>> 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*p i*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,time resol] = 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,time resol] = 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)./n ormfactor)
> >>>>>> xlabel('Time (seconds)');
> >>>>>> ylabel('Amplitude');
> >>>>>> legend('input signal', 'wavelet result')
> >>>>>>
> >>>>>> ______________________________ _________________
> >>>>>> Eeglablist page: http://sccn.ucsd.edu/eeglab/ee glabmail.html
> <http://sccn.ucsd.edu/eeglab/eeglabmail.html>
> >>>>>> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
> sd.edu <eeglablist-unsubscribe at sccn.ucsd.edu>
> >>>>>> For digest mode, send an email with the subject "set digest mime" to
> >>>>>> eeglablist-request at sccn.ucsd.e du
> <eeglablist-request at sccn.ucsd.edu>
> >>>>>
> >>>>> ______________________________ _________________
> >>>>> Eeglablist page: http://sccn.ucsd.edu/eeglab/ee glabmail.html
> <http://sccn.ucsd.edu/eeglab/eeglabmail.html>
> >>>>> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
> sd.edu <eeglablist-unsubscribe at sccn.ucsd.edu>
> >>>>> For digest mode, send an email with the subject "set digest mime" to
> >>>>> eeglablist-request at sccn.ucsd.e du <eeglablist-request at sccn.ucsd.edu>
> >>>> ______________________________ _________________
> >>>> Eeglablist page: http://sccn.ucsd.edu/eeglab/ee glabmail.html
> <http://sccn.ucsd.edu/eeglab/eeglabmail.html>
> >>>> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
> sd.edu <eeglablist-unsubscribe at sccn.ucsd.edu>
> >>>> For digest mode, send an email with the subject "set digest mime" to
> >>>> eeglablist-request at sccn.ucsd.e du <eeglablist-request at sccn.ucsd.edu>
> > ______________________________ _________________
> > Eeglablist page: http://sccn.ucsd.edu/eeglab/ee glabmail.html
> <http://sccn.ucsd.edu/eeglab/eeglabmail.html>
> > To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
> sd.edu <eeglablist-unsubscribe at sccn.ucsd.edu>
> > For digest mode, send an email with the subject "set digest mime" to eeglablist-request at sccn.ucsd.e
> du <eeglablist-request at sccn.ucsd.edu>
>
> ______________________________ _________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/ee glabmail.html
> <http://sccn.ucsd.edu/eeglab/eeglabmail.html>
> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
> sd.edu <eeglablist-unsubscribe at sccn.ucsd.edu>
> For digest mode, send an email with the subject "set digest mime" to eeglablist-request at sccn.ucsd.e
> du <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
>
>
>
>
> --
> Mike X Cohen, PhD
> mikexcohen.com
>
>
>
>
> --
> 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
>
>
>


-- 
Mike X Cohen, PhD
mikexcohen.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20160826/da7365b8/attachment.html>