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 |
