数字信号处理实验指导书15162.docx
《数字信号处理实验指导书15162.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书15162.docx（46页珍藏版）》请在冰点文库上搜索。
数字信号处理实验指导书15162
数字信号处理
实验讲义
二O一六年三月
目录
实验一离散时间信号的时域分析3
实验二离散时间系统的时域分析6
实验三离散时间信号的频域分析9
实验四线性时不变离散时间系统的频域分析13
实验五IIR数字滤波器的设计17
实验六FIR数字滤波器的设计24
附录AMATLAB系统的常用概念28
附录B信号处理工具箱函数32
前言
数字信号处理研究数字序列信号的表示方法，并对信号进行运算，以提取包含在其中的特殊信息。
近几十年来，由于在研究及应用两方面均取得了进展，数字信号处理领域已日趋成熟。
由于计算机的大量使用，从而很容易向学生提供实际环境，以验证所学的概念和算法。
本指导书编程语言是MATLAB，它广泛应用于高性能数值计算和可视化。
本书假定读者已具备MATLAB基础知识。
前面的一些实验帮助学生理解信号处理的重要概念，后面以设计性实验项目为主，有利于加强对原理的理解并且加强对技术的应用。
附录中给出了本书中用到的MATALB函数及简要解释。
实验一离散时间信号的时域分析
一．实验目的
熟悉MATLAB中产生信号和绘制信号的基本命令；熟悉序列的简单运算，如：
加法、标量乘法、时间反转、延时、乘法等。
二．实验相关知识准备
1．用到的MATLAB命令
运算符号：
：
.+-*/;%
基本矩阵：
ionespirandrandnzeros
基本函数：
cosexpimagreal
数据分析：
sum
二维图形：
axisgridlegendplotsubplotstairsstemtitlexlableylableclf
工具箱：
sawtoothsquare
离散时间信号用数字序列x[n]来表示，常用的信号有单位冲激信号，单位阶跃信号，正弦信号，指数信号等
三．实验内容
1．离散时间序列的产生
（1）单位抽样序列（单位冲激信号）的产生和绘制
%program1
clf;%清除图形窗口
n=-10:
20;%产生向量n，取值-10-20，间隔为1
u=[zeros（1,10）1zeros（1,20）];%生成单位冲激信号，前面10个0，后面20个0
stem（n,u）;%以n为横坐标，u为纵坐标画杆状图
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;%定义横纵坐标轴名称
title（'UnitSampleSequence'）;%标注图形名称
axis（[-102001.2]）;%定义坐标轴范围
（2）单位阶跃信号的产生用S=ones（1,N）（表示产生长度为N的一维行向量）,请采用与程序1相似的过程产生单位阶跃信号并画图。
（3）指数信号的产生可以用命令.^和exp。
%program2%生成一个复指数信号
clf;%清除图形窗口
c=-（1/12）+（pi/6）*i;%生成一个复数
k=2;%复数的模
n=0:
40;%横坐标（41个点）
x=k*exp（c*n）;%生成一个复指数信号；
subplot（2,1,1）;%子图1
stem（n,real（x））;%复数的实部图形
xlabel（‘时间序列n’）；%坐标轴横轴
ylabel（‘振幅’）;%坐标轴纵轴
title（‘实部’）；%标注图形名称
subplot（2,1,2）;%子图2
stem（n,imag（x））;%复数的虚部图形
xlabel（‘时间序列n’）；
ylabel（‘振幅’）;
title（‘虚部’）；
%program3%生成一个实指数信号
clf;
n=0:
35;%横坐标（36个点）
a=1.2;%指数函数的底
k=0.2;%指数函数的系数
x=k*a.^n;%指数函数表达式
stem（n,x）;%画出实指数函数的杆状图
xlabel（‘时间序列n’）；%横坐标
ylabel（‘振幅’）;%纵坐标
（4）正余弦信号的产生和绘制，可使用函数sin和cos
%program4%产生余弦信号
n=0:
40;%横坐标（41个点）
f=0.1;%余弦信号的频率
phase=0;%余弦信号的初相
A=1.5;%余弦信号的振幅
arg=2*pi*f*n-phase;%余弦信号的相位
x=A*cos（arg）;%余弦信号
clf;%Clearoldgraph
stem（n,x）;%画出余弦信号
axis（[040-22]）;%定义坐标轴的范围
grid;%生成网格
title（'CosineSequence'）;
xlabel（'Timeindexn'）;
ylabel（'Amplitude'）;
（5）随机信号的产生可以使用x=rand（1,N）（表示在区间[0,1]中均匀分布的长度为N的随机信号）和x=randn（1,N）（表示长度为N的，且具有零均值和单位方差的正态分布的随机信号）
要求：
请自己编写程序用函数rand（1,N）和randn（1,N），生成N=20的随机信号，并画图，要求标注横纵坐标和图形名称。
（6）其它类型的信号如方波信号和锯齿信号可以用函数squre和sawtooth产生。
要求：
请自己编写程序用函数squre和sawtooth，生成N=20的方波信号和锯齿信号，并画图，要求标注横纵坐标和图形名称。
2．离散时间序列的简单运算
离散时间信号的简单运算包括加法（+）、标量乘除法（*,/）、乘除法（.*,./）、时间延迟和反转等。
例1：
数字信号处理应用的一个常见例子是从被加性噪声污染的信号中移除噪声。
假定信号s[n]被噪声d[n]所污染，得到一个含有噪声的信x[n]=s[n]+d[n]。
我们需要对x[n]进行运算，产生一个合理的逼近s[n]，对时刻n的样本求平均，产生输出信号是一种简单有效的方法。
如：
三点滑动平均的信号。
实现三点滑动平均的信号运算：
y[n]=（x（n-1）+x（n）+x（n+1））/3，程序如下:
%program5%SignalSmoothingbyAveraging
clf;
R=51;
d=0.8*（rand（1,R）-0.5）;%产生随机噪声，51个元素的行向量
m=0:
R-1;%横坐标，R个元素
s=2*m.*（0.9.^m）;%未被污染的信号，注意（行）向量之间用点乘
x=s+d;%被噪声污染的信号
subplot（2,1,1）;%子图1
plot（m,d,'r-',m,s,'g--',m,x,'b-.'）;%画出噪声（红色），信号（绿色）和它们的叠加（蓝色）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
legend（'d[n]','s[n]','x[n]'）;
x1=[00x];x2=[0x0];x3=[x00];%采用三点滑动
y=（x1+x2+x3）/3;%求三点滑动平均值
subplot（2,1,2）;%子图2
plot（m,y（2:
R+1）,'r-',m,s,'g--'）;%画出两种方法得出的图形，并比较（红色为y绿色为s）
legend（'y[n]','s[n]'）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
例2：
复杂信号的产生。
更复杂的信号可以通过在简单信号上执行基本的运算来产生，如振幅调制信号可以用低频调制信号来调制高频信号来得到。
%program6%产生振幅调制信号
n=0:
100;
m=0.4;
fh=0.1;
fl=0.01;
xh=sin（2*pi*fh*n）;%产生高频信号
xl=sin（2*pi*fl*n）;%产生低频信号（包络）
y=（1+m*xl）.*xh;%产生调制信号
Stem（n,y）;%画调制信号图
Xlabel（‘时间序号n’）;
ylabel（‘振幅’）
四．实验报告要求
1.按照实验要求完成相关程序的调试，并记录实验程序和结果。
2.总结本次实验结果，按照实验报告格式要求，书写实验报告。
实验二离散时间系统的时域分析
一．实验目的
通过MATLAB仿真一些简单的离散时间系统，并研究它们的时域特性。
二．实验相关知识准备
1．MATLAB命令
语言构造与调试：
breakendforifinput
基本函数：
absnum2str
多项式和内插函数：
conv
工具箱：
filterimpz
2.线性时不变系统（LTI）满足线性特性和时不变特性。
数字滤波器对单位冲激信号和单位阶跃信号的相应分别称为单位冲激响应和单位阶跃响应。
若离散时间系统的冲激响应是有限长的，则称该系统为有限冲激响应系统（FIR）；若离散时间系统的冲激响应是无限长的，则称该系统为无限冲激响应系统（IIR）。
Matlab中用函数filter来求系统的响应，impz来求系统的冲激响应。
y=filter（num,den,x）产生的输出向量与输入向量长度相同，而且初始值均为0.y=filter（num,den,x,ic）也可以计算系统的输出，其中ic=[y[-1],y[-2],...,y[-N]]为初始值，用[Y,fc]=filter（num,den,x,ic）可得到最终值。
Y=impz（num,den,N）可以计算线性时不变系统的冲激响应的前N个样本。
三．实验内容
1．线性离散时间系统
例1、设系统为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）=a.x1（n）+b.x2（n），计算并求出相应的输出响应y1[n],y2[n]和y[n]。
%program1%计算线性离散时间系统的输出响应
clf;
n=0:
40;%时间序列（41个元素）
a=2;b=-3;%信号的叠加系数
x1=cos（2*pi*0.1*n）;%输入信号1
x2=cos（2*pi*0.4*n）;%输入信号2
x=a*x1+b*x2;%输入信号1和2的叠加
num=[2.24032.49082.2403];%表示线性系统的差分方差x系数向量
den=[1-0.40.75];%表示线性系统的差分方差y系数向量
ic=[00];%设定初始条件
y1=filter（num,den,x1,ic）;%计算信号x1的响应y1[n]
y2=filter（num,den,x2,ic）;%计算信号x2的响应y2[n]
y=filter（num,den,x,ic）;%计算叠加信号的响应y[n]
yt=a*y1+b*y2;%计算响应y1[n]和y2[n]的叠加
d=y-yt;%计算两种响应的差d[n]
subplot（3,1,1）%Plottheoutputsandthedifferencesignal
stem（n,y）;%画出叠加信号的响应y[n]
ylabel（'Amplitude'）;
title（'OutputDuetoWeightedInput:
a\cdotx_{1}[n]+b\cdotx_{2}[n]'）;
subplot（3,1,2）
stem（n,yt）;%画出信号响应的叠加yt[n]
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'）;
请大家仔细观察两种响应，可以看出对于线性系统（满足叠加原理）两种响应相同。
2．时不变系统和时变系统
例2、因果系统的时不变特性研究，要求编写MATLAB程序仿真例1给出的系统，以产生两个不同的输入序列x（n）（在例1中x（n）=a.x1（n）+b.x2（n））和x（n+D），计算并画出相应的输出序列y1[n],y2[n+D]和y1[n]-y2[n+D]
参数：
n=0:
40;D=10;a=3.0;b=-2;
3.线性时不变系统的冲激响应
设系统为y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2]用y=impz（num,den,N）计算系统的冲激响应的前N个样本。
%program2%计算冲激响应
clf;
N=40;
num=[2.24032.49082.2403];
den=[1-0.40.75];
y=impz（num,den,N）;
stem（y）;
xlabel（'Timeindexn'）;
ylabel（'Amplitude'）;
title（'冲激响应'）;
grid;
4.卷积
假设待卷积的两个信号序列都是有限长的序列，卷积运算在Matlab中可以用函数conv实现。
例如，可以把系统的冲激响应与给定的有限长输入序列进行卷积，得到有限冲激系统的输出序列，下面为Matlab程序
%program3%计算两个信号的卷积
clf;
h=[321-210-403];%有限冲激系统的冲激响应
x=[1-23-4321];%有限长（7个元素）的输入序列
y=conv（h,x）;%计算卷积，求系统的响应
n=0:
14;%卷积运算结果的长度L等于两信号长度之和减1
Subplot（2,1,1）;
Stem（n,y）;%画出卷积计算出的系统响应
xlabel（'Timeindexn'）;
ylabel（'Amplitude'）;
title（'卷积得到的输出'）;grid;
x1=[xzeros（1,8）];%有限长（15个元素）的输入序列
y1=filter（h,1,x1）;%计算输入x1的响应，即num=h,den=1
Subplot（2,1,2）;
Stem（n,y1）;%画出滤波filter计算出的系统响应
xlabel（'Timeindexn'）;
ylabel（'Amplitude'）;
title（'由滤波生成的输出'）;grid;
要求：
请大家仔细比较两种方法计算出的系统响应
四．实验报告要求
1.按照实验要求完成相关程序的调试，并记录实验程序和结果。
2.总结本次实验结果，按照实验报告格式要求，书写实验报告。
实验三离散时间信号的频域分析
一．实验目的
在学习离散时间信号时域分析的基础上，对这些信号在频域上进行分析，从而进一步研究它们的性质。
熟悉离散时间序列频域分析的3种表示方法：
离散时间傅立叶变换（DTFT）,离散傅立叶变换（DFT）和Z变换。
二．实验相关知识准备
1.MATLAB命令
运算符和特殊字符：
<>.*^.^
语言构造与调试：
errorfunctionpause
基本函数：
angleconjrem
数据分析和傅立叶变换函数：
fftifftmaxmin
工具箱：
freqzimpzresiduezzplane
三．实验内容
1.离散时间傅立叶变换（DTFT）
序列x[n]的离散时间傅立叶变换
是
的连续函数。
在Matlab中可以用函数freqz计算序列的离散时间傅立叶变换在给定的L个离散频率点上的抽样值，我们需要尽可能大地取L的值，以使得plot命令画出的图形和真实离散傅立叶变换的图形尽可能一致。
假设
可以表示为
则freqz函数有如下几种调用方式：
（1）[H,w]=freqz（b,a,N）,其中b和a分别表示
的分子和分母多项式的系数向量。
此函数在单位圆上半部等间隔地计算N个点处的频率响应，返回该系统的N点频率矢量w和N点复数频率响应矢量H。
（2）H=freqz（b,a,w），它返回频率矢量w指定的那些频率点上的频率响应，频率范围为0~pi。
（3）[H,w]=freqz（b,a,N,‘whole’），它在整个单位圆上等间隔地计算N点频率响应，频率的范围为0~2pi。
也可以用函数abs,angle,real,image等来计算DTFT的幅度、相位以及实部和虚部。
例1：
计算离散傅立叶变换的频率样本
%program1%计算离散傅立叶变换的频率样本
clf;
w=-4*pi:
8*pi/512:
4*pi;%给定频率范围
b=[21];a=[1-0.6];%信号做DTFT后的分子和分母多项式的系数向量
h=freqz（b,a,w）;%在给定频率点w的频率响应
plot（w/pi,real（h））;%画出频率响应H的实部
grid;
title（‘H（e^{j\omega}）的实部’）；
xlabel（‘\omega/\pi’）;
ylabel（‘振幅’）；
要求：
请大家写出画频率响应H的虚部、振幅和相位的程序
例2：
已知系统y[n]-0.85y[n-1]=0.5x[n],请画出
的幅度响应和相位响应。
%program2%计算幅度响应和相位响应
clf;
b=[0.5];a=[1-0.85];%离散时间系统分子和分母多项式的系数向量
[h,w]=freqz（b,a,200,‘whole’）;%在整个单位圆上，求200个点的频率响应
magH=abs（h（1:
101））;%幅度响应
phaH=angle（h（1:
101））;%相位响应
w=w（1:
101）;%取单位圆上一半的频率点
subplot（2,1,1）;
plot（w/pi,magH）;%画出幅度响应
grid;
title（‘Magnituderesponse’）；
xlabel（‘\omega/\pi’）;
ylabel（‘振幅’）；
subplot（2,1,2）;
plot（w/pi,phaH）;%画出相位响应
grid;
title（‘phaseresponse’）；
xlabel（‘\omega/\pi’）;
ylabel（‘相位’）；
2.离散傅立叶变换（DFT）
周期序列的离散傅立叶级数（变换）定义为：
在MATLAB中，使用fft可以很容易地计算有限长序列x[n]的离散傅立叶变换。
此函数有两种形式：
y=fft（x）；y=fft（x,n）
求出时域信号x的离散傅立叶变换，n为规定的点数，n的默认值为所给x的长度。
当n取2的整数幂时变换速度最快。
通常取大于又最靠近x的幂次。
（即一般在使用fft函数前用n=2^nextpow2（length（x））得到最合适的n）。
当x的长度小于n时，fft函数在x的尾部补0，以构成长为n点数据。
当x的长度大于n时，fft函数将序列x截断，取前n点。
一般情况下，fft求出的函数多为复数，可用abs及angle分别求其幅度和相位。
注意：
栅栏效应，截断效应（频谱泄露和谱间干扰），混叠失真。
例3：
fft函数最通常的应用是计算信号的频谱。
考虑一个由100hz和200hz正弦信号构成的信号，受零均值随机信号的干扰，数据采样频率为1000hz。
通过fft函数来分析其信号频率成分。
%program3%傅立叶频率分析
t=0:
0.001:
1;%采样周期为0.001s，即采样频率为1000hz
x=sin（2*pi*100*t）+sin（2*pi*200*t）+1.5*rand（1,length（t））;%产生受噪声污染的正弦波信号
subplot（2,1,1）;
plot（x（1:
50））;%画出时域内的信号（前50个点）
y=fft（x,512）;%对信号x进行512点的fft
f=1000*（0:
256）/512;%设置频率轴（横轴）坐标（257个点），1000为采样频率
subplot（2,1,2）;%频率轴归一化后再乘采样频率即得实际信号频率
plot（f,y（1:
257））;%画出频域内的信号，若取512个点，则会频域重复
例4：
混叠现象与频率分辨率
假设现有含有三种频率成分的信号x（n）=sin（2nπf1/fs）+sin（2nπf2/fs）+sin（2nπf3/fs）,其中f1=20hz,f2=20.5hz,f3=40hz,采样频率fs=100hz.
（1）求x（n）在0～128之间的DFT的x（k）；
（2）求把
（1）中的x（n）以补0方式使其加长到0～512之间后的DFT的x（k）；
（3）求x（n）在0～512之间的DFT的x（k）。
最后分析三种形式的DFT结果图给出结论。
提示：
根据采样定理，在信号采样的过程中，如果采样频率小于信号最高频率的两倍，那么信号的各次调制频谱就会相互重叠起来，有些频率部分的幅值就与原始情况不同，因而不能分开和恢复这些部分，取样造成了信号的损失，这种频谱重叠的出现，就成为“混叠现象”。
令
为采样频率，
为信号的最高频率，为了避免混叠现象，则必须使
，而频谱分辨率，即在频率轴上所能得到的最小频率间隔F又为
。
所以，在
和
两个参数中，保持其中一个不变而增加另一个的唯一方法，就是增加在一记录长度内的点数N。
如果
与
都已给定，则N必须满足
3.一维逆快速傅立叶变换
y=ifft（x）;y=ifft（x,n）
例5：
频域采样定理的验证,要求写出Matlab程序
（1）产生一个三角波序列x（n）
（2）对M=40，计算x（n）的64点DFT，并图示x（n）和X（k）=DFT[x（n）],k=0,1,…,63
（3）对
（2）中所得X（k）在[0,2pi]上进行32点抽样得X1（k）=X（2k）,k=0,1,…,31。
并图示。
（4）求X1（k）的32点IDFT,并图示。
即x1（n）=IDFT[X1（k）],k=0,1,…,31
（5）采用周期延拓的方法绘出
的波形，评述
与x（n）的关系，并根据频域采样理论加以解释。
4.线性调频Z变换
离散傅立叶变换（DFT）可以看作信号在Z域上沿单位圆的均匀采样。
但在实际应用中，并非整个单位圆上的频谱都有意义。
一些情况下，如对于窄带信号，只希望分析信号所在的一段频带等，采样点的轨迹是一条弧线或圆周。
这种需求，就导致了线性调频Z变换（Chirpz变换）的出现。
Chirpz变换与DFT计算整个频谱的算法不同，它是一种更为灵活的计算频谱的算法，可以用来计算单位圆上任一段曲线的Z变换，作频谱分析时输入的点数和输出的点数可以不相等，从而达到频域“细化”的目的。
y=czt（x,m,w,a），
表示信号x的线性调频z变换的值y，它是信号x沿着w和a定义的螺旋线进行的z变换。
其中M是信号的长度（点数），W是z平面上感兴趣的那部分螺旋线上抽样点之间的比值，a是螺
升级会员 