数字信号处理课程设计

一、课程设计要求

1.熟练掌握MATLAB语言的编程方法；

2.熟悉用于一维数字信号处理的MATLAB主要函数的应用； 3.记录实验结果（包括波形和数据），撰写课程设计报告。 二、课程设计内容

1.序列的产生

（1）编写程序产生以下序列x1(n)=(0.8)nu(n),x2(n)=u(n+2)-u(n-2),x3(n)=δ(n-4), X4(n)=R4(n),并画出波形；

（2）求卷积x1(n)* x2(n),.x1(n)*x3(n),画出波形 程序如下 clear all n=-20:20; a1=(0.8.^n) x1=a1.*(n>=0); figure(1) subplot(2,1,1) stem(n,x1);

xlabel('n');ylabel('x1(n)'); title('单位采样序列'); axis([-10 20 0 1.2]); grid;

x2=(n+2>=0)-(n-2>=0); subplot(2,1,2) stem(n,x2);

xlabel('n');ylabel('x2(n)'); title('单位阶跃序列'); axis([-10 20 0 1.2]); grid; x3=(n==4); figure(2) subplot(2,1,1) stem(n,x3);

xlabel('n');ylabel('x2(n)'); title('单位脉冲序列'); axis([-10 20 0 1.2]); grid;

x4=(n>=0)-(n-5>=0); subplot(2,1,2) stem(n,x4);

xlabel('n');ylabel('x4(n)'); title('矩形窗函数'); axis([-10 20 0 1.2]); grid;

y1=conv(x1,x2); M=length(y1)-1; n=0:1:M; figure(3) subplot(2,1,1) stem(n,y1); y2=conv(x1,x3); subplot(2,1,2) M=length(y2)-1; n=0:1:M; stem(n,y2);

2、序列的傅立叶变化

（1）设x(n)是有限长的因果序列，编写求x(n)傅里叶变换X（ejw）的函数：function[X]=dtft(x,w)

（2）验证傅里叶变换的时移和频移性质。

function[X]=dtft(x,w) N=128

n=[0:1:(N-1)]; x=0.8.^n;

w=[0:(2*pi/N):(2*pi)]; k=n'*w; X=x*exp(-j*k); figure(1) subplot(1,2,1)

plot(w,abs(X),'linewidth',2); xlabel('w/rad');ylabel('幅度'); title('幅频特性'); subplot(1,2,2)

plot(w,angle(X),'linewidth',2); xlabel('w/rad');ylabel('相位'); title('相频特性'); x1=0.8.^(n+2); X1=x1*exp(-j*k); figure(2) subplot(2,2,1)

plot(w,abs(X),'-',w,abs(X1),'-.','linewidth',2); xlabel('w/rad');ylabel('幅度'); title('幅频特性');

legend('x(n)','x(n+5)',2); subplot(2,2,2)

plot(w,angle(X),'-',w,angle(X1),'-.','linewidth',2); xlabel('w/rad');ylabel('相位'); title('相频特性');

legend('x(n)','x(n+5)',2); w0=pi/2;

x2=exp(j*w0*n).*x; X2=x2*exp(-j*k); subplot(2,2,3)

plot(w,abs(X),'-',w,abs(X2),'-.','linewidth',2);

xlabel('w/rad');ylabel('幅度'); title('幅频特性'); legend('X(w)','X(w-w0)'); subplot(2,2,4)

plot(w,angle(X),'-',w,angle(X2),'-.','linewidth',2); xlabel('w/rad');ylabel('相位'); title('相频特性');

legend('X(w)','X(w-w0)',4);

3、运用DFT分析信号的频谱

（1）对矩形窗序列x1(n)=Rs(n)进行频谱分析，截取长度N=64，记录幅频特性曲线和相频特性曲线；

（2）选取截取长度N为周期序列周期的整数倍，对以下周期序列进行频谱分析：x2(n)=cos(π/2*n)+cos(π/4*n)+cos(π/8*n),请画出x2(n)的幅频特性曲线，记下最大峰点的高度和位置；改变N值，观测峰点位置及高度变化，并与理论结果相比较。

（3）截取长度为N改为周期的非整数倍，再次分析x2(n)幅频特性，与（2）的结果作比较，说明有何区别。

