1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
|
function cwt = CWT_Wavelab(x,nvoice,wavelet,oct,scale)
% CWT -- Continuous Wavelet Transform
% Usage
% cwt = CWT_Wavelab(x,nvoice,wavelet,oct,scale)
% Inputs
% x signal, dyadic length n=2^J, real-valued
% nvoice number of voices/octave
% wavelet string 'Gauss', 'DerGauss','Sombrero', 'Morlet'
% octave Default=2
% scale Default=4
% Outputs
% cwt matrix n by nscale where
% nscale = nvoice .* noctave
%
% Description
%
%
if nargin<4,
oct = 2;
scale = 4;
end
% preparation
x = ShapeAsRow(x);
n = length(x);
xhat = fft(x);
xi = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*pi/n);
% root
omega0 = 5;
% noctave = floor(log2(n))-2;
% noctave = floor(log2(n))-1;
noctave = floor(log2(n))-oct;
nscale = nvoice .* noctave;
cwt = zeros(n,nscale);
kscale = 1;
% scale = 4;
% scale = 16;
for jo = 1:noctave,
for jv = 1:nvoice,
qscale = scale .* (2^(jv/nvoice));
omega = n .* xi ./ qscale ;
if strcmp(wavelet,'Gauss'),
window = exp(-omega.^2 ./2);
elseif strcmp(wavelet,'DerGauss'),
window = i.*omega.*exp(-omega.^2 ./2);
elseif strcmp(wavelet,'Sombrero'),
window = (omega.^2) .* exp(-omega.^2 ./2);
elseif strcmp(wavelet,'Morlet'),
window = exp(-(omega - omega0).^2 ./2) - exp(-(omega.^2 + omega0.^2)/2);
end
% Renormalization
window = window ./ sqrt(qscale);
what = window .* xhat;
w = ifft(what);
cwt(1:n,kscale) = real(w)';
kscale = kscale+1;
end
scale = scale .*2;
end |
Partager