%% 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');