（4）观测截取长度N=250时，在分别对x2(n)加矩形窗和加海明窗两种情况下，x2(n)的扶贫特性曲线，记录曲线的大致形状，并分析两种情况下的频谱泄漏的程度。 N=64; n=-50:1:500; x1=(n>=0)-(n-8>=0); X1=fft(x1,N); figure(1) subplot(1,2,1) k=0:1:(N-1);

plot(k,abs(X1),'linewidth',2); xlabel('k');ylabel('幅度'); title('幅频特性'); subplot(1,2,2)

plot(k,angle(X1),'linewidth',2); xlabel('k');ylabel('相位');

title('相频特性'); %R8的幅频和相频特性 x2=cos(pi/2*n)+cos(pi/4*n)+cos(pi/8*n); X2=fft(x2,N); k=0:1:(N-1); figure(2) subplot(3,1,1)

plot(k,abs(X2),'linewidth',2); xlabel('k');ylabel('幅度');

title('幅频特性 N=64'); %N=64 N=128; X3=fft(x2,N); k=0:1:(N-1); subplot(3,1,2)

plot(k,abs(X3),'linewidth',2); xlabel('k');ylabel('幅度');

title('幅频特性 N=128'); %N=128 N=70;

X4=fft(x2,N); k=0:1:(N-1); subplot(3,1,3)

plot(k,abs(X4),'linewidth',2); xlabel('k');ylabel('幅度');

title('幅频特性 N=70'); %N=70 N=250; W1=boxcar(N); W2=hamming(N); n=0:1:(N-1);

x2=cos(pi/2*n)+cos(pi/4*n)+cos(pi/8*n); y1=W1'.*x2; y2=W2'.*x2; k=length(y1); Y1=fft(y1,k); Y2=fft(y2,k); K=0:1:(k-1); figure(3)

plot(K,abs(fftshift(Y1)),'-',K,abs(fftshift(Y2)),'-.','linewidth',2); xlabel('k');ylabel('幅度'); title('加窗后幅频特性 N=250'); legend('矩形窗','汉明窗')；

4、取样定理的验证

（1）自选一费周期性的模拟信号Xa(t),画出该模拟信号的波形及幅度频谱。

（2）对Xa(t)进行取样，构成x(n),画出当取样频率取大小不同的值时，x(n)的幅度谱，分析不同采样频率时频谱混叠程度，并对采样定理进行验证。 t=-10:0.1:40; N=1000; w=0:N-1;

xa=(t>=0)-(t-10>=0); figure(1) subplot(2,1,1)

plot(t,xa,'linewidth',2); xlabel('t');ylabel('xa'); title('信号波形'); axis([-10 40 -0.5 1.5]); s = exp(-j*2*pi/length(w));

