[Eeglablist] Units in output of timefreq - wavelet normalization

[Norman Forschack]

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