%% Example 4.18 Daubechies 9/7 filters
%
% I follow the construction of these filters given in Section 4.12 #11.
%

syms z y w
B = 1+4*y+10*y^2+20*y^3

%
% I find the roots of B.

r = roots([20 10 4 1]')

%
% I see the real root is the third one. The first two are a conjugate pair.

p1 = y-r(3)
p2 = (y-r(1))*(y-r(2))

%

H0 = sqrt(2)*z^(-2)*(1+z)^4*subs(p2,y,(1-z)^2/(-4*z))/(16*subs(p2,y,0));

%

G0 = sqrt(2)*(1+z)^4*subs(p1,y,(1-z)^2/(-4*z))/(16*subs(p1,y,0));

%

H1 = z*subs(G0,z,-z);
G1 = subs(H0,z,-z)/z;

%
% I substitute z=e^{i\omega } into the absolute values of H_0 \text{ },H_1 \text{ },G_0 \text{ },G_1.
%

absH0 = abs(subs(H0,z,exp(i*w)))
absH1 = abs(subs(H1,z,exp(i*w)))
absG0 = abs(subs(G0,z,exp(i*w)))
absG1 = abs(subs(G1,z,exp(i*w)))

%
% Here are the graphs of the frequency response of the analysis filters.
%

figure
fplot(absH0,[0,pi])
hold on
fplot(absH1, [0 pi])
axis([0 pi 0 1.6])
xticks([0 pi/4 pi/2 3*pi/4 pi])
xticklabels({'0','.25','.5', '.75', '1'})
title('Frequency response of Daub 9/7 analysis filters')
xlabel('frequency \omega (multiples of \pi)')
annotation('textarrow',[0.31 0.38],...
    [0.71 0.8],'String','low-pass filter');
annotation('textarrow',[0.7 0.75],...
    [0.73 0.85],'String','high-pass filter');
annotation('textarrow',[0.45 0.53],...
    [0.62 0.63],'String','crossover point');

%
% Here are the graphs for the synthesis filters.
%

figure
fplot(absG0,[0,pi])
hold on
fplot(absG1, [0 pi])
axis([0 pi 0 1.6])
xticks([0 pi/4 pi/2 3*pi/4 pi])
xticklabels({'0','.25','.5', '.75', '1'})
title('Frequency response of Daub 9/7 synthesis filters')
xlabel('frequency \omega (multiples of \pi)')
annotation('textarrow',[0.31 0.38],...
    [0.71 0.8],'String','low-pass filter');
annotation('textarrow',[0.67 0.75],...
    [0.67 0.8],'String','high-pass filter');
annotation('textarrow',[0.37 0.48],...
    [0.57 0.63],'String','crossover point');