小波变换第五次作业专业：信息与通信工程 学号：406130714098 姓名：徐标1. 设计一 CQMFB ，低通滤波器 ()0H z 来自一半带滤波器。

该半带滤波器的长度为47，通带截止频率0.42p ωπ=，试给出()0H z ，()0G z ，()1H z ，()0z G 的幅频响应，单位抽样响应。

2.产生一信号()x n ，它由两个正弦加白噪声组成，一个在低频，一个在高频，正弦的频率及与白噪声的信噪比自己给定。

试用所设计的滤波器组对该信号进行分解和重建。

比较重建后的效果。

设计思路：参照课本181页（1） 首先设计一个半带滤波器()LF H z ，N=47，p ω=0.42π。

根据第六章半带滤波器的设计思路，先要用Chebyshew 最佳一致逼近法设计一个单带滤波器G （z ），令其通带截止频率为2p ω=0.84π，s ωπ=，长度为2J=24。

由此单带滤波器，可得半带滤波器()LF H z =()()1/2212N G z z --⎡⎤+⎣⎦，可以通过时域对g （n ）作二倍的插值，并令插值后的序列的中心点位0.5。

结果如下：（2）对半带滤波器()LF H z 进行处理，得到幅频响应非负的半带滤波器()P z 。

方法：令中间过度的滤波器：()()H jwjw LF LF eH e σ++=+，假定()LF H z 为零相位，实现上式的简单办法是令：()()()n 0,0LF LFLFh n h n h n n σ+≠⎧⎪=⎨+=⎪⎩，，再令()()0.5H 0.5LF P z z σ+=+，则()P z 是一个半带滤波器，()jw P e 是非负的。

产生半带滤波器()P z 的幅频响应：零极分析如图：由上面左图图可以看出，()P z 共有46对极零点，其中11个零点在单位圆内，11个零点在单位圆外，其余24个零点在单位圆上，对()P z 做谱分解，因为对于CQFMB 而言，谱分解的最佳选择是使()0H z 成为最小相位系统，取单位圆内的11个零点以及12个单位圆上的零点赋予()0H z ，从而构造出符合要求的()0H z ，并得到其幅频响应、单位抽样响应，同时，根据()()()()0101H z H z G z G z 之间的对应关系，很容易得到其他3个滤波器的幅频响应、单位抽样响应。

结果如下：()0jw H e 和()1jw H e 的对数幅频响应如下图 ：下面产生一个信号来检验刚才所设计的滤波器组的实际性能，产生一个信号()x n ，包含两个正弦信号，频率分别为119,30f Hz f Hz == ，还包含一个最大幅值 0.2 的高斯白噪。

计算得()x n 的信噪比 SNR=19.6040dB 。

