-
Notifications
You must be signed in to change notification settings - Fork 7
/
qmf_examples.m
114 lines (105 loc) · 3.54 KB
/
qmf_examples.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
clear all; close all; clc
%% QMF: critically decimated DFT filterbank with modulation seqs [1;1] and [1;-1]
fb = FilterBankStruct( );
fb.T = 2;
fb.B = 2;
Lh = 32;
Lg = Lh;
fb.symmetry = [1;0;0];
fb.tau0 = Lg-1; % We will have two QMF flavors: fb.i=1 for odd latency; fb.i=0 for even latency
eta = 1e4;
fb.w_cut = 0.5*pi;
lambda = 0;
best_cost = inf;
best_fb = fb;
for num_trial = 1 : 100
fb.h = rand(Lh,1); fb.g = rand(Lg,1);
[fb, cost, recon_err, iter] = FilterBankDesign(fb, eta, lambda, 100);
fprintf('Trial: %g; cost: %g; reconstruction error: %g; iterations %g\n', num_trial, cost, recon_err, iter)
if cost < best_cost
best_cost = cost;
best_fb = fb;
end
end
[fb, cost, recon_err, iter] = FilterBankDesign(best_fb, eta, lambda, 1000);
fprintf('Refinement. Cost: %g; reconstruction error: %g; iterations %g\n', cost, recon_err, iter)
figure;
subplot(3,1,1)
plot(conv(fb.h, fb.g));
subplot(3,1,2)
fft_size = 32768;
H = 20*log10(abs(fft(conv(fb.h, fb.g), fft_size)));
plot(pi*(0:fft_size/2-1)/(fft_size/2), H(1:end/2))
fb_qmf = fb;
%% We design a nonsubsampled version of the above wavelet
fb = FilterBankStruct( );
fb.T = 2;
fb.B = 1; % no decimation
Lh = 32;
Lg = Lh;
fb.tau0 = Lh - 1; % Also two flavors: fb.i=1 for odd latency; fb.i=0 for even latency
eta = 1e4;
fb.w_cut = 0.5*pi;
lambda = 0;
best_cost = inf;
best_fb = fb;
for num_trial = 1 : 100
fb.h = rand(Lh,1); fb.g = rand(Lg,1);
[fb, cost, recon_err, iter] = FilterBankDesign(fb, eta, lambda, 100);
fprintf('Trial: %g; cost: %g; reconstruction error: %g; iterations %g\n', num_trial, cost, recon_err, iter)
if cost < best_cost
best_cost = cost;
best_fb = fb;
end
end
[fb, cost, recon_err, iter] = FilterBankDesign(best_fb, eta, lambda, 1000);
fprintf('Refinement. Cost: %g; reconstruction error: %g; iterations %g\n', cost, recon_err, iter)
%figure;
subplot(3,1,1)
hold on; plot(conv(fb.h, fb.g));
title('(a) Time domain LP filter')
legend('Wavelet/QMF', 'Nonsubsampled version')
xlim('tight')
ylim('tight')
xlabel('Time')
ylabel('Impulse response')
subplot(3,1,2)
fft_size = 32768;
H = 20*log10(abs(fft(conv(fb.h, fb.g), fft_size)));
hold on; plot(pi*(0:fft_size/2-1)/(fft_size/2), H(1:end/2))
xlabel('\omega')
xlim('tight')
ylim([-60, 20])
ylabel('Magnitude in dB')
legend('Wavelet/QMF', 'Nonsubsampled version')
title('(b) Frequency domain LP filter')
fb_nonsubsampled = fb;
% Nearly perfect reconstruction (NPR) check
%
% Important: fb.i determines the phase of the modulation sequence of high pass analysis filters.
% If fb.i = 0, the modulation sequence of high pass analysis filters is
% [1, -1, 1, -1, ...]
% If fb.i = 1, the modulation sequence of high pass analysis filters is
% [-1, 1, -1, 1, ...]
%
% The same for QMF. We can design two flavors of QMFs
%
subplot(3,1,3)
fb = fb_nonsubsampled;
h_highpass = fb.h;
if fb.i==1
% the analysis modulation sequence is [-1, 1, -1, 1, ...]
h_highpass(1:2:end) = -h_highpass(1:2:end);
else % fb.i = 0
% the analysis modulation sequence is [1, -1, 1, -1, ...]
h_highpass(2:2:end) = -h_highpass(2:2:end);
end
% always have fb.j=0
% so the modulation sequence for high pass synthesis filter always is [1, -1, 1, -1, ...]
g_highpass = fb.g; g_highpass(2:2:end) = -g_highpass(2:2:end);
stem((conv(fb.h, fb.g) + conv(h_highpass, g_highpass))/2, '.')
xlim('tight')
ylim([-0.1, 1.1])
xlabel('Time')
ylabel('LP + HP impulse response', 'FontSize', 8)
title('(c) NPR check for the nonsubsampled version')