数字信号处理实验指导书.docx
《数字信号处理实验指导书.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书.docx（46页珍藏版）》请在冰豆网上搜索。
数字信号处理实验指导书
《数字信号处理》
余根坚编撰
闽江学院计算机实验教学中心印制
目录
目录1
前言信号处理实验简介2
实验一离散信号的时域运算与变换4
实验二因果离散线性系统的时域分析9
实验三DFT变换的性质及应用13
实验四数字低通巴特沃斯滤波器的设计18
实验五窗函数设计FIR滤波器23
实验六基于FFT的图象压缩与放大26
附录31
附录1：
MATLAB编程及绘图方法31
附录2：
MATLAB矩阵及矩阵的运算38
前言信号处理实验简介
对于一个信号处理系统来说，可以将流程分为三个阶段，首先是定义输入序列，第二是对输入的序列进行运算，第三是保存输出序列，并显示结果。
一、实验课程任务与要求
1.实验目的
信号处理实验教学是为了将学生的计算机操作能力、分析能力、软件设计能力与应用实践结合起来，引导学生由浅入深地掌握信号处理理论与开发工具，具备实际应用的信号处理软件开发与制作基础。
2.实验基本要求（以软件设计为主要表现形式）
（1）上机前应准备好实验的程序设计算法描述与关键分析内容；
（2）准备好程序测试数据和设备操作步骤，上机调试、运行；
（3）完成每个实验后进行数据与程序对比分析；
（4）递交实验结果的可执行程序、源程序并演示实验结果；
（5）写出实验报告，
二、实验学时安排
1．实验一：
离散信号的时域运算与变换（2学时）
2．实验二：
因果离散线性系统的时域分析（2学时）
3．实验三：
DFT变换的性质及应用（2学时）
4．实验四：
数字低通巴特沃斯滤波器的设计（2学时）
5.实验五：
窗函数设计FIR滤波器（2学时）
6.实验六：
基于FFT的图象压缩与放大（2学时）
三、实验报告格式
实验报告
姓名：
________________学号：
_______________实验日期：
_______________
实验题目：
实验目的：
实验内容：
实验地点：
实验结果：
（包括列出实验编写的所有文件及各项实验结果的曲线，并加注必要的说明）
结果分析：
总结：
四、实验考核
（1）实验预习报告；
（2）实验签到；
（3）上机实际操作；
（4）实验设计报告；
五、实验仪器设备要求
（1）有快速的较高性能微机和较大内存与硬盘的设备；
（2）设备数量能适应学生人数；
（3）有Matlab程序设计环境；
六、教材及参考书
1．张志勇等.精通MATLAB6.5.北京航空航天大学出版社.2003年3月
2.邹鲲等.MATLAB6.x信号处理.清华大学出版社.2002年5月
3．陈怀琛等.MATLAB及在电子信息课程中的应用.电子工业出版社.2002
4．程佩清.数字信号处理教程[M].清华大学出版社.2003
实验一离散信号的时域运算与变换
实验目的
1.熟悉MATLAB编程特点
2.了解离散序列的延迟、相加、相乘及平移、反折、及倒相变换
实验内容
1.设计一个实现序列移位的函数
将序列x（n）＝｛1,2,3,4,0,7｝，其中n＝0：
5的每一个样本都移动3个周期，移位后的序列y（n）=x（n-3）
2.序列的奇偶分解
将序列x（n）＝｛0,1,2,3,4,3,2,1,0｝，其中n＝-3：
5进行奇偶分解
用函数stem显示其奇偶序列
3.序列的加法运算
设x1（n）＝｛1,0.5,0.3,0.4｝其中n=-1：
2；x2（n）＝｛0.2,0.3,0.4,0.5,0.8,1｝其中n=-2：
3
4.实现2中序列的翻转
实验报告
1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.对本实验结果做理论计算，解释实验结果
3.总结实验体会及实验中存在的问题
思考：
1.在matlab中应如何实现序列的相乘运算
2.在matlab中应如何实现序列的倒相运算？
1.实验原理：
熟悉MATLAB
首先，在用MATLAB表示离散序列并将其可视化时，我们还应注意以下几点：
第一、离散时间序列无法用符号运算来表示，要用矩阵的形式；
第二、由于在MATLAB中矩阵的元素个数有限，因此，MATLAB无法表示无限长的序列；
第三、在绘制离散信号波形的函数stem命令，而不是plot命令。
下面是MATLAB绘制单位抽样序列及其移位序列的函数，n0是单位抽样序列的位移量，n1,n2是序列的起止时刻，，调用该函数就可以绘出单位抽样序列及其移位序列的波形图。
function[x,n]=impuls（n0,n1,n2）
%Generatesx（n）=delta（n-n0）;n=n0处建立一个单位抽样序列
%[x,n]=impuls（n0,n1,n2）
if（（n0n2）|（n1>n2））
error（'argumentsmustsatisfyn1<=n0<=n2'）
end
n=[n1:
n2];
x=[zeros（1,（n0-n1））,1,zeros（1,（n2-n0））];
将上述文件存为：
impuls.m，在命令窗口输入
>>n0=1,n1=-20,n2=30;
[x,n]=impuls（n0,n1,n2）;figure
（1），stem（n,x）
上图为右移了一个单位的单位抽样序列
其次，对于离散序列来说，序列相加、相乘是将两序列对应时间序号的值逐项相加或相乘，平移、反折、及倒相变换与连续信号的定义完全相同，但需要注意，与连续信号不同的是，在MATLAB中，离散序列的时域运算和变换不能用符号运算来实现，必须用向量表示的方法，即在MATLAB中离散序列的相加、相乘需表示成两个向量的相加、相乘，因而参加运算的两序列向量必须具有相同的维数。
•在MATLAB中，用一个向量即可表示一个有限长度的序列。
但是，这样的向量并没有包含对应的时间序号的信息。
所以，要完整地表示离散信号要用两个向量
•如：
序列:
f（k）={1,2,-1,3,4,5}
•k=0
•在MATLAB中应表示为：
•f=[1,2,-1,3,4,5]
•k=[-3,-2,-1,0,1,2]或是k＝-3：
2
2.MATLAB常用信号生成函数：
熟悉MATLAB
•★ZEROS
•功能：
产生全零阵列
•调用格式：
X=ZEROS（N）%产生N行N列的全零矩阵
•X=ZEROS（M,N）%产生M行N列的全零矩阵
•★ ONES
•功能：
产生全1阵列
•调用格式：
X=ONES（N）%产生N行N列的全1矩阵
•X=ONES（M,N）%产生M行N列的全1矩阵
•★ SINC
•功能：
辛格函数
•调用格式：
Y=SINC（X）%
•★ RECTPULS
•功能：
产生矩形脉冲信号
•调用格式：
Y=RECTPULS（T）%产生高度为1、宽度为1、关于T＝0对称的矩形脉冲
•Y=RECTPULS（T，W）%产生高度为1、宽度为W、关于T＝0对称的矩形脉冲
•★ RAND
•功能：
产生伪随机序列
•调用格式：
Y=RAND（1,N）%产生[0,1]上均匀分布的随机序列
•Y=RANDN（1,N）%产生均值为0，方差为1的白噪声序列
•★ SAWTOOTH
•功能：
产生周期锯齿波或三角波
•调用格式：
Y=SAWTOOTH（T）%产生幅值为＋1，－1，以2为周期的方波
•Y=SAWTOOTH（T,WIDTH）%产生幅值为＋1，－1，以WIDTH*2为周期的方波
•★ SQUARE
•功能：
产生方波
•调用格式：
Y=SQUARE（T）%产生幅值为＋1，－1，以2为周期的锯齿波
•Y=SQUARE（T,DUTY）%产生幅值为＋1，－1，以占空比为DUTY的方波
•例：
t=0:
.0001:
.0625;
•y=SQUARE（2*pi*30*t，80）;plot（t,y）％产生一个占空比为80%的方波
•★ FLIPLR
•功能：
序列左右翻转
•调用格式：
Y=FLIPLR（X）
•%X=123翻转后321
•456654
•★ CUMSUM、SUM
•功能：
计算序列累加
•调用格式：
Y=CUMSUM（X）%向量X元素累加，记录每一次的累加结果,而SUM只记录最后的结果
3.参考程序
（1）实现序列移位的函数
•function[y,n]=segshift（x,n,n0）
•%功能y（n）=x（n-n0）
•%[y,n]=segshift（x,n,n0）
•n=n+n0;y=x;
（2）实现两序列相加
•function[y,n]=sigadd（x1,n1,x2,n2）
•%实现y（n）=x1（n）+x2（n）
•%-----------------------------
•%[y,n]=sigadd（x1,n1,x2,n2）
•%y=在包含n1和n2的n点上求序列和,
•%x1=在n1上的第一序列
•%x2=在n2上的第二序列（n2可与n1不等）
•n=min（min（n1）,min（n2））:
max（max（n1）,max（n2））;%y（n）的长度
•y1=zeros（1,length（n））;y2=y1;%初始化
•y1（find（（n>=min（n1））&（n<=max（n1））==1））=x1;%具有y的长度的x1
•y2（find（（n>=min（n2））&（n<=max（n2））==1））=x2;%具有y的长度的x2
•y=y1+y2;%序列相加.
（3）实现序列奇偶分解
•m=-fliplr（n）;
•m1=min（[m,n]）;m2=max（[m,n]）;m=m1:
m2;
•x1=[zeros（1,（length（m）-length（n）））,x];
•xe=0.5*（x1+fliplr（x1））
•xo=0.5*（x1-fliplr（x1））
•subplot（1,2,1）,stem（m,xe）,ylabel（'xe'）,
•subplot（1,2,2）,stem（m,xo）,ylabel（'xo'）,
奇偶分解结果入下图所示：
（4）正弦序列翻转
n=-1:
10;
x=sin（0.4*pi*n）;
y=fliplr（x）;
n1=-fliplr（n）;
subplot（2,1,1）,stem（n,x）
subplot（2,1,2）,stem（n1,y）
效果如右图所示：
实验二因果离散线性系统的时域分析
实验目的
实现由差分方程构成的数字滤波器：
实验内容
1.用两种不同的程序计算数字滤波器的单位抽样响应
，
给定差分方程：
y（n）-y（n-1）+0.9y（n-2）=x（n）
a.编制文件himpz.m，实现数字滤波器的单位抽样响应
b.编制函数文件hfilter.m，实现数字滤波器的单位抽样响应
2.给定
计算此数字滤波器的单位阶跃响应
，并判断系统的稳定性。
用stem（n,y）画出相应的图形。
参考程序gfilter.m
实验报告
1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.对本实验结果做理论计算，解释实验结果
3.总结实验体会及实验中存在的问题
思考：
1.在编写程序gfilter.m时，调用himpz.m和hfilter.m哪一个比较合理？
2.实验中给定的滤波器是什么性质的滤波器？
1.实验原理：
熟悉MATLAB
LTI离散系统的单位抽样响应定义为：
当输入为单位抽样序列是系统产生的零状态响应，用h（n）表示。
LTI离散系统
（1）MATLAB提供了函数impz（）求离散系统的单位抽样响应，并绘制时域波形。
★IMPZ
功能：
函数将绘出由矢量a和b表示的数字滤波器在指定取样时刻范围内的抽样响应h（n）的时域波形，并求其数值解。
调用格式：
H=IMPZ（B,A）%计算离散系统冲激响应，取样点数n由函数自动选取
或H=IMPZ（B,A,N）％计算指定时刻范围内（0:
N-1）的离散系统抽样响应
说明：
由矢量a和b构成的数字滤波器的差分方程为
其中：
B=[b0,b1,…,bM,bM-1],A=[a0,a1,…,aN,aN-1]，要求A的首项不为0
（2）由于系统的单位冲激响应h（n）也就是系统输入为时系统的零状态响应，除了用上述的impz求解外。
还可以调用filter函数求h（n），此时系统的输入为单位抽样序列。
离散系统差分方程求解的filter函数
调用格式
（1）y=filter（b,a,x）
说明：
b,a是式
的系数组成的向量b=[b0,b1,b2,…bM-1]及a=[a0,a1a2,…aN-1]，x是输入向量数组，y是输出向量数组和x的长度相同，而且系数a0要保证不为零。
调用格式
（2）y=filter（b,a,x,xic）
说明：
xic是初始条件等效的输入序列，MATLAB提供了函数filtic（）来确定xic
xic调用格式xic=filtic（b,a,x,y）
说明：
b,a是差分方程的系数数组，y,x是初始条件数组，分别由y（n）,x（n）的初始条件确定，即y=[y（-1）,y（-2）,….]及x=[x（-1）,x（-2）,…]
（3）要判断系统稳定性，一个方法是必须对h（n）求累加，如果结果不为inf（无穷大），则系统稳定；另一个方法是求系统的极点（即差分方程等号右端为零时方程的解），判断极点是否都在单位园内，即模值是否小于1。
MATLAB提供了一个函数ROOTS（C）求齐次方程的根
★函数ROOTS
功能：
求多项式的根
调用格式：
X＝ROOTS（C）%计算形如＝0的根，
％其中多项式系数矩阵为C，C有N+1个元素，
2.参考程序：
熟悉MATLAB
（1）himpz.m
a=[1,-1,0.9];b=1;
impz（b,a）
title（'单位抽样响应'）
（2）hfilter.m
a=[1,-1,0.9];b=1;
%求h（n）
x=impuls（0,-10,100）;n=[-10:
100];
h=filter（b,a,x）;
subplot（2,1,1）;stem（n,h）
axis（[-10,100,-1.1,1.1]）
title（‘单位冲激响应’）;xlabel（‘n’）;ylabel（‘h（n）’）
（3）gfilter.m
%求单位阶跃响应
a=[1,-1,0.9];b=1;
x=[zeros（1,10）,ones（1,101）];
n=-10:
100;
s=filter（b,a,x）;
subplot（2,1,2）;stem（n,s）
axis（[-10,100,-.5,2.5]）
title（‘单位阶跃响应’）;xlabel（‘n’）;ylabel（‘s（n）’）
%判断系统稳定
sum（abs（h））
z=roots（a）;
magz=abs（z）
程序运行结果：
ans=14.8337
magz=
0.9487
h（n）累加的结果是一个实数，说明系统稳定；magz说明极点的模小于1，所以系统稳定。
实验三DFT变换的性质及应用
实验目的
1.实现信号的DFT变换
2.了解DFT应用：
（1）用DFT计算卷积
（2）用DFT对序列进行谱分析
实验内容
1.用三种不同的DFT程序计算
的傅立叶变换X（k），并比较三种程序的计算机运行时间
步骤：
a.用for循环语句编制函数文件dft1.m，实现循环计算X（k）；
b.编写矩阵运算的函数文件dft2.m，实现矩阵计算X（k）；
根据定义：
为序列x（n）的DFT
则
c.调用FFT函数直接计算X（K），如程序函数dft3.m
d.分别利用上述三种不同的方式编写的DFT程序计算序列x（n）的DFT变换X（k），并画出幅频和相频特性，并比较3个程序的运行时间。
2.给定
，
利用DFT实现两序列的线性卷积运算，并研究DFT的点数与混叠的关系，并用stem（n,y）画出相应的图形。
参考程序dft4.m
3.讨论序列补零及增加数据长度对信号频谱的影响
（1）求出序列x（n）=cos（0.48n）+cos（0.52n）基于有限个样点n=10的频谱
（2）求n=100时，取x（n）的前10个，后90个设为零，得到x（n）的频谱
（3）增加x（n）有效的样点数，取100个样点得到x（n）的频谱
实验报告
参考程序dft5.m
1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.对本实验结果做理论计算，解释实验结果
3.总结实验体会及实验中存在的问题
思考：
1.直接计算N=256电DFT、FFT两者在理论上的计算量有何差别？
2.补零对x（n）计算出的谱线有什么影响，而增加有效数据长度对x（n）计算出的谱线有什么影响，这两种方式说明了什么问题？
参考程序：
熟悉MATLAB
1.dft.m
function[am,pha]=dft1（x）
N=length（x）;
w=exp（-j*2*pi/N）;
fork=1:
N
sum=0;
forn=1:
N
sum=sum+x（n）*w^（（k-1）*（n-1））;
end
am（k）=abs（sum）;
pha（k）=angle（sum）;
end
2.dft2.m
function[am,pha]=dft2（x）
N=length（x）;
n=[0:
N-1];
k=[0:
N-1];
w=exp（-j*2*pi/N）;
nk=n’*k;
wnk=w.^（nk）;
Xk=x*wnk;
am=abs（Xk）;
pha=angle（Xk）;
3.dft3.m
function[amfft,phafft]=dft3（x）
N=length（x）;
Xk=fft（x）;
amfft=abs（Xk）;
phafft=angle（Xk）;
4.dft4.m
%N1+N2-1=23<32
N=32;
x=[0:
15];
xx=[x,zeros（1,16）];
h=[ones（1,8）,zeros（1,24）];
Xk=fft（xx,N）;
Hk=fft（h,N）;
Yk=Xk.*Hk;
y=ifft（Yk,N）;
n=0:
N-1;
stem（n,y）;
holdon
%N=N1=16
N1=16;
x1=[0:
15];
h1=[ones（1,8）,zeros（1,8）];
Xk1=fft（x1,N1）;
Hk1=fft（h1,N1）;
Yk1=Xk1.*Hk1;
y1=ifft（Yk1,N1）;
n1=0:
N1-1;
stem（n1,y1,’.’,’m’）;
5.dft5.m
％x（n）基于10个样点的频谱
figure
（1）
n=[0:
1:
99];
x=cos（0.48*pi*n）+cos（0.52*pi*n）;
n1=[0:
1:
9];y1=x（1:
1:
10）;
subplot（2,1,1）;stem（n1,y1）;title（'signalx（n）,0<=n<=9'）;xlabel（'n'）
axis（[0,10,-2.5,2.5]）
Y1=fft（y1）;magY1=abs（Y1（1:
1:
6））;
k1=0:
1:
5;w1=2*pi/10*k1;
subplot（2,1,2）;stem（w1/pi,magY1）;title（‘10点DFT'）;
xlabel（'w/pi'），axis（[0,1,0,10]）
%在10个样点的基础上添90个零，得到密度高的频谱
figure
（2）
n3=[0:
1:
99];y3=[x（1:
1:
10）zeros（1,90）];％添90个零。
得到100个数据
subplot（2,1,1）;stem（n3,y3）;title（'signalx（n）,0<=n<=9+90zeros'）;xlabel（'n'）
axis（[0,100,-2.5,2.5]）
Y3=fft（y3）;magY3=abs（Y3（1:
1:
51））;
k3=0:
1:
50;w3=2*pi/100*k3;
subplot（2,1,2）;stem（w3/pi,magY3）;title（'100点DFT'）;xlabel（'w/pi'）
axis（[0,1,0,10]）
%增加x（n）有效的样点数，取100个样点
figure（3）
n=[0:
1:
99];
x=cos（0.48*pi*n）+cos（0.52*pi*n）;
subplot（2,1,1）;stem（n,x）;title（'signalx（n）,0<=n<=99'）;xlabel（'n'）
axis（[0,100,-2.5,2.5]）
X=fft（x）;magX=abs（X（1:
1:
51））;
k=0:
1:
50;w=2*pi/100*k;
subplot（2,1,2）;stem（w/pi,magX）;title（'100点DFT）;xlabel（'w/pi'）
axis（[0,1,0,60]）
figure
（1）
figure
（2）
figure（3）
实验四数字低通巴特沃斯滤波器的设计
实验目的
1.掌握IIR数字滤波器的设计方法
实验原理
IIR数字滤波器的设计，主要采用间接法，即：
首先设计出低通模拟滤波器H（S）；进行频率变换，将其转换为高通、带通、带阻滤波器；再用脉冲响应不变法或双线性变换法从模拟滤波器转换为数字滤波器。
对单极点的N阶H（S）用部分展开式：
冲激响应不变法取H（S）的单阶极点Sk的指数函数
作为H（Z）的极点
双线性变换法是用
代换H（S）中的S得到H（Z）,双线性变换法可完全消除频率混叠失真但存在非线性频率失真，而冲激响应不变法存在混叠失真。
在不同的设计阶段MATLAB的信号处理工具箱都给出了相应的滤波器设计函数，这些函数代表了不同类型的逼近函数的滤波器，常用的有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器类型。
设计流程图
下图所示为IIR滤波器设计过程及所用到的MATLAB工具箱中的函数。
模拟滤波器转数字滤波器（bilinear，impinvar函数）
合成一步的设计函数
（butter，cheb1，cheb2，ellip，besself）
滤波器系数
A,B
实验内容
1.设计一个低通巴特沃斯模拟滤波器，绘
升级会员 