/
wfiltfn.m
132 lines (125 loc) · 5.11 KB
/
wfiltfn.m
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
% function psihfn = wfiltfn(type, opt)
%
% Wavelet transform function of the wavelet filter in question,
% fourier domain.
%
% Input:
% type: string (see below)
% opt: options structure, e.g. struct('s',1/6,'mu',2)
%
% Output:
% anonymous function @(xi) psihfn(xi)
%
% Filter types Use for synsq? Parameters (default)
%
% gauss yes s (1/6), mu (2)
% mhat no s (1)
% cmhat yes s (1), mu (1)
% morlet yes mu (2*pi)
% shannon no -- (NOT recommended for analysis)
% hshannon yes -- (NOT recommended for analysis)
% hhat no
% hhhat yes mu (5)
% bump yes s (1), mu (5)
%
% Example:
% psihfn = wfiltfn('bump', struct('mu',1,'s',.5));
% plot(psihfn(-5:.01:5));
%
%---------------------------------------------------------------------------------
% Synchrosqueezing Toolbox
% Authors: Eugene Brevdo, Gaurav Thakur
%---------------------------------------------------------------------------------
function psihfn = wfiltfn(type, opt, derivative)
if nargin<3, derivative=false; end
switch type
%bump window should be optimal for most purposes
case 'bump',
if ~isfield(opt,'mu'), mu = 1; else mu = opt.mu; end
if ~isfield(opt,'s'), s = 1; else s = opt.s; end
if ~isfield(opt,'om'), om = 0; else om = opt.om; end
psihfnorig = @(w) (abs(w)<0.999).*exp(-1./(1-(w.*(abs(w)<0.999)).^2))/0.443993816053287;
if ~(derivative)
psihfn = @(w) exp(2*pi*i*om*w).*psihfnorig((w-mu)/s)/s;
else
psihfn = @(w) exp(2*pi*i*om*w).*psihfnorig((w-mu)/s)/s .* (2*pi*i*om - 2*((w-mu)/s^2)./(1-((w-mu)/s).^2).^2);
end
case 'gauss',
% Similar to morlet, but can control bandwidth
% can be used with synsq for large enough mu/s ratio
if ~isfield(opt,'s'), s = 1/pi; else s = opt.s; end
if ~isfield(opt,'mu'), mu = 0; else mu = opt.mu; end
if ~(derivative)
psihfn = @(w) (2*pi*s^2)^(-1/2)*exp(-(w-mu).^2/(2*s^2));
else
psihfn = @(w) (2*pi*s^2)^(-1/2).*exp(-(w-mu).^2/(2*s^2)).*-(w-mu)/s^2;
end
case 'mhat', % mexican hat
if ~isfield(opt,'s'), s = 1; else s = opt.s; end
psihfn = @(w) -sqrt(8)*s^(5/2)*pi^(1/4)/sqrt(3)*w.^2.*exp(-s^2*w.^2/2);
case 'cmhat',
% complex mexican hat: hilbert analytic function of sombrero
% can be used with synsq
if ~isfield(opt,'s'), s = 1; else s = opt.s; end
if ~isfield(opt,'mu'), mu = 1; else mu = opt.mu; end
psihfnshift = @(w)2*sqrt(2/3)*pi^(-1/4)*s^(5/2)*w.^2.*exp(-s^2*w.^2/2).*(w>=0);
psihfn = @(w) psihfnshift(w-mu);
case 'morlet',
% can be used with synsq for large enough s (e.g. >5)
if ~isfield(opt,'mu'), mu = 2*pi; else mu = opt.mu; end
cs = (1+exp(-mu^2)-2*exp(-3/4*mu^2)).^(-1/2);
ks = exp(-1/2*mu^2);
psihfn = @(w)cs*pi^(-1/4)*(exp(-1/2*(mu-w).^2)-ks*exp(-1/2*w.^2));
case 'shannon',
psihfn = @(w)exp(-i*w/2).*(abs(w)>=pi & abs(w)<=2*pi);
case 'hshannon',
% hilbert analytic function of shannon transform
% time decay is too slow to be of any use in synsq transform
if ~isfield(opt,'mu'), mu = 0; else mu = opt.mu; end
psihfnshift = @(w)exp(-i*w/2).*(w>=pi & w<=2*pi).*(1+sign(w));
psihfn = @(w) psihfnshift(w-mu);
case 'hhat', % hermitian hat
psihfn = @(w)2/sqrt(5)*pi^(-1/4)*w.*(1+w).*exp(-1/2*w.^2);
case 'hhhat',
% hilbert analytic function of hermitian hat
% can be used with synsq
if ~isfield(opt,'mu'), mu = 5; else mu = opt.mu; end
psihfnshift = @(w)2/sqrt(5)*pi^(-1/4)*(w.*(1+w).*exp(-1/2*w.^2)) .* (1+sign(w));
psihfn = @(w) psihfnshift(w-mu);
otherwise
error('Unknown wavelet type: %s', type);
end % switch type
end
% case 'daub',
% %
% % Construction of such filters follows from Mallat
% % Eq. (7.140), where g is given by Mallat Eqs. (7.53) and
% % (7.68) . Note that these are DTFT equations.
% %
% if nargin<4,
% a = 2; % Dyadic scales
% p = 2; % Daubechies-2 wavelets (2 vanishing moments)
% j = 2*(ceil(log(N)/log(2))-1); % Which scale?
% else
% a = opt.a;
% p = opt.p;
% j = opt.j;
% end
% L = ceil(log(N)/log(2));
% assert(j>=L);
% [h,g,rh,rg] = daub(2*p);
% hmat = repmat(h(end:-1:1)', 1, N);
% % h^(w) and g^(w)
% hh = @(wp) sum(hmat .* exp(-i*[0:2*p-1]'*wp));
% gh = @(wp) exp(-i*wp) .* conj(hh(wp+pi));
% % Calculate g^(a^{j-L-1} w) prod_{k=0}^{j-L-2} h^(a^k w)
% psih = gh(a^(j-L-1) * w);
% for k=0:(j-L-2), psih = psih .* hh(a^k * w); end
% % TODO -- why is this necessary?
% %psih = -psih;
% return;
% case 'symm', % TODO - make a wfilt_qmf for this
% if nargin<4, p = 2; else p = opt.p; end
% [h,g,rh,rg] = symmlets(p);
% psih = zeros(size(W));
% for pn=1:p, psih = psih + rg(pn)*exp(-i*(pn-1)*W); end