数字信号处理实验指导书49418doc.docx
《数字信号处理实验指导书49418doc.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书49418doc.docx（41页珍藏版）》请在冰豆网上搜索。
数字信号处理实验指导书49418doc
数字信号处理
实验指导书
武汉理工大学教材中心
2012年7月
实验一时域离散信号的产生
一、实验目的
1、了解常用时域离散信号及其特点；
2、掌握MATLAB程序的编程方法；
3、熟悉MATLAB函数的调用方法。
二、实验原理
在时间轴上的离散点取值的信号，称为离散时间信号。
离散时间信号只在某些离散的瞬时给出函数值，而在其他时刻无定义。
它是时间上不连续按一定先后次序排列的一组数的集合，称为时间序列，用x（n）表示，n取整数代表时间的离散时刻。
在MATLAB中用向量来表示一个有限长度的序列。
常用离散信号：
1、单位抽样序列
2、单位阶跃序列
3、实指数序列
4、复指数序列
5、正（余）弦序列
6、随机序列
在利用计算机进行系统的研究时，经常需要产生随机信号，MATLAB提供一个工具函数rand来产生随机信号。
7、周期序列
三、实验用函数
1、stem
功能：
绘制二维图形。
调用格式：
stem（n,x）；n为横轴，x为纵轴的线性图形。
2、length
功能：
计算某一变量的长度或采样点数。
调用格式：
N=length（t）；计算时间向量t的个数并赋给变量N。
3、axis
功能：
限定图形坐标的范围。
调用格式：
axis（[x1,x2,y1,y2]）；横坐标从x1—x2，纵坐标从y1—y2。
4、zeros
功能：
产生一个全0序列。
调用格式：
x=zeros（1,n）；产生n个0的序列。
5、ones
功能：
产生一个全1序列。
调用格式：
y=ones（1,n）；产生n个1的序列。
四、参考实例
例1.1用Matlab产生单位抽样序列。
%先建立函数impseq（n1,n2,n0）
function[x,n]=impseq（n1,n2,n0）
n=[n1:
n2];
x=[（n-n0）==0];
%编写主程序调用该函数
[x,n]=impseq（-2,8,2）;
stem（n,x）
程序运行结果如图1-1所示：
图1-1单位抽样序列
例1.2实数指数序列（运算符“.^”）
Matlab程序如下：
n=[0:
10];
x=0.9.^n;
stem（n,x）
程序运行结果如图1-2所示
图1-2实数指数序列
例1.3复数指数序列（
）
Matlab程序如下：
n=[-10:
10];alpha=-0.1+0.3*j;x=exp（alpha*n）;
real_x=real（x）;image_x=imag（x）;
mag_x=abs（x）;phase_x=angle（x）;
subplot（2,2,1）;stem（n,real_x）
subplot（2,2,2）;stem（n,image_x）
subplot（2,2,3）;stem（n,mag_x）
subplot（2,2,4）;stem（n,phase_x）
程序运行结果如图1-3所示
图1-3复数指数序列
例1.4正、余弦序列（
）
Matlab程序如下：
n=[0:
10];
x=3*cos（0.1*pi*n+pi/3）;
stem（n,x）
程序运行结果如图1-4所示
图1-4正、余弦序列
例1.5随机序列
rand（1,N）产生其元素在[0，1]之间均匀分布长度为N的随机序列
randn（1,N）产生均值为0，方差为1，长度为N的高斯随机序列
例1.6周期序列
如何生成周期序列
1、将一个周期复制p次；
2、借助矩阵运算、matlab下标能力。
先生成一个包含p列x（n）值的矩阵，然后用结构（:
）来把p列串接成一个长周期序列。
因为这个结构只能用于列向，最后还需要做矩阵转置获得所需序列。
Matlab程序如下：
x=[1,2,3];%一个x（n）
xn=x'*ones（1,3）%生成p列x（n）
xn=xn（:
）'%将p列串接成长列序列并转置
stem（xn）
程序运行的结果如图1-5所示
图1-5周期序列
五、实验任务
1、调试部分例题程序，掌握Matlab基本操作方法。
2、编写程序，完成下列函数波形：
1）利用zeros函数生成单位抽样序列；
2）利用zeros函数和ones函数生成单位阶跃序列；
六、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线。
实验二离散序列的基本运算
一、实验目的
1、加强MATLAB运用。
2、了解离散时间序列在时域中的基本运算。
3、熟悉相关函数的使用方法，掌握离散序列运算程序的编写方法。
二、实验原理
离散序列的时域运算包括信号的相加、相乘，信号的时域变换包括信号的移位、反折、倒相及尺度变换等。
在MATLAB中，序列的相加和相乘运算是两个向量之间的运算，因此参加运算的两个序列必须具有相同的长度，否则不能直接进行运算，需要进行相应的处理后再进行运算。
三、实验用函数
1、find
功能：
寻找非零元素的索引号。
调用格式：
find（（n>=min（n1））&（n<=max（n1）））：
在符号关系运算条件的范围内寻找非零元素的索引号。
2、fliplr
功能：
对矩阵行元素进行左右翻转。
调用格式：
x1=fliplr（x）：
将x的行元素进行左右翻转，赋给变量x1。
四、实例
1、信号的时域变换
1）序列移位
将一个离散序列进行移位，形成新的序列：
x1（n）=x（n-m）。
当m>0时，原序列向右移m位，当m<0时，原序列向左移。
%建立移位函数（sigshift（x,m,n0））
function[y,n]=sigshift（x,m,n0）
n=m+n0;
y=x;
2）序列反折
在这个运算中，x（n）以n=0为基准点，以纵轴为对称轴反折得到一个新的序列。
y（n）=|x（-n）|
在MATLAB中提供了fliplr函数实现序列反折。
%建立反折函数（sigfold（x,n））
function[y,n]=sigfold（x,n）
y=fliplr（x）;
n=-fliplr（n）;
3）序列倒相
是求一个与原序列的向量值相反，对应的时间向量不变的新序列。
4）序列的尺度变换
通过对时间轴的放大或压缩形成新的序列。
2、序列的算术运算
1）序列相加
序列相加是指两个序列中相同序号的序列值逐项对应相加，形成新的序列。
参加运算的两个序列的维数不同时
%建立通用函数
function[y,n]=sigadd（x1,n1,x2,n2）
n=min（min（n1）,min（n2））:
max（max（n1）,max（n2））;
y1=zeros（1,length（n））;
y2=y1;
y1（find（（n>=min（n1））&（n<=max（n1））==1））=x1;
y2（find（（n>=min（n2））&（n<=max（n2））==1））=x2;
y=y1+y2;
1）序列相乘
序列相加是指两个序列中相同序号的序列值逐项对应相乘，形成新的序列。
参加运算的两个序列的维数不同时处理方法与序列相加相同。
五、实验任务
1、理解序列运算的性质，了解函数语句的意义。
2、利用例题函数完成下列序列运算
1）已知x1（n）=u（n+1）（-3x2（n）=u（n-3）（-4求：
x（n）=x1（n）+x2（n）
2）已知x1（n）=3e-0.25n（-2x2（n）=u（n+1）（-3求：
x（n）=x1（n）*x2（n）
六、实验报告
1、简述实验目的和原理。
2、列写上机调试通过的程序，并描绘其波形曲线。
实验三离散卷积的原理及应用
一、实验目的
1、通过实验进一步理解卷积定理，了解卷积过程；
2、掌握应用线性卷积求解离散时间系统响应的基本方法。
二、实验原理
对于线性移不变离散系统，任意的输入信号x（n）可以用
及其位移的线性组合来表示，即
当输入为
时，系统的输出y（n）=h（n），由系统的线性移不变性质可以得到系统对x（n）的响应y（n）为
称为离散系统的线性卷积，简写为
y（n）=x（n）*h（n）
也就是说，如果已知系统的冲激响应，将输入信号与系统的冲激响应进行卷积运算，即可求得系统的响应。
三、实验用函数
1、卷积函数conv
功能：
进行两个序列的卷积运算。
调用格式：
y=conv（x,h）；用于求解两有限长序列的卷积。
2、sum
功能：
求各元素之和。
调用格式：
y=sum（x）；求序列x中各元素之和。
3、hold
功能：
控制当前图形窗口是否刷新的双向切换开关。
调用格式：
holdon：
使当前图形窗口中的图形保持且不被刷新，准备接受绘制新的图形。
holdoff：
使当前图形窗口中的图形不具备不被刷新的性质。
4、impz
功能：
求解数字系统的冲激响应。
调用格式：
[h,t]=impz（b,a）；求解数字系统的冲激响应h，取样点数为缺省值。
[h,t]=impz（b,a,n）；求解数字系统的冲激响应h，取样点数由n确定。
impz（b,a）；在当前窗口用stem（t,h）函数绘制图形。
5、dstep
功能：
求解数字系统的阶跃响应。
调用格式：
[h,t]=dstep（b,a）；求解数字系统的冲激响应h，取样点数为缺省值。
[h,t]=dstep（b,a,n）；求解数字系统的冲激响应h，取样点数由n确定。
dstep（b,a）；在当前窗口用stairs（t,h）函数绘制图形。
四、参考实例
在利用Matlab提供的卷积函数进行卷积运算时，主要是确定卷积结果的时间区间。
conv函数默认两信号的时间序列从n=0开始，卷积结果对应的时间序列也从n=0开始。
如果信号不是从0开始，则编程时必须用两个数组确定一个信号，其中，一个数组是信号波形的幅度样值，另一个数组是其对应的时间向量。
%建立一个适用于信号从任意时间开始的通用函数
function[y,ny]=sconv（x,h,nx,nh,p）
y=conv（x,h）;
n1=nx
（1）+nh
（1）;%计算y的非零样值的起点位置
n2=nx（length（x））+nh（length（h））;%计算y的非零样值的宽度
ny=n1:
p:
n2;%确定y的非零样值的时间向量
五、实验任务
已知一个IIR数字低通滤波器的系统函数公式为
输入一个矩形信号序列x=square（n/5）（-2），求该系统的响应。
六、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线。
实验四离散傅立叶级数
一、实验目的
1、加深对离散周期序列傅里叶级数基本概念的理解。
2、掌握MATLAB求解周期序列傅里叶级数变换和逆变换的方法。
3、观察离散周期序列的重复周期数对频谱特性的影响。
二、实验原理
离散时间序列x（n）满足x（n）=x（n+rN），称为离散周期序列，其中N为周期，x（n）为主值序列。
周期序列可用离散傅里叶级数表示成
n=0,1,…,N-1
其中，
是周期序列离散傅里叶级数第K次谐波分量的系数，也称为周期序列的频谱，可表示为
k=0,1,…,N-1
上面两式是周期序列的一对傅里叶级数变换对。
令
，以上两式可简写为：
三、实验用函数
1、mod
功能：
模除求余。
调用格式：
mod（x,m）：
x整除m取正余数。
2、floor
功能：
向-
舍入为整数。
调用格式：
floor（x）：
将x向-
舍入为整数。
四、实例
1、周期序列的傅里叶变换和逆变换
依据变换公式编写通用函数
1）离散傅里叶级数正变换通用函数
functionxk=dfs（xn,N）
n=[0:
1:
N-1];%n的行向量
k=n;%k的行向量
WN=exp（-j*2*pi/N）;%WN因子
nk=n'*k;%产生一个含nk值的N乘N维矩阵
WNnk=WN.^nk;%DFS矩阵
xk=xn*WNnk;%DFS系数行向量
2）离散傅里叶级数逆变换通用函数
functionxn=idfs（xk,N）
n=[0:
1:
N-1];%n的行向量
k=n;%k的行向量
WN=exp（-j*2*pi/N）;%WN因子
nk=n'*k;%产生一个含nk值的N乘N维矩阵
WNnk=WN.^（-nk）;%DFS矩阵
xn=xk*WNnk/N;%DFS系数行向量
例：
已知一个周期性矩形序列的脉冲宽度占整个周期的1/4，一个周期的采样点为16点。
用傅里叶级数求信号的幅度和相位频谱；求傅里叶级数逆变换的图形，与原信号图形进行比较。
MATLAB程序
N=16;
xn=[ones（1,N/4）,zeros（1,3*N/4）];
n=0:
N-1;
xk=dfs（xn,N）;
xn1=idfs（xk,N）;
subplot（2,2,1）;stem（n,xn）;title（'x（n）'）;
subplot（2,2,2）;stem（n,abs（xn1））;title（'idfs（|X（k）|）'）;
subplot（2,2,3）;stem（n,abs（xk））;title（'|X（k）|'）;
subplot（2,2,4）;stem（n,angle（xk））;title（'arg|X（k）|'）;
程序运行结果如图4-1
图4-1
2、周期重复次数对序列频谱的影响
理论上讲，周期序列不满足绝对可积条件，要对周期序列进行分析，可以先取K个周期进行处理，然后让K无限增大，研究其极限情况。
这样可以观察信号序列由非周期到周期变换时，频谱由连续谱逐渐向离散谱过渡的过程。
例：
已知一个矩形序列的脉冲宽度占整个周期的1/2，一个周期的采样点数为10，用傅立叶级数变换求信号的重复周期数分别为1、4、7、10时的幅度频谱。
MATLAB程序：
xn=[ones（1,5）,zeros（1,5）];
Nx=length（xn）;
Nw=1000;dw=2*pi/Nw;
k=floor（（-Nw/2+0.5）:
（Nw/2+0.5））;
forr=0:
3;
K=3*r+1;
nx=0:
（K*Nx-1）;
x=xn（mod（nx,Nx）+1）;
Xk=x*（exp（-j*dw*nx'*k））/K;
subplot（4,2,2*r+1）;stem（nx,x）
axis（[0,K*Nx-1,0,1.1]）;ylabel（'x（n）'）;
subplot（4,2,2*r+2）;plot（k*dw,abs（Xk））
axis（[-4,4,0,1.1*max（abs（Xk））]）;ylabel（'X（k）'）;
end
程序运行结果如图4-2
图4-2
从上图可以看出，信号序列的周期数越多，则频谱越是向几个频点集中，当信号周期数趋于无穷大时，频谱转化为离散谱。
五、实验任务
1、输入并运行例题程序，熟悉基本指令的使用。
2、已知一个信号序列的主值为x（n）=[0,1,2,3,2,1,0]，显示两个周期的信号序列波形。
求解：
用傅里叶级数求信号的幅度和相位频谱；求傅里叶级数逆变换的图形，与原信号图形进行比较。
六、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线。
实验五离散傅里叶变换
一、实验目的
1、加深对离散傅里叶变换基本概念的理解。
2、了解有限长序列傅里叶变换与周期序列傅里叶级数的联系。
3、熟悉相关函数的使用方法。
二、实验原理
有限长序列的傅里叶变换和逆变换
对于非周期序列，在实际中常常使用有限长序列。
有限长序列x（n）表示为
x（n）是非周期序列，但可以理解为某一周期序列的主值序列。
由离散傅立叶级数DFS和IDFS引出有限长序列的离散傅立叶正、逆变换关系式。
DFT与DFS的关系
比较两者的变换对，可以看出两者的区别仅仅是将周期序列换成了有限长序列。
有限长序列x（n）可以看作是周期序列
的一个周期；反之周期序列
可以看作是有限长序列x（n）以N为周期的周期延拓。
由于公式非常相似，在程序编写上也基本一致。
三、实例
1、已知有限长序列x（n）为：
x（n）=[0,1,2,3,4,5,6,7,8,9]，求x（n）的DFT和IDFT。
要求
1）画出序列傅里叶变换对应的|X（k）|和arg[X（k）]图形。
2）画出原信号与傅里叶逆变换IDFT[X（k）]图形进行比较。
MATLAB程序：
xn=[0,1,2,3,4,5,6,7,8,9];N=length（xn）;
n=0:
N-1;k=0:
N-1;
xk=xn*exp（-j*2*pi/N）.^（n’*k）;
x=（xk*exp（j*2*pi/N）.^（n’*k））/N;
subplot（2,2,1）;stem（n,xn）;title（‘x（n）’）;
subplot（2,2,2）;stem（n,abs（x））;title（‘IDFT|X（k）|’）;
subplot（2,2,3）;stem（k,abs（xk））;title（‘|X（k）|’）;
subplot（2,2,4）;stem（k,angle（Xk））;title（‘arg|X（k）|’）;
程序运行结果如图5-1：
图5-1
由上图可看出，与周期序列不同的是，有限长序列本身是仅有N点的离散序列，相当于周期序列的主值部分。
因此，其频谱也对应序列的主值部分，是含N点的离散序列。
2、有限长序列DFT与周期序列DFS的联系
已知周期序列的主值x（n）=[0,1,2,3,4,5]，求x（n）周期重复次数为4次时的DFS。
要求
1）画出原主值序列和信号周期序列；
2）画出序列傅里叶变换对的图形。
MATLAB程序：
xn=[0,1,2,3,4,5];N=length（xn）;
n=0:
4*N-1;k=0:
4*N-1;
xn1=xn（mod（n,N）+1）;
xk=xn1*exp（-j*2*pi/N）.^（n'*k）;
subplot（2,2,1）;stem（xn）;title（'原主值信号x（n）'）;
subplot（2,2,2）;stem（n,xn1）;title（'周期序列信号'）;
subplot（2,2,3）;stem（k,abs（xk））;title（'|X（k）|'）;
subplot（2,2,4）;stem（k,angle（xk））;title（'arg|X（k）|'）;
程序运行结果如图5-2：
图5-2
与上一个例题比较，有限长序列x（n）可以看成是周期序列的一个周期，反之，周期序列可以看成是有限长序列以N为周期的周期延拓。
频域上的情况也是相同的。
从这个意义上说，周期序列只是有限个序列值有意义。
四、实验任务
1、输入并运行例题程序，熟悉基本指令的使用。
2、验证离散傅里叶变换的线性性质。
有两个有限长序列分别为x1（n）和x2（n），长度分别为N1和N2，且y（n）=ax1（n）+bx2（n），（a，b均为常数），则该y（n）的N点DFT为
Y（k）=DFT[y（n）]=aX1（k）+bX2（k）（0<=k<=N-1）
其中：
N=max（N1,N2），X1（k）和X2（k）分别为x1（n）和x2（n）的N点DFT。
已知序列：
x1（n）=[0,1,2,4]
x2（n）=[1,0,1,0,1]
五、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线。
实验六快速傅里叶变换
一、实验目的
1、加深对快速傅里叶变换基本理论的理解。
2、了解用MATLAB语言进行快速傅里叶变换的方法。
3、掌握常用函数的调用方法。
二、实验原理
DFT是在时域和频域均为离散序列的变换方法，它适用于有限长序列。
但如果按照变换公式进行运算的话，当序列长度很大时，将占用很大的内存空间，运算时间也会很长，无法实时处理问题。
快速傅里叶变换是用于提高DFT运算的高速运算方法的统称，FFT是其中的一种，FFT不是一种新的变换形式，它仅仅只是一种快速算法。
FFT主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N的序列分解成多个短序列再进行运算，如基2算法、基4算法等等，从而可以大大缩短运算时间。
三、实验用函数
1、fft
功能：
一维基2快速傅里叶变换。
调用格式：
y=fft（x）：
利用FFT算法计算矢量x的离散傅里叶变换。
当x为2的幂次方时，采用高速基2FFT算法，否则为稍慢的混合算法。
y=fft（x,N）：
采用n点FFT。
当x的长度小于N时，FFT函数在x的尾部补零，以构成N点数据，当x的长度大于N时，FFT函数在x的尾部截断，以构成N点数据。
2、ifft
功能：
一维基2快速傅里叶逆变换。
调用格式：
与FFT调用方法相同，只需改换函数名。
3、fftshift
功能：
对FFT的输出进行重新排列，将零频分量移到频谱的中心。
调用格式：
y=fftshift（x）：
对FFT的输出进行重新排列，将零频分量移到频谱的中心。
四、实例
1、用MATLAB工具箱函数FFT进行频谱分析时需要注意：
1）函数fft的返回值y的数据结构的对称性
若已知序列x=[4,3,2,6,7,8,9,0]，求X（k）=DFT[x（n）]
利用函数fft计算，其MATLAB程序如下：
N=8;
n=0:
N-1;
xn=[4,3,2,6,7,8,9,0];
xk=fft（xn）’
程序运行结果如下：
xk=
39.0000
-10.7782-6.2929i
0+5.0000i
4.7782+7.7071i
5.0000
4.7782-7.7071i
0-5.0000i
-10.7782+6.2929i
由程序运行结果可见，xk的第一行元素对应频率值为0，第五行元素对应频率为莱奎斯特（Nyquist）频率，即标准频率值为1。
因此第一行至第五行对应的标准频率为0~1。
而第五行至第八行对应的是负频率，其K（x）值是以Nyquist频率为轴对称。
一般情况，对于N点的x（n）序列的FFT是N点的复数序列，其点n=N/2+1对应Nyquist频率，作谱分析时仅取序列X（k）的前一半即可，其后一半序列和前一半是对称的。
2）频率计算
若N点序列x（n）（n=0,1,…,N-1）是在采样频率fs（Hz）下获得的。
它的FFT也是N点序列，即X（k）（k=0,1,…,N-1），则第K点对应实际频率值为:
f=k*fs/N
3）作FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。
2、已知信号由15Hz幅值0.5的正弦信号和40Hz幅值2的正弦信号组成，数据采样频率为100Hz,试绘制N=128点DFT的幅频图。
MATLAB程序如下：
fs=100;
N=128;
n=0:
N-1;
t=n/fs;
x=0.5*si
升级会员 