数字信号处理实验指导.docx
《数字信号处理实验指导.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导.docx（26页珍藏版）》请在冰豆网上搜索。
数字信号处理实验指导
数字信号处理
实验指导书
信息与控制工程学院电子信息工程系
前言
一、实验目的和基本要求
《数字信号处理实验》是和《数字信号处理》课程同步开设的非独立设课实验，是理论教学的深化和补充。
通过实验，使学生巩固和加深对自动控制原理理论知识的理解，进一步培养学生独立分析问题和解决问题的能力，同时注意培养学生综合设计能力、创新能力和实事求是、严谨认真的科学作风以及良好的实验习惯，为今后工作学习打下良好的基础。
通过实验学生应达到以下基本要求：
（1）通过实验验证《数字信号处理》的基本理论，并进一步巩固和加深对基本知识的理解。
（2）能根据实验指导以及相关资料，综合运用所学知识，深入钻研有关问题，学会自己独立设计实验，分析问题、解决问题，培养一定的实验研究能力和创新能力。
（3）能正确使用实验设备，掌握实验原理，熟练运用计算机处理问题。
（4）能独立撰写实验报告，准确分析实验结果，及时发现及解决实验中的问题。
二、实验项目汇总
序号
实验项目名称
学时
实验
类型
实验
要求
每组
学生
1
离散信号的MATLAB产生和显示
2
验证
必修
1
2
离散系统的差分方程及冲激响应
2
验证
必修
1
3
离散系统的频率响应分析
2
验证
必修
1
4
用FFT作信号的频谱分析
2
验证
选修
1
5
IIR数字滤波器的设计
2
综合
选修
1
6
FIR数字滤波器的设计
2
综合
选修
1
三、实验报告与考核方式
要求学生每人独立完成实验，实验结束后按照学院标准格式，自行完成实验报告并上交。
按照学校教务处对学生实验考核有关文件精神以及实验过程考勤、操作技能、实验结果和实验报告综合考核。
实验成绩占该课程平时成绩的1/3计入总成绩。
实验一常见离散信号的MATLAB产生和图形显示
实验目的：
1.熟悉MATLAB使用环境及使用方法。
2.加深对常用离散时间信号的理解。
实验原理：
1．单位抽样序列
在MATLAB中可以利用zeros（）函数实现。
注意：
MATLAB中数组下标从1开始
如果
在时间轴上延迟了k个单位，得到
即：
2．单位阶跃序列
在MATLAB中可以利用ones（）函数实现。
3．正弦序列
在MATLAB中
4．复正弦序列
在MATLAB中
5．指数序列
在MATLAB中
实验内容：
编制程序产生上述离散时间信号，绘出其图形，并回答问题。
参考MATLAB程序清单：
%ProgramP1.1
%GenerationofaUnitSampleSequence
clf;
n=-10:
20;%Generateavectorfrom-10to20
u=[zeros（1,10）1zeros（1,20）];%Generatetheunitsamplesequence
stem（n,u）;%Plottheunitsamplesequence
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'UnitSampleSequence'）;
axis（[-102001.2]）;
%ProgramP1.2（参见教材P_9页的实指数序列MATLAB程序）
%Generationofarealexponentialsequence
clf;
n=0:
35;a=1.2;K=0.2;
x=K*a.^n;
stem（n,x）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
%ProgramP1.3
%Generationofasinusoidalsequence
n=0:
40;
f=0.1;
phase=0;
A=1.5;
arg=2*pi*f*n-phase;
x=A*cos（arg）;
clf;%Clearoldgraph
stem（n,x）;%Plotthegeneratedsequence
axis（[040-22]）;
grid;
title（'SinusoidalSequence'）;
xlabel（'Timeindexn'）;
ylabel（'Amplitude'）;
axis;
实验思考题：
1.运行程序P1.1,以产生单位抽样序列。
a）其中命令clf,axis,title,xlabel和ylabel的作用是什么？
b）修改程序，产生带有延时11个样本的延迟单位抽样序列。
2.运行程序P1.2，以产生实指数序列。
a）哪个参数控制该序列的增长或衰减率？
哪个参数控制该序列的振幅？
b）若参数a小于1，会发生什么情况？
将参数a改为0.9，将参数k改为20，再次运行程序P1.2.
3.运行程序P1.3，以产生正弦序列。
a）该序列的频率是多少？
怎样改变它？
哪个参数控制该序列的相位？
哪个参数控制该序列的振幅？
该序列的周期是多少？
b）用plot命令代替stem命令，重新运行程序，显示的图形与原图形有什么区别？
实验二离散系统的差分方程、冲激响应和卷积分析
实验目的：
1.学习用差分方程判断系统是否为线性的方法。
2.进一步加深对离散系统冲激响应的理解。
3.熟悉离散系统的卷积分析方法。
实验原理：
离散系统
其输入、输出关系可用以下差分方程描述：
输入信号分解为冲激信号，
。
系统单位冲激响应
，则系统响应为如下的卷积计算式：
当
时，h[n]是有限长度的（n：
[0，M]），称系统为FIR系统；反之，称系统为IIR系统。
在MATLAB中，可以用函数y=Filter（p,d,x）求解差分方程，也可以用函数y=Conv（x,h）计算卷积。
实验内容：
一、判断如下系统是否是线性系统
假定系统由下列差分方程描述：
y（n）-0.4y（n-1）+0.75y（n-2）=2.2403x（n）+2.4908x（n-1）+2.2403x（n-2）
试用MATLAB程序仿真该系统，输入三个不同的输入序列x1（n）,x2（n）和x（n）=ax1（n）+bx2（n）,计算相应的输出响应y1（n）,y2（n）和y（n）.从而判断该系统是否是线性系统。
%ProgramP2.1
%Generatetheinputsequences
clf;
n=0:
40;
a=2;b=-3;
x1=cos（2*pi*0.1*n）;
x2=cos（2*pi*0.4*n）;
x=a*x1+b*x2;
num=[2.24032.49082.2403];
den=[1-0.40.75];
ic=[00];%Setzeroinitialconditions
y1=filter（num,den,x1,ic）;%Computetheoutputy1[n]
y2=filter（num,den,x2,ic）;%Computetheoutputy2[n]
y=filter（num,den,x,ic）;%Computetheoutputy[n]
yt=a*y1+b*y2;
d=y-yt;%Computethedifferenceoutputd[n]
%Plottheoutputsandthedifferencesignal
subplot（3,1,1）
stem（n,y）;
ylabel（'Amplitude'）;
title（'OutputDuetoWeightedInput:
a\cdotx_{1}[n]+b\cdotx_{2}[n]'）;
subplot（3,1,2）
stem（n,yt）;
ylabel（'Amplitude'）;
title（'WeightedOutput:
a\cdoty_{1}[n]+b\cdoty_{2}[n]'）;
subplot（3,1,3）
stem（n,d）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'DifferenceSignal'）;
二、线性时不变离散时间系统的冲激响应的计算
系统：
y（n）-0.4y（n-1）+0.75y（n-2）=2.2403x（n）+2.4908x（n-1）+2.2403x（n-2）
用MATLAB函数y=impz（num,den,N）来计算系统的冲激响应。
参考程序清单：
%ProgramP2.2
%Computetheimpulseresponsey
clf;
N=40;
num=[2.24032.49082.2403];
den=[1-0.40.75];
y=impz（num,den,N）;
%Plottheimpulseresponse
stem（y）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'ImpulseResponse'）;grid;
运行程序P2.2，计算上例系统的冲激响应。
运行结果：
三、用conv命令计算离散卷积
%ProgramP2.3
clf;
h=[321-210-403];%impulseresponse
x=[1-23-4321];%inputsequence
y=conv（h,x）;
n=0:
14;
subplot（2,1,1）;
stem（n,y）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'OutputObtainedbyConvolution'）;grid;
x1=[xzeros（1,8）];
y1=filter（h,1,x1）;
subplot（2,1,2）;
stem（n,y1）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'OutputGeneratedbyFiltering'）;grid;
运行结果：
参考:
计算任意两个有限长序列的卷积和：
%Program2.4
%IllustrationofConvolution
%
a=input（'Typeinthefirstsequence='）;
b=input（'Typeinthesecondsequence='）;
c=conv（a,b）;
M=length（c）-1;
n=0:
1:
M;
disp（'outputsequence='）;disp（c）
stem（n,c）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
运行结果：
Typeinthefirstsequence=[-202-13]
Typeinthesecondsequence=[120-1]
outputsequence=
-2-425141-3
实验思考题：
1．运行程序P2.1，对由加权输入得到的y（n）与在相同系数下输出y1（n）和y2（n）相加得到的yt（n）进行比较,这两个序列是否相等?
该系统是线性系统吗?
用几组不同的权系数a和b的值以及几组不同的输入频率,重新运行程序进一步验证.
2.运行程序2.2,若想画出系统冲激响应的前45个样本,怎样修改程序?
实验三离散系统的频率响应分析和零、极点分布
实验目的：
1.加深对离散系统的频率响应分析理解。
2.加深对离散系统的零、极点分布的概念理解。
实验原理：
离散系统的时域方程为
其变换域分析方法如下：
频域
系统的频率响应为
Z域:
系统的转移函数为
分解因式
，其中
和
称为零、极点。
在MATLAB中，可以用函数[z，p，K]=tf2zp（num，den）求得有理分式形式的系统函数的零、极点；
用函数zplane（z,p）绘出零、极点分布图；
用函数zplane（num,den）直接绘出有理分式形式的系统函数的零、极点分布图。
用函数[r，p，k]=residuez（num，den）完成部分分式展开计算；
用函数sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统（即“二阶节”）的级联。
实验内容：
求系统
的零、极点，并用因式分解形式表示H（z）,并确定收敛域。
实验要求：
编程实现系统参数输入，绘出零、极点分布图和幅度频率响应曲线。
程序清单：
%Program3.1
num=input（'Typeinthenumeratorcoefficients='）;
den=input（'Typeinthedenominatorcoefficients='）;
[z,p,k]=tf2zp（num,den）;
m=abs（p）;
disp（'Zerosareat'）;disp（z）;
disp（'Polesareat'）;disp（p）;
disp（'Gainconstant'）;disp（k）;
disp（'Radiusofpoles'）;disp（m）;
sos=zp2sos（z,p,k）;
disp（'Second-ordersections'）;disp（real（sos））;
subplot（2,1,1）;
zplane（num,den）;
w=0:
pi/255:
pi;
h=freqz（num,den,w）;
subplot（2,1,2）;
plot（w/pi,abs（h））;grid
xlabel（'w/pi'）;
ylabel（'Amplitude'）;
运行结果：
Typeinthenumeratorcoefficients=[…]
Typeinthedenominatorcoefficients=[…]
Zerosareat…
Polesareat…
Gainconstant…
Radiusofpoles…
Second-ordersections…
其零极点图为：
…
实验思考题：
运行程序P3.1,解释程序的运行结果，并利用运行结果，将该系统的系统函数写成几个“二阶节”级联的形式，并画出其级联形式的结构流图。
实验四用FFT作信号的频谱分析
实验目的：
1.加深对离散信号的DTFT和DFT及其相互关系的理解。
2.对快速傅里叶变换FFT算法的理解。
3.理解用FFT对典型信号的频谱分析方法。
实验原理：
序列x[n]的DTFT定义：
N点序列x[n]的DFT定义：
在MATLAB中，对形式为
的DTDFT可以用函数H=Freqz（num，den，w）计算；可以用函数U=fft（u，N）和u=ifft（U，N）计算N点序列的DFT正、反变换。
实验内容：
1、计算N点序列u（n）的M点DFT。
程序要求输入数据N，M。
M必须大于N。
这里我们取N=8，M=16。
然后再取M=32再运行一遍程序。
比较两次的结果。
程序清单：
%Program3.1
%IllustrationofDFTComputation
%
%ReadinthelengthNofsequenceandthedesired
%lengthMoftheDFT
N=input（'Typeinthelengthofthesequence='）;
M=input（'TypeinthelengthoftheDFT='）;
%Generatethelength-Ntime-domainsequence
u=[ones（1,N）];
%ComputeitsM-pointDFT
U=fft（u,M）;
%Plotthetime-domainsequenceanditsDFT
t=0:
1:
N-1;
stem（t,u）
title（'Originaltime-domainsequence'）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
pause
subplot（2,1,1）
k=0:
1:
M-1;
stem（k,abs（U））
title（'MagnitudeoftheDFTsamples'）
xlabel（'Frequencyindexk'）;ylabel（'Magnitude'）
subplot（2,1,2）
stem（k,angle（U））
title（'PhaseoftheDFTsamples'）
xlabel（'Frequencyindexk'）;ylabel（'Phase'）;
运行结果：
原时间序列（8点）：
对其进行16点DFT：
对其进行32点DFT：
2、求正弦序列的DFT：
%Program3.2
closeall;
t=linspace（1e-3,100e-3,10）;
xn=sin（100*2*pi*t）;
N=length（xn）;
WNnk=dftmtx（N）;
Xk=xn*WNnk;
figure;
subplot（1,2,1）,stem（1:
N,xn）,title（'时域离散序列x（n）'）;
subplot（1,2,2）,stem（1:
N,abs（Xk））,title（'x（n）的DFT变换结果'）;
运行结果：
……
3．求实指数序列x（n）=an的DFT:
%Program3.3
clear
N=8;
a=0.7;
n=[0:
7];
xn=a.^n;
Xk=fft（xn,N）;
subplot（3,1,1）
stem（n,xn,'.k'）;
xlabel（'n'）;ylabel（'x（n）'）;title（'实指数序列'）;
axis（[0,8,0,1.5]）
subplot（3,1,2）
stem（n,abs（Xk）,'.k'）;
xlabel（'k'）;ylabel（'abs[X（k）]'）;title（'实指数序列8点DFT幅度'）;
axis（[0,8,0,5]）
subplot（3,1,3）
stem（n,angle（Xk）,'.k'）;
xlabel（'k'）;ylabel（'angle[X（k）]'）;title（'实指数序列8点DFT相位'）;
axis（[0,8,-1.5,1.5]）
……
实验思考题：
1.简要回答离散信号的DTFT与DFT之间的关系.
2.运行程序P3.3,通过实验结果来分析实数序列的DFT的圆周共轭对称性质.
实验五IIR数字滤波器设计
实验目的：
1.加深对IIR型数字滤波器的理解。
2.掌握IIR数字滤波器的结构和设计方法。
实验原理：
MATLAB中滤波器的分析和实现可用以下函数完成：
1、freqs函数：
模拟滤波器的频率响应；
2、freqz函数：
数字滤波器的频率响应；
3、buttord函数：
ButterWorth滤波器阶数的选择；
4、butter函数：
ButterWorth滤波器设计；
5、cheb1ord函数：
ChebyshevⅠ型滤波器阶数计算；
6、cheby1函数：
ChebyshevⅠ型滤波器设计；
7、impinvar函数：
脉冲响应不变法设计数字滤波器；
8、bilinear函数：
双线性变换法设计数字滤波器。
实验内容：
1、模拟滤波器的频率响应（freqs函数）
系统传递函数为
的模拟滤波器，在MATLAB中可以用以下程序来实现：
%Program4.1
b=[0.20.31];
a=[10.41];
w=logspace（-1,1）;%产生从10-1到101之间等间距点，即50个频率点
freqs（b,a,w）%根据输入的参数绘制幅度谱和相位谱
2、数字滤波器的频率响应（freqz函数）:
系统传递函数为
的数字滤波器在MATLAB中可以用以下程序来实现：
%Program4.2
a=[10.41];
b=[0.20.31];
freqz（b,a,128）%根据输入的参数绘制幅度谱和相位谱，得到0到π之间128个点处的频率响应
3、ButterWorth模拟和数字滤波器
（1）buttord函数：
ButterWorth滤波器阶数的选择。
调用格式：
[n,Wc]=buttord（Wp,Ws,Rp,Rs），在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内最大衰减Rp和阻带内最小衰减Rs），计算ButterWorth滤波器的阶数n和截止频率Wn。
相同参数条件下的模拟滤波器则调用格式为：
[n,Wc]=buttord（Wp,Ws,Rp,Rs,’s’）
（2）butter函数：
ButterWorth滤波器设计。
调用格式：
[b,a]=butter（n,Wc），根据阶数n和截止频率Wc计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。
相同参数条件下的模拟滤波器则调用格式为：
[b,a]=butter（n,Wc,’s’）
采样频率为1Hz，通带临界频率fp=0.2Hz，通带内衰减小于1dB（Rp=1）；阻带临界频率fs=0.3Hz，阻带内衰减大于25dB（As=25）。
设计一个数字滤波器满足以上参数。
%Program4.3
[n,Wc]=buttord（0.2,0.3,1,25）;
[b,a]=butter（n,Wc）;
freqz（b,a,512,1）;
4、Chebyshev模拟和数字滤波器
（1）cheb1ord函数：
ChebyshevⅠ型Ⅱ滤波器阶数计算。
调用格式：
[n,Wc]=cheb1ord（Wp,Ws,Rp,Rs），在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内波纹Rp和阻带内衰减Rs），选择ChebyshevⅠ型滤波器的最小阶n和截止频率Wc。
（2）cheby1函数：
ChebyshevⅠ型滤波器设计。
调用格式：
[b,a]=cheby1（n,Rp,Wc），根据阶数n、通带内波纹Rp和截止频率Wc计算ChebyshevⅠ型滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。
用ChebyshevⅠ型实现上例中的滤波器
%Program4.4
[n,Wn]=cheb1ord（0.2,0.3,1,25）;
[b,a]=cheby1（n,1,Wc
升级会员 