数字信号处理课程设计报告书.docx
《数字信号处理课程设计报告书.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计报告书.docx（28页珍藏版）》请在冰豆网上搜索。
数字信号处理课程设计报告书
一、课程设计要求
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）,画出波形
程序如下
clearall
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（[-102001.2]）;
grid;
x2=（n+2>=0）-（n-2>=0）;
subplot（2,1,2）
stem（n,x2）;
xlabel（'n'）;ylabel（'x2（n）'）;
title（'单位阶跃序列'）;
axis（[-102001.2]）;
grid;
x3=（n==4）;
figure
（2）
subplot（2,1,1）
stem（n,x3）;
xlabel（'n'）;ylabel（'x2（n）'）;
title（'单位脉冲序列'）;
axis（[-102001.2]）;
grid;
x4=（n>=0）-（n-5>=0）;
subplot（2,1,2）
stem（n,x4）;
xlabel（'n'）;ylabel（'x4（n）'）;
title（'矩形窗函数'）;
axis（[-102001.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（[-1040-0.51.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;%采样周期为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;%采样周期为1.5
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';%对原函数进行傅里叶变换,
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.00040.00220.00560.00750.00560.00220.0004];
b=[1-3.18364.6223-3.77951.8136-0.480.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（[050-0.10.14]）;
figure（3）
zplane（a,b）;
x=[-4-20-4-6-4-2-4-6-6-4-4-6-6-261280-16-38-60-84-90-66-32-4-2-4812121066400000-2-200-2-2-2-20];
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（[050-5050]）;
6、IIR滤波器的设计
（1）分别设计切比雪夫II型和巴特沃什低通数字滤波器，满足一下技术指标：
通带边界角频率Wp=0.2π处，幅度最大衰减为1dB；阻带截止角频率Ws=0.5π处，幅度最小衰减为70dB。
求滤波器的最小阶数、通带（3dB）截止角频率和滤波器的系数，并画出幅频响应曲线；
（2）设计一个二阶的巴特沃斯数字低通滤波器，要求其3dB处截止角频率为π/2（rad）,取样频率为4KHZ。
分别应用脉冲响应不变法和双线性变换法设计，求出数字滤波器的系数并画出其幅频响应曲线，同时比较采用哪种变换法可以更好地满足技术指标。
（3）模拟低通滤波器的系统函数为，运用双线性变换法设计3dB截止角频率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=[11.411];
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（[0401.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（[0401.2]）
legend（'hanning'）;
三、总结
经过本次课程设计，让我熟悉了数字信号处理的基本知识，把课上的理论知识运用到实际中去，更近一步地巩固了课堂上所学的理论知识，并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。
因为学过数字信号处理这门课，但这只是理论知识，通过实验我们才能真正理解其意义。
经过这次的课程设计，让我有机会将自己学到的理论知识运用到实际中，提高了自己的动手能力和思维能力。
在课程设计中发现自己的不足，所以在今后的学习和生活中我们要更加努力，学习好我们的专业知识并要能运用到实际。
升级会员 