%---------------------------------------------------------------------------------------
% exa010201b.m,  for  example 1.2.3
% 說明信號的子帶分解；
%----------------------------------------------------------------------------------------
clear;

%以下內容涉及到涉及對信號分解的半帶濾波器，具體方法見本書第7章；
%-------------------------------------------------------------------------------
N=25;
Wp=0.46
F=[0 2*Wp 1 1];
A=[1 1 0 0];
W=[1 1];
G=REMEZ(N,F,A,W);
for j=1:N+1
   tmp(2*j-1)=0.5*G(j);
   tmp(2*j)=0;
end
P=tmp(1:(2*N+1));
P(N+1)=P(N+1)+0.5;
[H,w1]=freqz(P,1);
err=max(abs(H))-1
P(N+1)=P(N+1)+err;

P=P*0.5/(0.5+err);
[h0,P1]=resolve(P)
h1=qmf(h0,1)
%-------------------------------------------------------------------------------

%產生實驗信號x（n）；
n=256;
w=pi/5;
pow=0.08;
u=randn(1,n);
umean=mean(u);
u=u-umean;
u_p=dot(u,u)/n;
u=u.*sqrt(pow/u_p);
   for j=1:n
      x(j)=cos(0.8*w*j)+cos(1.5*w*j);
   end;
   x=2*x+u;
   
% 信號分解，并顯示分解后的信號；
x0=filter(h0,1,x);
x1=filter(h1,1,x);
figure(1)
subplot(221)
plot(x0);grid on;
ylabel(' x0(n)');
subplot(222)
plot(x1);grid on;
ylabel(' x1(n)');

% 顯示分解后信號的頻譜；
w=0:1/n:0.5-1/n;
subplot(223)
x00=fft(x0);
x00=20*log10(abs(x00));
plot(w,x00(1:n/2));grid on;
ylabel(' the Spectrum of x0(dB)');
subplot(224)
x11=fft(x1);
x11=20*log10(abs(x11));
plot(w,x11(1:n/2));grid on;
ylabel(' the Spectrum of x1(dB)');