www.pudn.com > tftb2002toolbox.rar > gaussn.m

function g = gaussn(f0,n) 
 
% function gn = gaussn(f0,n) : generates the order n derivative of the  
% gaussian window, centered at frequency f0 
% The wavelet gn is real, but it is its analytic form that is  
% synthetized first. The real part is then normalized by the analytical 
% value of its energy (so, L2 normalization!) 
 
% alpha : std of the initial Gaussian  
% g_0(t) = sqrt(2) (alpha/2*pi)^(1/4) exp(-alpha t^2) 
 
alpha = 2*pi^2*f0^2 / n ; 
 
% N = number of points, depends on "n". The support of the wavelet  
% at tolerance 10^(-3) is determined experimentaly and modeled by 
% a 4th order polynome (with reference frequency f0 is 0.1)  
 
PolyCoef = [-0.00000001420054   0.00001199042535  -0.00403466080616 ... 
      0.92639865727446 20.64720217778273] ; 
N = ceil((0.1/f0).*polyval(PolyCoef,n)) ; 
 
f = linspace(0,1,2*N+1) ; f = f(1:2*N) ; 
 
% Synthesis of the "analytic" wavelet Gn in the frequency domain  
 
G = (2*pi/alpha)^(1/4) * (2*i*pi*f).^n .* exp(-(pi^2*f.^2)./alpha) ; 
 
% Calculus of the wavelet energy (theoretical expression, cf. Gratshteyn, 
% sec. 3.461, form. 2) 
 
p = (2*pi^2)/alpha ; 
E = (2*pi/alpha)^(1/2) * (2*pi).^(2*n) * ... 
    (fact2(2*n-1))/(2*(2*p)^n)*sqrt(pi/p) ; 
if isinf(E) | isnan(E) | E == 0 
  E = integ(abs(G).^2,f) ; 
end 
 
% L2 Normalized wavelet in time domain  
 
g = real ( fftshift(ifft(G))./(sqrt(E/2)) ) ; 
g = g(2:2*N) ; 
t = -(N-1):(N-1) ; 
 
% subplot(211) ; 
% plot(f,abs(G)) ; title ('Wavelet Spectrum') ; 
% xlabel('Frequency') ; grid  
% subplot(212) 
% plot(t,g) 
% title('Wavelet in time') 
% xlabel('Time') ; grid