小波变换第五次作业.doc

小波变换第五次作业专业：信息与通信工程 学号：406130714098 姓名：徐标1. 设计一 CQMFB，低通滤波器 来自一半带滤波器。该半带滤波器的长度为47，通带截止频率，试给出，的幅频响应，单位抽样响应。2. 产生一信号，它由两个正弦加白噪声组成，一个在低频，一个在高频，正弦的频率及与白噪声的信噪比自己给定。试用所设计的滤波器组对该信号进行分解和重建。比较重建后的效果。设计思路：参照课本181页（1） 首先设计一个半带滤波器，N=47，=0.42。根据第六章半带滤波器的设计思路，先要用Chebyshew最佳一致逼近法设计一个单带滤波器G（z），令其通带截止频率为2=0.84，长度为2J=24。由此单带滤波器，可得半带滤波器=，可以通过时域对g（n）作二倍的插值，并令插值后的序列的中心点位0.5。结果如下：（2）对半带滤波器进行处理，得到幅频响应非负的半带滤波器。 方法：令中间过度的滤波器：，假定为零相位，实现上式的简单办法是令：，再令，则是一个半带滤波器，是非负的。产生半带滤波器的幅频响应： 零极分析如图：由上面左图图可以看出，共有46对极零点，其中11个零点在单位圆内，11个零点在单位圆外，其余24个零点在单位圆上，对 做谱分解，因为对于CQFMB而言，谱分解的最佳选择是使成为最小相位系统，取单位圆内的11个零点以及12个单位圆上的零点赋予 ，从而构造出符合要求的，并得到其幅频响应、单位抽样响应，同时，根据 之间的对应关系，很容易得到其他3个滤波器的幅频响应、单位抽样响应。结果如下：和的对数幅频响应如下图 ：下面产生一个信号来检验刚才所设计的滤波器组的实际性能，产生一个信号 ，包含两个正弦信号，频率分别为 ，还包含一个最大幅值 0.2 的高斯白噪。计算得 的信噪比 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 curve figure(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 on plot(0:N),zeros(N+1),b);grid on; % To obtain half-band filter P(z) from G(z), % F(z)=G(z2)+(N-1)/2/2 N1=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 on plot(0:N1),zeros(N1+1),b);grid on; %Modify P(z) and prepare to do spectrum decomposition minP=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 Analysis A0=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.0001 temp=real(temp); else disp(Spectrum Decomposition Failed!); end KK=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); %H0 H0w=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 of H_0(z) ;hold on; plot(0:N),zeros(N+1),b);grid on; %H1 H1w=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 of H_1(z) hold on; plot(0:N),zeros(N+1),b);grid on; %G0 G0w=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 of G_0(z) hold on; plot(0:N),zeros(N+1),b);grid on; %G1 G1w=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 of G_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(ejw)|/dB,|H_1(ejw)|/dB); sumh=absH0.*absH0+absH1.*absH1; sum=10*log10(sumh); save H_filters h0n disp(Press any key to process a signal using the filters.); pause; clear all; close all; load H_filters h0=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 decimation x0=filter(h0,1,Sig); v0=decimate(x0,2); %low pass x1=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(g1,1,u1