数字信号处理编程作业-机电

1、数字信号处理编程作业 姓名： 学号： 单位： 专业： 导师： 实验一 离散时间系统及离散卷积1、 离散时间系统的单位脉冲响应（1）选择一个离散时间系统； y(n)=x(n)+2y(n-1)y(n)=0,n0;h(n)=0,n yy = 1 2 4 8 16 32 64 128 256 512与笔算结果完全符合。1-2 参考实例验证：离散时间系统的单位脉冲响应%单位取样序列函数function x,n=impseq(n0,n1,n2)n=n1:n2;x=(n-n0)=0;%单位阶跃序列函数function x,n=stepseq(n0,n1,n2)n=n1:n2;x=(n-n0)=0;(1) 参

2、考实例验证 x,n=impseq(0,-20,120); a=1,-1,0.9;b=1; h=filter(b,a,x); subplot(2,1,1),stem(n,h);title(冲激响应);xlabel(n);ylabel(h(n); subplot(2,1,2),zplane(b,a);title(系统零极点图);2、 离散系统的幅频、相频的分析方法(1) 参考实例验证b=0.0181,0.0543,0.0543,0.0181;a=1.000,-1.76,1.1829,-0.2781;function work1(h,N)sum=0;w=0:0.1:pi; for n=1:N sum

3、=sum+h(n).*exp(-j*n*w); endHm=abs(sum);Ha=angle(sum);subplot(2,1,1); plot(w/pi,Hm); title(幅度响应); grid;ylabel(幅度);xlabel(以pi为单位的频率);subplot(2,1,2); plot(w/pi,Ha);title(相位响应);grid;ylabel(相位/pi);xlabel(以pi为单位的频率);b=0.0181,0.0543,0.0543,0.0181;a=1.000,-1.76,1.1829,-0.2781; x,n=impseq(0,0,100);h=filter(b

4、,a,x);work1(h); 3、 离散卷积的计算(1) 参考实例验证 %矩形序列函数function x,n=rectseq(n0,n1,n2,n3)n=n2:n3;x=(n-n0=0)&(n-n1 x,n0=rectseq(0,10,-5,50);y,n1=stepseq(0,-5,50);h=y.*(0.9.n1);c=conv(x,h);subplot(3,1,1);stem(n0,x);title(输入序列); xlabel();ylabel(x(n);axis(-5,50,0,2); subplot(3,1,2);stem(n1,h);title(冲激响应); xlabel(n)

5、;ylabel(h(n);axis(-5,50,0,2); subplot(3,1,3);stem(-10:100,c);title(输出); xlabel(n);ylabel(y(n);axis(-10,100,0,8);实验二 离散傅立叶变换与快速傅立叶变换一、实验内容1、 用离散傅立叶变换程序处理时间抽样信号，并根据实序列离散傅立叶变换的对称性，初步判定程序的正确性；function Xk=dft(xn,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk; F=50;n=1:64;T=0.000625;

6、 x=cos(2*pi*F*n*T);X=dft(x,64); subplot(2,1,1);stem(n,x);xlabel(n);ylabel(x(n); subplot(2,1,2);stem(n,X);xlabel(k);ylabel(X(k);2、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列和的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。 三角波序列 反三角波序列function xc=xc()%三角波序列for n=0:7 i=n+1; if(n=0 & n=4) xc(i)=8-n; else xc(i)=0; end

7、endendN=8;XC=fft(xc,N);n=0:7;k=0:7;subplot(2,2,1);stem(n,xc);title(序列xc);xlabel(n);ylabel(xc(n);subplot(2,2,2);stem(k,XC);title(FFT变换);xlabel(k);ylabel(XC(k);subplot(2,2,3);stem(k,abs(XC);title(幅度);xlabel(k);ylabel(|XC(k)|);subplot(2,2,4);stem(k,angle(XC);title(相位);xlabel(k);xc;function xd=xd()%反三角波

8、序列for n=0:7 i=n+1; if(n=0 & n=4) xd(i)=n-3; else xd(i)=0; end endendN=8;XD=fft(xd,N);n=0:7;k=0:7;subplot(2,2,1);stem(n,xd);title(序列xd);xlabel(n);ylabel(xd(n);subplot(2,2,2);stem(k,XD);title(FFT变换);xlabel(k);ylabel(XD(k);subplot(2,2,3);stem(k,abs(XD);title(幅度);xlabel(k);ylabel(|XD(k)|);subplot(2,2,4)

9、;stem(k,angle(XD);title(相位);xlabel(k);xd;从图中的幅频相位图可以看出，它们的幅度特性是偶对称的，相位特性是奇对称。实验中的信号序列和，在单位圆上的Z变换频谱和会相同吗？如果不同，你能说出哪一个低频分量更多一些吗？为什么？function work2(xc,xd,N)XC=0;XD=0;w=0:0.1:pi; for n=1:N XC=XC+xc(n).*exp(-j*(n-1)*w); XD=XD+xd(n).*exp(-j*(n-1)*w); endHmc=abs(XC);Hmd=abs(XD);subplot(2,1,1);plot(w/pi,Hmc

10、); title(Xc(w)幅度响应); grid;subplot(2,1,2);plot(w/pi,Hmd); title(Xd(w)幅度响应); grid; xc0=1,2,3,4,4,3,2,1;xd0=4,3,2,1,1,2,3,4; work2(xc0,xd0,8);3、已知余弦信号如下 当信号频率，采样间隔，采样长度时，对该信号进行傅立叶变换。用FFT程序分析正弦信号，分别在以下情况下进行，并且分析比较结果（1） F=50, N=32,T=0.000625；（2） F=50, N=32,T=0.005；（3） F=50, N=32,T=0.0046875；（4） F=50, N=3

11、2,T=0.004；（5） F=50, N=64,T=0.000625function x,X=work3(F,N,T)n=1:N;x=cos(2*pi*F*n*T);X=fft(x,N);function work3_1(x1,x2,x3,x4,x5,X1,X2,X3,X4,X5)n1=1:32;n2=1:64;subplot(5,2,1);stem(n1,x1);xlabel(n);ylabel(x1(n);subplot(5,2,3);stem(n1,x2);xlabel(n);ylabel(x2(n);subplot(5,2,5);stem(n1,x3);xlabel(n);ylabe

12、l(x3(n);subplot(5,2,7);stem(n1,x4);xlabel(n);ylabel(x4(n);subplot(5,2,9);stem(n2,x5);xlabel(n);ylabel(x5(n); subplot(5,2,2);stem(n1,X1);xlabel(k);ylabel(X1(k);subplot(5,2,4);stem(n1,X2);xlabel(k);ylabel(X2(k);subplot(5,2,6);stem(n1,X3);xlabel(k);ylabel(X3(k);subplot(5,2,8);stem(n1,X4);xlabel(k);ylab

13、el(X4(k);subplot(5,2,10);stem(n2,X5);xlabel(k);ylabel(X5(k); x1,X1=work3(50,32,0.000625); x2,X2=work3(50,32,0.005); x3,X3=work3(50,32,0.0046875); x4,X4=work3(50,32,0.004); x5,X5=work3(50,64,0.000625); work3_1(x1,x2,x3,x4,x5,X1,X2,X3,X4,X5);4、选定某一时间信号进行N=64点离散傅立叶变换，并且，对同一信号进行快速傅立叶变换，并比较它们的速度。 format

14、long; tic;n=1:64; x=cos(2*pi/32*n);X1=dft(x,64);toc;Elapsed time is 0.000000 seconds. format long; tic;n=1:64; x=cos(2*pi/32*n);X2=fft(x,64);toc;Elapsed time is 0.000000 seconds.程序太小，程序运算时间约为0.一、 实验要求1、 调试实验程序，并且，给参考程序加注释；2、 完成实验内容2，并对结果进行分析。实验中的信号序列和，在单位圆上的Z变换频谱和会相同吗？如果不同，你能说出哪一个低频分量更多一些吗？为什么？3、 完成

15、实验内容3，并对结果进行分析；4、利用编制的计算卷积的计算程序，分别给出一下三组函数的卷积结果。（1） （2） （3） 、function y1=work4()n1=0:14; x1=1.n1; h1=(0.8).n1; %(1)n21=0:9; n22=0:19; x2=1.n21; h2=0.5*sin(0.5*n22); %(2)n3=0:9; x3=1-(0.1)*n3; h3=(0.1)*n3; %(3) %(1)l1=length(x1)+length(h1)-1;X1=fft(x1,l1); H1=fft(h1,l1);y1=ifft(X1.*H1);%(2)l2=length(

16、x2)+length(h2)-1;X2=fft(x2,l2); H2=fft(h2,l2);y2=ifft(X2.*H2);%(3)l3=length(x3)+length(h3)-1;X3=fft(x3,l3); H3=fft(h3,l3);y3=ifft(X3.*H3); % show resultfigure(1);subplot(3,1,1);stem(n1,x1);xlabel(n);ylabel(x1(n);subplot(3,1,2);stem(n1,h1);xlabel(n);ylabel(h1(n);subplot(3,1,3);stem(1:l1,real(y1);xlab

17、el(n);ylabel(y1(n);figure(2);subplot(3,1,1);stem(n21,x2);xlabel(n);ylabel(x2(n);subplot(3,1,2);stem(n22,h2);xlabel(n);ylabel(h2(n);subplot(3,1,3);stem(1:l2,real(y2);xlabel(n);ylabel(y2(n);figure(3);subplot(3,1,1);stem(n3,x3);xlabel(n);ylabel(x3(n);subplot(3,1,2);stem(n3,h3);xlabel(n);ylabel(h3(n);su

18、bplot(3,1,3);stem(1:l3,real(y3);xlabel(n);ylabel(y3(n);（1） （2） （3） 、实验三 IIR数字滤波器设计一、 实验内容如下：1、 设计一个巴特沃思数字低通滤波器，设计指标如下：通带内幅度衰减不大于1dB；阻带幅度衰减不小于15dB；%非归一化巴特沃斯模拟滤波器原型函数function b,a=buttap_o(N,omega)z,p,k=buttap(N);p=p*omega;k=k*omegaN;B=real(poly(z);b0=k;b=k*B;a=real(poly(p);%非归一化巴特沃斯模拟滤波器设计函数function b

19、,a=afd_buttap(wp,ws,Rp,As)if wp=0 error(通带边缘必须大于0);endif ws=0 error(阻带边缘必须大于通带边缘);endif(Rp=0)|(As wp=0.2*pi; ws=0.35*pi; Rp=1; As=15;T=1;op=wp/T;os=ws/T;cs,ds=afd_buttap(op,os,Rp,As);b,a=imp_invr(cs,ds,T);db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag); title(幅度响应); grid;ylabel(幅度);xlabe

20、l(以pi为单位的频率);subplot(2,2,2);plot(w/pi,pha); title(相位响应); grid;ylabel(相位);xlabel(以pi为单位的频率);subplot(2,2,3);plot(w/pi,db); title(幅度响应（dB）); grid;ylabel(对数幅度/dB);xlabel(以pi为单位的频率);subplot(2,2,4);plot(w/pi,grd); title(群延迟); grid;ylabel(样本);xlabel(以pi为单位的频率);*Butterworth Filter Order= 5二、 实验要求1、 调试实验程序，并

21、且，给实验程序加注释；2、 根据实验结果，给出自己设计的数字滤波器的幅度特性和相位特性；3、 用所设计的滤波器对不同频率的正弦波信号进行滤波，以说明其特性；4、 fp=0.2KHz,Rp=1dB, fs=0.3KHz,As=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。冲激响应不变法： fp=200; fs=300; Rp=1; As=25;T=0.001;wp=2*pi*fp*T;ws=2*pi*fs*T;op=wp/T;os=ws/T;cs,

22、ds=afd_buttap(op,os,Rp,As);b,a=imp_invr(cs,ds,T);db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag); title(幅度响应); grid;ylabel(幅度);xlabel(以pi为单位的频率);subplot(2,2,2);plot(w/pi,pha); title(相位响应); grid;ylabel(相位);xlabel(以pi为单位的频率);subplot(2,2,3);plot(w/pi,db); title(幅度响应（dB）); grid;ylabel(对数幅度/d

23、B);xlabel(以pi为单位的频率);subplot(2,2,4);plot(w/pi,grd); title(群延迟); grid;ylabel(样本);xlabel(以pi为单位的频率);*Butterworth Filter Order= 9双线性变换法： fp=200; fs=300; Rp=1; As=25;T=0.001;wp=2*pi*fp*T;ws=2*pi*fs*T;op=(2/T)*tan(wp/2);os=(2/T)*tan(ws/2);cs,ds=afd_buttap(op,os,Rp,As);b,a=bilinear(cs,ds,1/T);db,mag,pha,g

24、rd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag); title(幅度响应); grid;ylabel(幅度);xlabel(以pi为单位的频率);subplot(2,2,2);plot(w/pi,pha); title(相位响应); grid;ylabel(相位);xlabel(以pi为单位的频率);subplot(2,2,3);plot(w/pi,db); title(幅度响应（dB）); grid;ylabel(对数幅度/dB);xlabel(以pi为单位的频率);subplot(2,2,4);plot(w/pi,grd); title(群延迟

25、); grid;ylabel(样本);xlabel(以pi为单位的频率);*Butterworth Filter Order= 6实验四 FIR数字滤波器设计一、实验内容1、 设计一个FIR数字滤波器，设计指标如下：通带内幅度衰减不大于1dB；阻带幅度衰减不小于15dB；%理想低通滤波器的计算函数function hd=ideal_lp(wc,M);alpha=(M-1)/2;n=0:1:(M-1);m=n-alpha+eps;hd=sin(wc*m)./(pi*m); wp=0.2*pi;ws=0.35*pi;N=20;n=0:1:N-1;wc=(ws+wp)/2; %理想低通截止频率hd=

26、ideal_lp(wc,N); %理想低通冲激响应w_ham=(hamming(N); %Hamming窗h=hd.*w_ham; %FIR滤波器的冲激响应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); % 实际的最小阻带衰减%plotssubplot(2,2,1);stem(n,hd);title(理想冲激响应);axis(0 N-1 -0.1 0.3);xlabel(n);ylabel(hd(

27、n);subplot(2,2,2);stem(n,w_ham);title(Hamming窗);axis(0 N-1 0 1.1);xlabel(n);ylabel(w(n);subplot(2,2,3);stem(n,h);title(实际冲激响应);axis(0 N-1 -0.1 0.3);xlabel(n);ylabel(h(n);subplot(2,2,4);plot(w/pi,db);title(幅度响应（dB）);grid;axis(0 0.8 -100 0);xlabel(以pi为单位的频率);ylabel(对数幅度/dB);二、实验要求1、 调试实验程序，并且，给实验程序加注释

28、；2、 根据实验结果，记录自己设计的数字滤波器的幅度特性和相位特性；3、 比较相同设计指标情况下， FIR滤波器和IIR滤波器的幅度特性和相位特性。4、 分别用矩形窗、汉宁窗、海明窗、三角窗和Blackman窗设计一的线性相位带通滤波器，观察它们的3dB和20dB带宽，并比较五种窗的特点。%矩形窗设计FIR数字滤波器 function Rp,As=work4_2_1( N )wp=0.3*pi;ws=0.5*pi;n=0:1:N-1;wc=(ws+wp)/2; %理想低通截止频率hd=ideal_lp(wc,N); %理想低通冲激响应for i=1:N w_tr(i)=1;endh=hd.*w

29、_tr; %FIR滤波器的冲激响应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); % 实际的最小阻带衰减%plotssubplot(2,2,1);stem(n,hd);title(理想冲激响应);axis(0 N-1 -0.2 0.5);xlabel(n);ylabel(hd(n);subplot(2,2,2);stem(n,h);title(实际冲激响应);axis(0 N-1 -0.2 0.5

30、);xlabel(n);ylabel(h(n);subplot(2,2,3);plot(w/pi,db);title(幅度响应（dB）);grid;axis(0 0.8 -100 0);xlabel(以pi为单位的频率);ylabel(对数幅度/dB);subplot(2,2,4);plot(w/pi,pha);title(相位响应);grid;axis(0 1 -4 4);xlabel(以pi为单位的频率);ylabel(相位);%Hanning窗设计FIR数字滤波器function Rp,As=work4_2_2( N )wp=0.3*pi;ws=0.5*pi;n=0:1:N-1;wc=(

31、ws+wp)/2; %理想低通截止频率hd=ideal_lp(wc,N); %理想低通冲激响应w_han=hanning(N);h=hd.*w_han; %FIR滤波器的冲激响应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); % 实际的最小阻带衰减%plotssubplot(2,2,1);stem(n,hd);title(理想冲激响应);axis(0 N-1 -0.2 0.5);xlabel(n)

32、;ylabel(hd(n);subplot(2,2,2);stem(n,h);title(实际冲激响应);axis(0 N-1 -0.2 0.5);xlabel(n);ylabel(h(n);subplot(2,2,3);plot(w/pi,db);title(幅度响应（dB）);grid;axis(0 0.8 -100 0);xlabel(以pi为单位的频率);ylabel(对数幅度/dB);subplot(2,2,4);plot(w/pi,pha);title(相位响应);grid;axis(0 1 -4 4);xlabel(以pi为单位的频率);ylabel(相位);%Hamming窗设

33、计FIR数字滤波器function Rp,As=work4_2_3( N )wp=0.3*pi;ws=0.5*pi;n=0:1:N-1;wc=(ws+wp)/2; %理想低通截止频率hd=ideal_lp(wc,N); %理想低通冲激响应w_ham=hamming(N);h=hd.*w_ham; %FIR滤波器的冲激响应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); % 实际的最小阻带衰减%plo

34、tssubplot(2,2,1);stem(n,hd);title(理想冲激响应);axis(0 N-1 -0.2 0.5);xlabel(n);ylabel(hd(n);subplot(2,2,2);stem(n,h);title(实际冲激响应);axis(0 N-1 -0.2 0.5);xlabel(n);ylabel(h(n);subplot(2,2,3);plot(w/pi,db);title(幅度响应（dB）);grid;axis(0 0.8 -100 0);xlabel(以pi为单位的频率);ylabel(对数幅度/dB);subplot(2,2,4);plot(w/pi,pha);title(相位响应);grid;axis(0 1 -4 4);xlabel(以pi为单位的频率);ylabel(相位);%Blackman窗设计FIR数字滤波器function Rp,As=work4_2_4( N )wp=0.3*pi;ws=0.5*pi;n=0:1:N-1;wc=(ws+wp)/2; %理想低通截止频率hd=ideal_lp(wc,N); %理想低通冲激响应w_blk=blackman(N);h=hd.*w_blk; %FIR滤波器的冲激响应db,mag,pha,