源程序 数字信号处理课程设计.doc

1 IIR数字巴特沃思滤波器设计（冲激不变法）function IIR_Imp% Impulse Invariance Transformation% Butterworth Lowpass Filter Design%wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;OmegaP=wp*T;OmegaS=ws*T;ep=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);% Analog Butterworth Prototype Filter Calculationcs,ds=afd_butt(OmegaP,OmegaS,Rp,As);b,a=imp_invr(cs,ds,T);C,B,A=dir2par(b,a);%figure(1);db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag);title(Magnitude Response);xlabel(1frequency in pi units);ylabel(/H/);axis(0,1,0,1.1);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,0,Attn,Ripple,1);grid;subplot(2,2,3);plot(w/pi,db);title(Magnitude in dB);xlabel(2frequency in pi units);ylabel(decibels);axis(0,1,-40,5);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,-50,-15,-1,0);grid;set(gca,YTickLabelMode,manual,YTickLabels,50,15,1,0);subplot(2,2,2);plot(w/pi,pha/pi);title(Phase Response);xlabel(3frequency in pi units);ylabel(pi units);axis(0,1,-1,1);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickmode,manual,YTick,-1,0,1);grid;subplot(2,2,4);plot(w/pi,grd);title(Group Delay);xlabel(4frequency in pi units);ylabel(Samples);axis(0,1,0,10);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,0:2:10);grid; 2 IIR数字巴特沃思滤波器设计（双线性变换法）function IIR_BIL%利用双线性法设计Butterworth数字滤波器%wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;Fs=1/T;OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);ep=sqrt(10(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20);%cs,ds=afd_butt(OmegaP,OmegaS,Rp,As);b,a=bilinear(cs,ds,T);C,B,A=dir2cas(b,a);figure(1);db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag);title(Magnitude Response)xlabel(frequency in pi units);ylabel(|H|);axis(0,1,0,1.1)set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,0,Attn,Ripple,1);gridsubplot(2,2,3);plot(w/pi,db);title(Magnitude in dB);xlabel(frequency in pi units);ylabel(decibels);axis(0,1,-40,5);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,-50,-15,-1,0);grid%set(gca,YTickLabelmode,manual,YTickLabels,50:15:1:0)subplot(2,2,2);plot(w/pi,pha/pi);title(Phase Response)xlabel(frequency in pi units);ylabel(pi units);axis(0,1,-1,1);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,-1,0,1);gridsubplot(2,2,4);plot(w/pi,grd);title(Group Delay)xlabel(frequency in pi units);ylabel(Samples);axis(0,1,0,10)set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickMode,manual,YTick,0:2:10);grid3 FIR数字低通滤波器设计function FIR_DF% FIR digital Lowpass Filter Design%wp=0.2*pi;ws=0.3*pi;tr_width=ws-wpM=ceil(6.6*pi/tr_width)+1n=0:1:M-1;wc=(ws+wp)/2hd=ideal_lp(wc,M);w_ham=(hamming(M);h=hd.*w_ham;db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1)As=-round(max(db(ws/delta_w+1:1:501)figure(1)subplot(2,2,1);stem(n,hd);title(Ideal Impulse Response)axis(0 M-1 -0.1 0.3);ylabel(hd(n)subplot(2,2,2);stem(n,w_ham);title(Hamming Window)axis(0 M-1 0 1.1);ylabel(w(n)subplot(2,2,3);stem(n,h);title(Actual Impulse Response)axis(0 M-1 -0.1 0.3);ylabel(h(n)subplot(2,2,4);plot(w/pi,db);title(Magnitude Response in db);gridaxis(0 1 -100 10);ylabel(Decibels)set(gca,XTickMode,manual,XTick,0,0.2,0.3,1)set(gca,YTickMode,manual,YTick,-50,0)set(gca,YTickLabelMode,manual,YTickLabels,50,0)4 实际功率谱function yuan%首先产生实验数据%方差为sigma2=0.1的白噪声u(n)N=256;Au=sqrt(0.1);u=Au*randn(1,N);%让u(n)通过指定的滤波器B=1 -0.1 0.09 0.648;y=filter(B,1,u);%filter(B,A,X) is a(1)*y(n)+.+a(na+1)*y(n-na)= b(1)*x(n)+.+b(nb+1)*x(n-nb)%加上三个指定的正炫波t=1:N;%根据需要的信噪比求出相应的正弦波幅度,SNR=(A/sqrt(2)/sigman(db)A1=sqrt(2)*Au*10(10/20);A2=sqrt(2)*Au*10(45/20);A3=sqrt(2)*Au*10(55/20);f1=0.1;f2=0.28;f3=0.3;xn=y+A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+A3*sin(2*pi*f3*t);%画实际功率谱h,f=freqz(B,1,N);p=10*log(abs(h);p(fix(0.1*N/0.5)=p(fix(0.1*N/0.5)+10;p(fix(0.28*N/0.5)=p(fix(0.28*N/0.5)+45;p(fix(0.3*N/0.5)=p(fix(0.3*N/0.5)+55;plot(t/N/2,p-max(p);xlabel(频率);ylabel(归一化 P(jw)/dB);title(真实功率谱曲线)5 周期图法谱估计function zhou%首先产生实验数据%方差为sigma2=0.1的白噪声u(n)N=256;Au=sqrt(0.1);u=Au*randn(1,N);%让u(n)通过指定的器滤波B=1 -0.1 0.09 0.648;y=filter(B,1,u);%filter(B,A,X) is a(1)*y(n)+.+a(na+1)*y(n-na)= b(1)*x(n)+.+b(nb+1)*x(n-nb) %加上三个指定的正炫波t=1:N;%根据需要的信噪比求出相应的正璇波幅度A1=sqrt(2)*Au*10(10/20);A2=sqrt(2)*Au*10(45/20);A3=sqrt(2)*Au*10(55/20);f1=0.1;f2=0.28;f3=0.3;xn=y+A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+A3*sin(2*pi*f3*t);%画实际功率谱N=length(xn)*2;Xk=fft(xn,N);p=abs(Xk).2/length(Xk);p=10*log10(p);t=1:1:length(Xk)/2/length(Xk);p=p-max(p);plot(t,p(1:length(Xk)/2);xlabel(频率);ylabel(P(jw)/dB);6 welch法谱估计function welchN=4096;Au=sqrt(0.1);u=Au*randn(1,N);B=1 -0.1 0.09 0.648;y=filter(B,1,u);Subplot(2,1,1);Csd(u,y,N,1,64,16,linear);Subplot(2,1,2);h,f=freqz(B,1,N);h=20*log10(abs(h);plot(f./pi/2,h);课程设计指导书下载含源代码2016/4/28 22:12:37阅读575次附件下载后，解压到最后一个硬盘根目录下。需要设置matlab的当前目录。文件dir2par内容如下：function C,B,A=dir2par(b,a);%DIRECT-form to PARALLEL-form conversion %M=length(b);N=length(a);r1,p1,C=residuez(b,a);p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);if K*2=N; for i=1:2:N-2 Brow=r(i:1:i+1,:); Arow=p(i:1:i+1,:); Brow,Arow=residuez(Brow,Arow,); B(fix(i+1)/2),:)=real(Brow); A(fix(i+1)/2),:)=re