现代信号处理大作业

现代信号处理课程大作业选择【二、三、四】二、试用用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：Fp＝1.7KHz,Fr＝2.3KHz,Fs＝8KHz,Armin≥70dB。1、设计步骤：(1)对Fp、Fr进行预畸(2)计算，判断是否等于1，即该互补滤波器是否为互补镜像滤波器(3)计算相关系数);();(''FsFrtgFsFptgrp'''*rpc'c偶数)N为(;21奇数)N为(;;lg/)16/1lg(;150152;1121;1;;])110)(110[(1213090500''02'''211－min1.0min1.0iiuqkNqqqqqkkqkkkkrpArAp;)2cos()1(21))12(sin()1(210)1(21'2mmmmmmmiuNmquNmqq;42NN;221NNN;)/1)(1(2'2'kkviii12'1212,1;12Niviii22'22,1;12Niviii(4)互补镜像滤波器的数字实现2、程序：functionfilter2()Fp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;forjj=1:Ma=0;form=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%Nisodd,u=jendab=0;form=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endbW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));endfori=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfori=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);;22iiiA;22iiiB1221,1;1)(NiZAZAZHiii22212,1;1)(NiZBZBZZHiii)];()([21)(21ZHZHZHLendfori=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfori=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LP=zeros(size(w));HP=zeros(size(w));forn=1:length(w)z=exp(j*w(n)*2*pi);H1=1;fori=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2));endH2=1/z;fori=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endLP(n)=abs((H1+H2)/2);HP(n)=abs((H1-H2)/2);endplot(w,LP,'k',w,HP,'m');%holdon;xlabel('数字频率');ylabel('幅度');3、实验结果：图1.两带滤波器4、四带滤波器组程序：functionfilterfourFp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;forjj=1:Ma=0;form=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%Nisodd,u=jendb=0;form=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));endfori=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfori=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfori=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfori=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LLP=zeros(size(w));LHP=zeros(size(w));HLP=zeros(size(w));HHP=zeros(size(w));forn=1:length(w)z=exp(j*w(n)*2*pi);H1=1;fori=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2));endH21=1;fori=1:N1H21=H21*(a(i)+z^(-4))/(1+a(i)*z^(-4));endH2=1/z;fori=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endH22=1/(z^2);fori=1:N2H22=H22*(b(i)+z^(-4))/(1+b(i)*z^(-4));endLP=((H1+H2)/2);HP=((H1-H2)/2);LLP(n)=abs((H21+H22)/2*LP);LHP(n)=abs((H21-H22)/2*LP);HHP(n)=abs((H21+H22)/2*HP);HLP(n)=abs((H21-H22)/2*HP);endplot(w,LLP,'k',w,LHP,'m',w,HLP,'g',w,HHP,'b');xlabel('数字频率');ylabel('幅度');5、实验结果：图2.四带滤波器组三、根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1）Levinson算法2）Burg算法3）ARMA模型法4）MUSIC算法1、Levinson算法分析：(1)、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：mNnxxNmnmxnxNmR10*1),()(1)(ˆ(2)、根据估计所得到的自相关函数，用下面的迭代公式估算AR模型参数：)()0(*1,2iRaRkiikkkikikkaikRaD00,,0),1(21kkkD22121)1(kkkkiaaaikkkikik,,2,1,*1,1,,111,1kkka(3)、对于AR（p）模型，按以上述递推公式迭代计算直到pk1时为止。将算出的模型参数代入下式即可得到功率谱估计值：21,21)(ˆpijwiippjxxeaeS程序：function[sigma2,a]=levinson(signal_source,p)%阶数由p确定N=length(signal_source);%确定自相关函数form=0:N-1R(m+1)=sum(conj(signal_source([1:N-m])).*signal_source([m+1:N]))/N;end%设置迭代初值a1=-R(2)/R(1);sigma2=(1-abs(a1)^2)*R(1);gamma=-a1;%开始迭代fork=2:psigma2(k)=R(1)+sum(a1.*conj(R([2:k])));D=R(k+1)+sum(a1.*R([k:-1:2]));gamma(k)=D/sigma2(k);sigma2(k)=(1-abs(gamma(k))^2)*sigma2(k);a1=[a1-gamma(k)*conj(fliplr(a1)),-gamma(k)];enda=[1a1];%计算功率谱估计值sigma2=real(sigma2);p=15;%p分别为15、30、45、60[sigma2,a]=Levinson(signal_in_complex,p);%计算功率谱f1=linspace(-0.5,0.5,512);%从-0.5到0.5生成512个等间隔数据fork=1:512S1(k)=10*log10(sigma2(end)/(abs(1+sum(a(2:end).*exp(-j*2*pi*f1(k)*[1:p])))^2));%公式(2.3.7)并以dB表示end；实验结果：图3.Levinson算法2、Burg算法分析：(1)、设输入数据序列为10)(Nnnx，，对前后向预测误差之和求偏导，得反射系数121211*11))1()(()1()(2NknkkNknkkknenenene前后向预测误差递推公式：)1()(11)()(11*nenenenekkkkkk1,,...,3,2,1,0,1,1,1,kikkkikikakiaaa(2)、重复以上步骤直至k=p，根据迭代得到的AR模型参数计算功率谱，计算功率谱的公式同上面算法。程序：function[sigma2,a]=burg(signal_source,p)N=length(signal_source);ef=signal_source;eb=signal_source;sigma2=sum(signal_source*signal_source')/N;a=[];fork=1:pefp=ef(k+1:end);ebp=eb(k:end-1);gamma(k)=2*efp*ebp'/(efp*efp'+ebp*ebp');sigma2(k+1)=(1-abs(gamma(k))^2)*sigma2(k);ef(k+1:end)=efp-gamma(k)*ebp;eb(k+1:end)=ebp-gamma(k)'*efp;a=[a-gamma(k)*conj(fliplr(a)),-gamma(k)];end;a=[1a];sigma2=real(sigma2);实验结果：图4.Burg算法3、ARMA算法分析：(1)、用x(n)通过A(z)，得到y(n)。(2)、用一无穷阶的AR模型近似MA模型。用Burg算法可得到此近似AR模型的参数以及激励白噪声的功率。一般此AR模型的阶数应大于MA模型阶数的两倍，以获得较好的近似效果。(3)、可以证明，将上一步求出的近似AR模型参数视为时间序列，则MA模型就可视为一线性预测滤波器，按最小均方误差准则就可以求出MA模型参数，方法同样可用Burg算法。这样，ARMA模型的参数就全部估计出来了，根据以下公式即可算出功率谱：2121211)(ˆpijwiiqijwiipjxxeaebeS程序：function[a,b,sigma2]=arma(signal_source,p,q)N=length(signal_source);M=N-1;r=xcorr