%-----------------------------------------------------------------------------------
%  exa070501.m    for example 7.5.1
%  To explain how to design a CQMFB;
%-----------------------------------------------------------------------------------
clear;

% 先設計一個單帶濾波器；
N0=512; N=21;Fs=1;
Wp=0.45
F=[0 2*Wp 1 1];
A=[1 1 0 0];
W=[1 1];
a=1;
hsb=remez(N,F,A,W);
[Hsb,w]=freqz(hsb,a,N0,Fs);
figure(1)
subplot(221)
plot(w,abs(Hsb));grid;
%title('single band filter hsb(n)')


% 由單帶濾波器得到半帶濾波器
for j=1:N+1
   tmp(2*j-1)=0.5*hsb(j);
   tmp(2*j)=0;
end
hlf=tmp(1:(2*N+1));
hlf(N+1)=hlf(N+1)+0.5;

wf=0:pi/N0:pi*(N0-1)/N0;
wff=0:0.5/N0:(N0-1)/N0/2;
a=1;
Gw= freqz(hlf,a,wf);
Ew=exp(i*2*N*wf/2);
Gr=Ew.*Gw;
% Gr 現在是零相位半帶濾波器的頻率特性，，其虛部應為零，實際上是濾波器的“增益”
mmax=max(real(Gr))
mmin=min(real(Gr))

subplot(222)
plot(wff,real(Gr));grid;
%title('half band filter HLF(z)')

% 獲得“增益”恒正的半帶濾波器， P+
hlf(N+1)=hlf(N+1)+abs(mmin);
p=hlf*0.5/(0.5+abs(mmin));
subplot(223)
stem(p,'.');grid;
%title('halfband filter p(n)');

Gw= freqz(p,a,wf);
Gr=Ew.*Gw;
subplot(224)
plot(wff,real(Gr));grid;
%title('halfband filter P(z)');
%save p.mat p;
%---------------------------------------------------------------------
% 譜分解部分，求出 h0(n);
l=2*N+1;
A=zeros(l,1);
A(1)=1;
A=A';

figure(2);
subplot(221);
zplane(p,A);
title('the zeros and poles before resolve') 

[Z,P,K]=tf2zp(p,A);
Z1=sort(Z)
ll=length(Z1)/2;
ZZ=Z1(1:ll);
temp=1;
for m=1:ll
    temp=(-ZZ(m))*temp;
end
if imag(temp)<0.0001
   temp=real(temp);
end
KK=sqrt(K/temp);
l2=length(ZZ);
PP=zeros(l2,1);
[h0,P1]=zp2tf(ZZ,PP,KK);
subplot(222)
zplane(h0,P1);
title('the zeros and poles after resolve') 
%-----------------------------------------------------------------------------
% 利用CQMFB四個濾波器之間的關系，得到 h1, g0, g1;
h1=qmf(h0,1);
g0=-wrev(h0);
g1=qmf(g0);

figure(3)
subplot(221)
stem(h0,'.');grid;
title(' filter h0(n)')
subplot(223)
stem(h1,'.');grid;
title(' filter h1(n)')
subplot(222)
stem(g0,'.');grid;
title(' filter g0(n)')
subplot(224)
stem(g1,'.');grid;
title(' filter g1(n)')

% 求出并畫出分析濾波器組的對數幅頻特性；
[H0,w]= freqz(h0,a,N0,Fs);
[H1,w]= freqz(h1,a,N0,Fs);
absH0=abs(H0);
absH1=abs(H1);
ah0=20*log10(absH0);
ah1=20*log10(absH1);
figure(4)
subplot(221)
plot(w,ah0,'k-',w,ah1,'b-');
grid on;

% 檢測兩個分析濾波器的功率互補情況；
sumh=absH0.*absH0+absH1.*absH1;
sum=10*log10(sumh);
subplot(222)
plot(w,sum);grid;

save h0.mat h0