[Eeglablist] About CWT morlet wavelet and energy

[Miguel Valencia Ustárroz]

Dear all, I'd like to do time-frequency analyze with complex morlet 
wavelet,
up to now I was interested of the fase values this decomposition gave me,
but now I'd like to observe the energy of the signals I'm working with.
Following is my short matlab program, I think I get the correct 
coefficients,
but I'd like to transform the T-F coefficients  map in  a T-F energy map
related to my signal...


<---- CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    function [f,wt] = p_morlet(signal,Nh0,fmin,fmax,nfreqs);
%
%   Usage:
%        [f,wt] = p_morlet(signal,Nh0,fmin,fmax,nfreqs,lin_freq);
%
%   Inputs:
%        signal :         1-d signal
%        Nh0:            morlet wavelet parameter
%        fmin:            minimun analysis frequency, normalized  (0,0.5]
%        fmax:           maximum analysis frequency, normalized  (0,0.5]
%        nfreqs:         number of frequencies
%
%   Outputs:
%        f:                  analyzing frequencies
%        wt:               wavelet coefficients
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Setting up the analyzing frequencies
f=logspace(log10(fmin),log10(fmax),nfreqs);

%Init the Coefficients matrix
wt = [];

%For every analyzing frequency do the decomposition
for i=1:nfreqs,
    SIGMAt = Nh0/(2*pi*f(i));
    Tsigma = -fix(4*SIGMAt) : fix(4*SIGMAt);
   
    % Ordinary Morlet wavelet => it does not have zero mean depending on 
Nh0  (mainly if  Nh0 < 5)
    %psi = exp(-0.5*(Tsigma/SIGMAt).^2).*exp(j*2*pi*f(i).*Tsigma);
    %psi = psi./(sqrt(SIGMAt*sqrt(pi)));               % Normalizing  => 
energy of psi = 1 despite of the parameters

    % This version is independent to the parameters (for example we can 
use Nh0 = 2)
    psi = 
exp(-0.5*(Tsigma/SIGMAt).^2).*(exp(j*2*pi*f(i).*Tsigma)-exp(-0.5*(2*pi*f(i).*SIGMAt).^2));
    psi = 
psi./sqrt((SIGMAt*sqrt(pi)*(1+exp(-Nh0^2)-2*exp((-3*Nh0^2)/4))));   
%Normalization factor for the "new wavelet"
    
    %Convolving with the signal
    tmp = conv(signal,psi);

    %Saving the data of interest
     wt(i,:) = tmp(fix(4*SIGMAt)+1:length(tmp)-fix(4*SIGMAt));  
  
end
 >---- CODE

Up to I know, I can stimate the energy by getting (abs(wt))²; if I have 
the signal 'signal' mesured in Volts,
every coefficient wt(i,j) represents the energy of a region [dt x df], 
so I'd like to translate this value to  V²*s/Hz,
and I think this will depend on the sampling frequency I'm using and on 
the frequency band I'm analyzing.
Could someone help me how to scale the coeffients?
Thanks very much.
Best wishes...

       Miguel

-- 


Miguel Valencia Ustarroz
Fundacion para la Investigacion Medica Aplicada

Ed. Clinica Universitaria, Pio 12, n.36
Servicio de Neurofisiología Clínica
31080 Pamplona Spain
Tfn: +34 948 255400 -4595
Fax: +34 948 296500
e-mail: mvustarroz at unav.es
URL: http://www.unav.es/