skn = s.^(w'*t); %F = xa*skn'; %subplot(2,1,2)

plot(w,abs(F),'linewidth',2); xlabel('w');ylabel('幅度'); title('幅频特性');

T=1; %t=-10:T:40; N=1000; w=0:N-1;

xa=(t>=0)-(t-10>=0); s = exp(-j*2*pi/length(w));

skn = s.^(w'*t); %F = xa*skn'; %figure(2) subplot(2,1,1)

plot(w,abs(F),'linewidth',2); xlabel('w');ylabel('幅度'); title('幅频特性 T=1');

T=1.5; %t=-10:T:40; N=1000; w=0:N-1;

代公式

对原函数进行傅里叶变换, 采样周期为1 代公式

对原函数进行傅里叶变换, 采样周期为1.5 xa=(t>=0)-(t-10>=0); s =exp(-j*2*pi/length(w));

skn = s.^(w'*t); %代公式

F = xa*skn'; %对原函数进行傅里叶变换, subplot(2,1,2)

plot(w,abs(F),'linewidth',2); xlabel('w');ylabel('幅度'); title('幅频特性 T=1.5');

5.离散LTI系统的分析

已知某离散系统的差分方程为：y(n)-3.1836y(n-1)+4.6223y(n-2)-3.7795y(n-3)+1.8136y(n-4)-0.48y(n-5)+0.0544y(n-6)=0.0004x(n)+0.0022x(n-1)+0.0056x(n-2)+0.0075x(n-3)+0.0056x(n-4)+0.0022x(n-5)+0.0004x(n-6), （1）画出该系统的幅频响应曲线；判断该系统有何功能？ （2）画出单位序列响应的波形，并记录其前十个点的值； （3）求出系统的零点，并画出零极点图，判断系统是否稳定； （4）已知受工业高频干扰的某人体心电图信号为：

X={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66, -32,-4,-2,-4,8,12,12,10,6,6,4,0,0,0,0,0,-2,-2,0,0,-2,-2,-2,-2,-0},求该信号通过上述系统后的响应y（n），并比较x(n)和y（n）的波形有何区别。 程序与图形： N=256;

a=[0.0004 0.0022 0.0056 0.0075 0.0056 0.0022 0.0004]; b=[1 -3.1836 4.6223 -3.7795 1.8136 -0.48 0.0544]; w=0:2*pi/N:2*pi; H=freqz(a,b,w); figure(1)

plot(w/(2*pi),abs(H));

xlabel('归一化频率');ylabel('幅度'); title('系统函数') grid; n=0:N-1; h=impz(a,b,N); figure(2) stem(n,h);

xlabel('n');ylabel('幅度'); title('单位脉冲响应'); axis([0 50 -0.1 0.14]); figure(3) zplane(a,b);

x=[-4 -2 0 -4 -6 -4 -2 -4 -6 -6 -4 -4 -6 -6 -2 6 12 8 0-16 -38 -60 -84 -90 -66 -32 -4 -2 -4 8 12 12 10 6 64 0 0 0 0 0 -2 -2 0 0 -2 -2 -2 -2 0]; k=0:length(x)-1; figure(4) subplot(2,1,1) stem(k,x);

xlabel('k');ylabel('幅度');

y=conv(x,h);

k1=0:(length(x)+length(h)-2); subplot(2,1,2) stem(k1,y);

xlabel('k');ylabel('幅度'); title('输出波形'); axis([0 50 -50 50]);

6、IIR滤波器的设计

（1）分别设计切比雪夫II型和巴特沃什低通数字滤波器，满足一下技术指标：通带边界角频率Wp=0.2π处，幅度最大衰减为1dB；阻带截止角频率Ws=0.5π处，幅度最小衰减为70dB。求滤波器的最小阶数、通带（3dB）截止角频率和滤波器的系数，并画出幅频响应曲线； （2）设计一个二阶的巴特沃斯数字低通滤波器，要求其3dB处截止角频率为π/2(rad),取样频率为4KHZ。分别应用脉冲响应不变法和双线性变换法设计，求出数字滤波器的系数并画出其幅频响应曲线，同时比较采用哪种变换法可以更好地满足技术指标。 （3）模拟低通滤波器的系统函数为Ha（s）=

1，运用双线性变换法设计3dB截止角

??2+ 2s+1频率Wc=0.4π（rad）的高通数字滤波器，取样间隔为1秒。求高通数字滤波器的系数，并画出其幅频响应曲线。（提示：运用1p2hp函数先将模拟低通滤波器转换为模拟高通滤波器）

程序与波形 (1) Wp=0.2*pi; Ws=0.5*pi; Rp=1; Rs=70;

[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s'); [z,p,k]=buttap(N); N Wn k

[B,A]=butter(N,Wn,'s'); w=[0:pi]; h=freqs(B,A,w); gain=20*log10(abs(h)); figure(1) plot(w,gain);

xlabel('w/rad');ylabel('幅度'); title('巴特沃斯低通滤波器'); [N,Wn]=cheb2ord(Wp,Ws,Rp,Rs,'s'); [z,p,k]=cheb2ap(N,Rs); N Wn k

[B,A]=cheby2(N,Rp,Wn,'s'); h=freqs(B,A,w); gain=20*log10(abs(h));

figure(2) plot(w,gain);

xlabel('w/rad');ylabel('幅度'); title('切比雪夫2型低通录波器');

（2）

[z,p,k]=buttap(2); fs=4000;

[B,A]=butter(2,fs/4,'s'); [num1,den1]=impinvar(B,A,4000); [h1,w]=freqz(num1,den1); Wn=2*4000*tan(fs/8); [B,A]=butter(2,Wn,'s');

[num2,den2]=bilinear(B,A,4000); [h2,w]=freqz(num2,den2); f=0:2*pi/(length(h1)-1):2*pi;

plot(f,abs(h2),'-.',f,abs(h1),'-','linewidth',2); xlabel('w/rad');ylabel('幅度'); title('2阶巴特沃斯低通滤波器'); legend('脉冲响应','双线性变换');

（3） num=[1]; den=[1 1.41 1]; wc=0.4*pi;

[B,A]=lp2hp(num,den,wc);

Wn=2/1*tan(0.4*pi/2); [num2,den2]=bilinear(B,A,1); [h2,w]=freqz(num2,den2); f=0:2*pi/(length(h2)-1):2*pi; plot(f,abs(h2),'linewidth',2); xlabel('w/rad');ylabel('幅度'); title('低通滤波器设计高通滤波器');

7、给予窗函数的FIR滤波器的设计

（1）设计一个30阶、通带截止角频率为0.6π的FIR低通数字滤波器，分别画出在加汉明窗、矩形窗、布莱克曼窗三种情况下滤波器的归一化幅频响应曲线，比较三种情况下幅频响应曲线的差异，冰球幅频响应曲线在阻带内最低点的dB值；

（2）设计一个30阶的带通数字滤波器，加汉明窗，通带截止角频率为[0.3π，0.6π]，画出单位序列响应和幅频响应曲线；

（3）设计30阶的高通数字滤波器，加汉明窗，通带截止角频率为0.4π，画出单位序列响应和幅频响应曲线。 程序与波性： (1) N=30;

wc=0.6*pi/(2*pi);

b=fir1(N,wc,hamming(N+1));

c=fir1(N,wc,boxcar(N+1)); d=fir1(N,wc,blackman(N+1)); [h1,w]=freqz(b); [h2,w]=freqz(c); [h3,w]=freqz(d);

w=0:2*pi/(length(h1)-1):2*pi;

plot(w,abs(h1),'-.',w,abs(h2),'-',w,abs(h3),'.-.','linewidth',2); xlabel('w/rad');ylabel('幅度'); title('窗口法设计FIR低通数字滤波器'); legend('hamming','boxcar','blackman');

（2） N=30;

wc=[0.3*pi/(2*pi) 0.6*pi/(2*pi)]; b=fir1(N,wc,hamming(N+1)); n=0:1:N; figure(1) stem(n,b);

xlabel('n');ylabel('幅度'); title('单位序列响应'); [h1,w]=freqz(b);

w=0:2*pi/(length(h1)-1):2*pi; figure(2)

plot(w,abs(h1),'-','linewidth',2); xlabel('w/rad');ylabel('幅度'); title('窗口法设计FIR带通数字滤波器'); axis([0 4 0 1.2]) legend('hamming');

（3） N=30;

wc=0.4*pi/(2*pi);

b=fir1(N,wc,'high',hanning(N+1)); n=0:1:N; figure(1) stem(n,b);

xlabel('n');ylabel('幅度'); title('单位序列响应'); [h1,w]=freqz(b);

w=0:2*pi/(length(h1)-1):2*pi; figure(2)

plot(w,abs(h1),'-','linewidth',2); xlabel('w/rad');ylabel('幅度'); title('窗口法设计FIR高通数字滤波器'); axis([0 4 0 1.2]) legend('hanning');

三、总结

经过本次课程设计，让我熟悉了数字信号处理的基本知识，把课上的理论知识运用到实际中去，更近一步地巩固了课堂上所学的理论知识，并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。因为学过数字信号处理这门课，但这只是理论知识，通过实验我们才能真正理解其意义。经过这次的课程设计，让我有机会将自己学到的理论知识运用到实际中，提高了自己的动手能力和思维能力。在课程设计中发现自己的不足，所以在今后的学习和生活中我们要更加努力，学习好我们的专业知识并要能运用到实际。

三、总结

经过本次课程设计，让我熟悉了数字信号处理的基本知识，把课上的理论知识运用到实际中去，更近一步地巩固了课堂上所学的理论知识，并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。因为学过数字信号处理这门课，但这只是理论知识，通过实验我们才能真正理解其意义。经过这次的课程设计，让我有机会将自己学到的理论知识运用到实际中，提高了自己的动手能力和思维能力。在课程设计中发现自己的不足，所以在今后的学习和生活中我们要更加努力，学习好我们的专业知识并要能运用到实际。

本文来源：https://www.bwwdw.com/article/vako.html