[Eeglablist] Fwd: Units in output of timefreq - wavelet normalization
Norman Forschack
forschack at cbs.mpg.de
Tue Aug 16 06:14:32 PDT 2016
Dear all,
below I forward a comment of Mike Cohen discussing wavelet normalization. He said, he will also upload a new lecturelet about the very topic on to his server probably by the end of the weekend.
Best
Norman
----- Forwarded Message -----
From: "Mike X Cohen" <mikexcohen at gmail.com>
To: forschack at cbs.mpg.de
Sent: Montag, 15. August 2016 19:43:01
Subject: Fwd: [Eeglablist] Units in output of timefreq - wavelet normalization
Dear Norman,
I'm not on the eeglab list, but someone just forwarded this message to me, so I thought I'd add a few more cents to the pot. Feel free to forward this message to the eeglab list.
In general, wavelet normalization in the time domain is an extremely difficult problem. It depends on many factors, including the frequency, the width of the Gaussian, the number of time points, and even the sampling rate. It is quite remarkable that everyone has gotten it wrong (mea culpa -- I didn't work through this issue in my book, so it's wrong in there as well). That factor of "A" that people often use comes from a normalization of the integral and simply doesn't translate to our sampled digital analysis environment. I think it hasn't been a huge problem because most people apply a baseline normalization (dB or percent change), so then the original scale doesn't matter.
Anyway, the solution is surprisingly simple: Normalize the wavelet energy in the frequency domain, not in the time domain. For wavelet convolution, it would work something like this:
cmwX = fft(cmw,nConv); % nConvolution-point FFT, assume 'cmw' is a complex Morlet wavelet
cmwX = cmwX./max(cmwX); % normalize to max-1 spectral energy
dataX = fft(data,nConv); % nConvolution-point FFT
convolution_result = ifft( dataX.*cmwX ); % convolution result
The power and real part of the convolution result now has the same units as the original signal (e.g., microvolts). If you want to have the time-domain Morlet wavelet, you could take the ifft of cmwX (using nfft according to the wavelet, not according to N+M-1 for convolution).
Hope that helps,
Mike
---------- Forwarded message ----------
From: Norman Forschack < forschack at cbs.mpg.de >
Date: Mon, Aug 15, 2016 at 6:19 PM
Subject: Re: [Eeglablist] Units in output of timefreq - wavelet normalization
To: Andreas Widmann < widmann at uni-leipzig.de >
Cc: eeglablist < eeglablist at sccn.ucsd.edu >
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
--
-----------------------------------------------------------------
Rasa Gulbinaite, PhD
Centre de Recherche Cerveau & Cognition (CerCo)
Toulouse (France)
e: rasa.gulbinaite at gmail.com
w: rasagulbinaite.com
--
Mike X Cohen, PhD
mikexcohen.com
More information about the eeglablist
mailing list