抽样频率 200Hz 下，用刚才所设计滤波器组进行处理，得到信号如下程序：clear all;close all;N0=512;%First step: To design a one-band filter G(z) by Chebyshev approximation% set the cutoff frequency of the G(w);F=[0 0.84 1 1];A=[1 1 0 0];N=23;B=firpm(N,F,A)a=1;wf=0:pi/N0:pi*(N0-1)/N0;wff=0:1/N0:(N0-1)/N0;Gw=freqz(B,a,wf);Ew=exp(i*N*wf/2);Gr=real(Gw.*Ew);gn=impz(B,a,N);%plot the response curvefigure(1);subplot(2,1,1)plot(wff,real(Gr));grid on;title('single-band filt G(jw)'); subplot(2,1,2)stem((1:N),gn,'filled')title('single-band filt g(n)') hold onplot((0:N),zeros(N+1),'b');grid on;% To obtain half-band filter P(z) from G(z),% F(z)=[G(z^2)+(N-1)/2]/2N1=47;s=length(B);B2=zeros(1,2*s-1);for k=1:s,B2(k*2-1)=B(k);end,N2=length(B2);B2(N+1)=1+B2(N+1); B2=B2/2;Pw=freqz(B2,a,wf);Ew=exp(i*N*wf);Pr=real(Pw.*Ew);pn=impz(B2,a,N1);figure(2);subplot(2,1,1)plot(wff,Pr);grid on;title('half-band filt H(jw)'); subplot(2,1,2)stem(pn,'filled');grid on;title('half-band filt h(n) ') ; hold onplot((0:N1),zeros(N1+1),'b');grid on;%Modify P(z) and prepare to do spectrum decompositionminP=min(Pr);B2(N+1)=B2(N+1)+abs(minP);B3=B2*0.5/(0.5+abs(minP));Pw0=freqz(B3,a,wf);Ew=exp(i*N*wf);Pr0=real(Pw0.*Ew);figure;plot(wff,Pr0);grid;title('Impulse Response of P(z)');%Spectrum AnalysisA0=zeros(1,N1);A0(1)=1;figure(3);subplot(121);zplane(B3,A0);title('the zeros and poles before resolve, P(z)')[Z,P,K]=tf2zp(B3,A0);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.0001temp=real(temp);elsedisp('Spectrum Decomposition Failed!');endKK=sqrt(K/temp);l2=length(ZZ);PP=zeros(l2,1);[H0Num,H0Den]=zp2tf(ZZ,PP,KK); subplot(122)zplane(H0Num,H0Den);title('the zeros and poles after resolve, H_0(z)')H1Num=qmf(H0Num,1);G0Num=-wrev(H0Num);G1Num=qmf(G0Num);%H0H0w=freqz(H0Num,H0Den,wf);h0n=impz(H0Num,H0Den);figure(4);subplot(2,1,1)plot(wff,abs(H0w));grid on;title('Amplitude-Frequency Response of H_0(z)');subplot(2,1,2)stem(h0n,'filled');grid on;title('Impulse Response ofH_0(z)') ;hold on;plot((0:N),zeros(N+1),'b');grid on;%H1H1w=freqz(H1Num,H0Den,wf);%Ew=exp(i*N*wf);%H1r=real(H1w.*Ew);h1n=impz(H1Num,H0Den);figure(5);subplot(2,1,1) plot(wff,abs(H1w));grid on;title('Amplitude-Frequency Response of H_1(z)');subplot(2,1,2)stem(h1n,'filled');grid on;title('Impulse Response ofH_1(z)')hold on;plot((0:N),zeros(N+1),'b');grid on;%G0G0w=freqz(G0Num,H0Den,wf);g0n=impz(G0Num,H0Den);figure(6);subplot(2,1,1)plot(wff,abs(G0w));grid on;title('Amplitude-Frequency Response of G_0(z)');subplot(2,1,2)stem(g0n,'filled');grid on;title('Impulse Response ofG_0(z)')hold on;plot((0:N),zeros(N+1),'b');grid on;%G1G1w=freqz(G1Num,H0Den,wf);g1n=impz(G1Num,H0Den);figure(7);subplot(2,1,1)plot(wff,abs(G1w));grid on;title('Amplitude-Frequency Response of G_1(z)');subplot(2,1,2)stem(g1n,'filled');grid on;title('Impulse Response ofG_1(z)')hold on;plot((0:N),zeros(N+1),'b');grid on;%features of filter group[H0,w]= freqz(h0n,a);[H1,w]= freqz(h1n,a);absH0=abs(H0);absH1=abs(H1);ah0=20*log10(absH0);ah1=20*log10(absH1);figure(8)plot(wff,ah0,'k-',wff,ah1,'b-'); grid on;xlabel('w/pi');ylabel('|H_0(e^j^ w)|/dB,|H_1(e^j^w)|/dB');sumh=absH0.*absH0+absH1.*absH1; sum=10*log10(sumh);save H_filters h0ndisp('Press any key to process a signal using the filters...'); pause;clear all;close all;load H_filtersh0=h0n;h1=qmf(h0,1);g0=-wrev(h0);g1=qmf(g0);fs=200;f1=8;f2=30;step=-15;N=1000;n=1:N;Sig1=2*sin(2*pi*f1*n/fs);Sig2=2*sin(2*pi*f2*n/fs);P1=var(Sig1+Sig2)Nos=0.2*randn(N,1);P2=var(Nos);snr=10*log10(P1/P2) Sig=Sig1'+Sig2'+Nos;%signal v0(n) and v1(n) after decompositon and decimationx0=filter(h0,1,Sig);v0=decimate(x0,2); %low passx1=filter(h1,1,Sig);v1=decimate(x1,2); % high pass%reconstructed signal by filter [g0(n) g1(n) ]u0=interp(v0,2);u1=interp(v1,2);x_recon=filter(g0,1,u0)+filter(g 1,1,u1);figure(9);subplot(311)plot(Sig(100:200));grid;ylabel('x(n)');subplot(312)plot(x_recon(100+step:200+step));grid;ylabel('x(n).approx'); subplot(313)plot(Sig(100:200),'k');grid;hold onplot(x_recon(100+step:200+step),'b');ylabel('x(n),x(n).approx');。