数字信号处理实验指导书-(校内讲义).doc

数字信号处理实验指导书 (校内讲义)

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

个人收集整理 勿做商业用途

《数字信号处理》
实验指导书

主编： 王春民 尹晶 郑海峰

长春理工大学光电信息学院

目 录

前 言 2
实验一 MATLAB简介 3
实验二 系统响应及系统稳定性 10
实验三 时域采样与频域采样 16
实验四 用FFT对信号作频谱分析 23
实验五 双线性变换法设计IIR数字滤波器 30
实验六 FIR数字滤波器设计与软件实现 37
实验七 数字信号处理在双音多频拨号系统中的应用 44
实验八 脉冲响应不变法设计IIR数字滤波器 51

前 言

《数字信号处理》是一门理论与实践并重的课程。它主要研究数字序列信号的表示方法，并对信号进行运算，以提取包含在其中的特殊信息.近年来，由于在研究及应用两方面均取得了进展，数字信号处理领域已经日趋成熟。为深入掌握课程内容，最好在学习理论的同时,做习题和上机实验。上机实验不仅可以帮助读者深入的理解和消化基本理论，将理论知识形象化；同时学习 MATLAB 在信号处理领域的使用,加强动手能力，积累实践经验；另外通过一些综合性实验达到综合分析、运用已学课程知识点的效果,锻炼初学者的独立解决问题的能力.本章共编写了八个实验，前五个实验属基础理论实验，后三个实验属于应用综合实验。

本课程实验:
实验一 MATLAB 简介
实验二 系统响应及系统稳定性
实验三 时域采样与频域采样
实验四 用FFT对信号作频谱分析
实验五 IIR数字滤波器设计及软件实现
实验六 FIR数字滤波器设计与软件实现
实验七 应用性实验数字信号处理在双音多频拨号系统中的应用
实验八 设计性实验脉冲响应不变法设计IIR数字滤波器

任课教师根据教学进度，安排学生上机进行实验.建议在学习完第一、二章后作实验一和实验二；在学习完第三、四章后作实验三和实验四；实验五IIR数字滤波器设计及软件实现在学习完第五章进行；实验六在学习完第六章后进行。实验七综合实验在学习完第七章或者再后些进行；实验八为设计性实验，要求同学自行设计实现。
本指导书适合电子信息类、通信类及其相关专业使用。

实验一 MATLAB简介
一、实验目的:
（1）熟悉Matlab环境.
（3）理解常用离散信号及系统
(2）掌握 Matlab 中一些基本函数的建立方法
二、实验原理：
Matlab是一套功能强大的工程计算及数据处理软件,广泛应用于工业、电子、医疗和建筑等众多领域。它是一种面向对象的，交互式程序设计语言,其结构完整又有优良的可移植性,在矩阵运算，数字信号处理方面有强大的功能。另外，Matlab提供了方便的绘图功能，便于用户直观地输出处理结果。本课程实验要求学生运用Matlab编程完成一些数字信号处理的基本功能，加深对教学内容的理解.
三、实验内容及步骤
Matlab的基本运算
1、矩阵的加、减法,矩阵的维数必须相同
>〉 A=［1 2;3 4]
A =
1 2
3 4
>> B=[5 6;7 8]
B =
5 6
7 8
〉〉 C=A+B
C =
6 8
10 12
2、矩阵的转置
>> A’
ans =
1 3
2 4
3、矩阵的乘法
〉> D=A*B
D =
19 22
43 50
4、矩阵的点乘
>〉 E=A。＊B
E =
5 12
21 32
matlab的绘图功能
1、图形窗口的建立
figure（H)：建立图形窗口函数，H为操作句柄,H为整数;
gcf：得到所建立图形的窗口句柄:H
get（H)：得到所建立的图形的各项属性的详细情况；
close（H)：关闭句柄为H的图形窗口。
2、绘制二维图形
plot（x)：若x为向量，则以x元素值为纵坐标，以相应元素下标为横坐标值，绘制连线图;若x为实数阵，则按列绘制每列元素值相对其下标的连线图,图中曲线等于x阵的列数；若x为复数阵，则分别以x实部阵和虚部阵的对应列元素为横坐标绘制多条连线图.
plot(x,y)：若x，y是同维向量，则绘制以x，y元素为横纵坐标的连线图.
3、同一窗口中不同曲线的绘制
a、为在同一窗口中绘制不同曲线，可以先将所绘制曲线的样值点数统一（通过对较少数点图形进行数值内插或补零等手段），然后将他们组合成一个矩阵,按上面的方法统一输出。
4、在一个窗口绘制不同样值点数的图形
hold on为打开；hold off则为取消。
5、分割图形窗口
subplot（m，n,k）。其中m和n表示图形窗口将分割成m行n列的子窗口，k表示将第k个子窗口作为当前的操作窗口。例：
>>x=0:pi/100:2＊pi；
>〉 figure(1);
〉> subplot(2，1,1);
>〉 plot（x，sin（x))；
〉〉 subplot（2，1，2)；
>〉 plot（x,cos（x））;
运行结果见右图.
Matlab的控制语句
1、for…end循环
格式:for 循环变量=表达式1：表达式2：表达式3
循环体
end
2、while…end循环
格式：while(条件表达式)
循环体
end
3、条件转移命令if…end
格式：if (条件表达式)
条件执行块
end
M文件的建立
用鼠标在命令窗口上点击File菜单中选择New中的m文件，如右图:
M文件有两个形式：命令文件和函数文件。
命令文件类似dos下的批处理文件，它包含一连串的MATLAB命令.
函数文件（function file)的第一句可执行语句是以function引导的定义语句，在函数文件中的变量都是局部变量.函数是一个“黑箱”,把一些数据送进去,经加工处理,把结果送出来。从形式上看，函数文件区别于命令文件之处是：命令文件的变量在文件执行完后保留在内存中；而函数文件内定义的变量仅在函数文件内部起作用，当函数文件执行完后，这些内部变量将被消除。
Matlab函数文件实际上包含五个部分：
1）、函数定义行；
2）、函数主体;
3）、函数说明；
4)、注解；
5）、H1行.
例：function y=comp（x） -—-——-——————-——-—--———函数定义
%comp(x）=(x+10）^2*x; --—--—---——-—————--—--—-H1说明
z=x+10;
y=z.^2＊x; —--—-———-———--—-----——-———-—----——--—函数主体.
%end of function and return y---——--—-———----———注解
〉> x=10
函数编制完成之后,执行以下命令
〉>x =
10
〉〉 y=comp（x)
y =
4000
函数文件就象Matlab提供的函数库一样，读者也可以建立自己的函数库及其说明文件，也可以为每一个函数编写一个使用说明,只要遵循Matlab的规定，就可以很简单的实现这种功能，这个规定如下：
1)、将自己的函数程序保存在同一个目录下；
2）、第一行执行指令的作用：指明该文件是函数文件、定义函数名、输入参数和输出参数。
3)、每个函数库的H1必须紧跟在定义行之后；
4）、每个函数库的使用说明紧跟在H1之后;
5）、在m文件前面，连续几行带‘%’的注释行有两个作用：一是随m文件全部显示打印时,直接起解释提示作用；二是提供help指令在线查询用。
常用信号与系统
1.理想采样信号序列
（1）首先产生信号x（n），O≤n≤50
n=0：50；%定义序列的长度是50
A=444.128;%设置信号有关的参数
a=50*sqrt（2）*pi;
T=0。001;
w=50*sqrt(2）＊pi;
x=A*exp(-a*n＊T）。*sin(w*n*T);%pi是Matlab定义的π，信号乘可采用”。＊”
close all;%清除已经绘制的x(n）图形
subplot(3，1,1);stem(x);%绘制x（n）的图形
title(‘理想采样信号序列’）; ％设置结果图形的标题
(2)绘制信号x（n）的幅度谱和相位谱
k=—25：25；
W=(pi/12。5)＊k；
x=x＊(exp(—j*pi/2.5）).^(n'＊k）；
magX=abs（x)；％绘制x（n)的幅度谱
subplot(3,1，2）；stem（magX);title(’序列的幅度谱’)；
angX=angle(x）；%绘制x（n)的相位谱
subplot(3，1，3）;stem（angX);title('序列的相位谱');
2。单位抽样序列
在Matlab中，这一函数可以用zeros函数实现。
n=1:50；%定义序列的长度是50
x=zeros(1，50);%注意MATLAB中数组下标是从1开始
x（1)=1;%x=［(n-k)==0］表示δ(n—k）
close all;
subplot(3,1，1);stem（x）;title（’单位冲激信号序列’）;
k=—25：25;
y=x＊(exp(-j*pi/12。5)).^（n'＊k）；
magy=abs(y);%绘制x(n)的幅度谱
subplot(3,1，2）；stem（magy)；title(’单位冲激信号的幅度谱'）；
angy=angle(y);%绘制x（n)的相位谱
subplot（3，1，3)；stem(angy）;title(’单位冲激信号的相位谱’）
3。矩形序列
n=1：50;
x=sign(sign(20—n）+1);
close all;
subplot（3，1，1)；stem（x）;title（’矩形序列'）；
k=-25：25；
X=x＊(exp(-j*pi/25))。^(n’*k）；
Magx=abs（X）;%绘制x（n)的幅度谱
subplot（3,1,2);stem(Magx）；title('矩形序列的幅度谱’)；
angx=angle(X);％绘制x（n）的相位谱
subplot(3,1,3)；stem(angx）;title（'矩形序列的相位谱')
4.特定冲激串
n=1：50；％定义序列的长度是50
x=zeros（1,50）；%注意：Matlab中数组下标从1开始
x（1)=1；
x(2）=2。5；
x（3)=2。5；
x(4）=1；
close all;
subplot(3，1,1)；stem（x);title（’特定冲激串序列'）；
k=—25:25;
X=x*(exp(-j＊pi/12.5）).^(n'＊k）；
magX=abs（x）；%绘制x（n）的幅度谱
subplot（3，1，2）;stem(magX);title('特定冲激串的幅度谱')；
angX=angle（X);%绘制x（n)的相位谱
subplot(3，1,3);stem(angX）;title(’特定冲激串的相位谱'）
5.其它序列
(1）单位阶跃序列：
在Matlab中可以利用ones(）函数实现。如：x=ones（1,N)
(2）正弦序列:x（n)=Asin（2πfn/fs+φ)
在Matlab中：n=0：N—1,x(n）=A*sin（2＊pi*f＊n/Fs+fai）
（3)复正弦序列：x(n)=ejwn
在Matlab中： n=0:N—1,x(n）=exp(j*w＊n）
（4）指数序列：x(n）=an
在Matlab中: n=0:n—1,x(n)=a。^n
常用MATLAB函数
表1.1 常用函数
序号
函数
功能
1
Abs
求绝对值
2
Angle
求相角
3
Conv
求卷积
4
Filter
直接滤波器实现
5
Filter2
二维数字滤波器
6
Freqs
模拟滤波器频率响应
7
Freqspace
画出频率响应曲线
8
Freqz
数字滤波器频率响应
9
freqzplot
画出频率响应曲线
10
impz
数字滤波器的单位冲激响应
11
zplane
离散系统零、极点图
12
unwrap
去除相位
Matlab的帮助
在matlab的命令窗口中输入:help 所需查询的函数名，即可得到所需查看的函数的详细使用资料。
信号处理工具箱函数
1）、波形产生
sawtooth：产生锯齿波或三角波；
square：产生方波；
sinc：产生sinc或函数；
diric：产生Dirichlet或周期sinc函数。
2）、滤波器分析和实现函数
abs：求绝对值（幅值）；
angle：求相角;
conv:：求卷积
fftfilt：重叠相加法FFT滤波器实现；
filter：直接滤波器实现；
filtic：filter函数初始条件选择；
freqs:模拟滤波器频率响应；
freqspace：频率响应中的频率间隔；
freqz：数字滤波器频率响应;
grpdelay:平均滤波器延迟;
impz:数字滤波器的冲激响应；
zplane：离散系统零极点图
3）、IIR滤波器设计
besself:Bessel（贝塞尔）滤波器设计；
butter：巴特沃斯滤波器设计
cheby1：切比雪夫1型滤波器设计；
cheby2：切比雪夫2型滤波器设计；
ellip：椭圆滤波器设计；
yulewalk：递归数字滤波器设计
4）、IIR滤波器的阶数选择
buttord：巴特沃斯滤波器阶数选择；
cheb1ord：切比雪夫1型滤波器阶数选择；
cheb2ord：切比雪夫2型滤波器阶数选择；
ellipord：椭圆滤波器阶数选择。
5)、FIR滤波器设计
fir1:基于窗口函数的FIR滤波器设计—-—标准响应；
fir2：基于窗口函数的FIR滤波器设计—-—任意响应；
firls：最小二乘FIR滤波器设计;
intfit：内插FIR滤波器设计；
6)、变换
czt：线性调频Z变换；
dct:离散余弦变换；
idct:逆离散余弦变换；
fft:一维快速傅立叶变换；
ifft：一维快速逆傅立叶变换;
7）、窗函数
boxcar：矩形窗;
triang:三角窗；
bartlert：巴特利特窗;
hamming:哈明窗；
hanning：汉宁窗；
blackman：布莱克曼窗；
Kaiser：凯塞窗。
8）、模拟滤波器设计
buttap：巴特沃思模拟低通滤波器设计；
cheb1ap：切比雪夫1型滤波器设计；
cheb2ap：切比雪夫2型滤波器设计；
9）、频率变换
lp2bp:低通到带通模拟滤波器变换；
lp2hp:低通到高通模拟滤波器变换；
lp2bs:低通到带阻模拟滤波器变化；
lp2lp：低通到低通模拟滤波器变换。
10)、滤波器离散化
bilinear：双线性变换;
impinvar：冲激响应不变法实现模拟到数字的滤波器变换。
实验二 系统响应及系统稳定性
一、实验目的
(1)掌握求系统响应的方法。
（2）掌握时域离散系统的时域特性。
（3）分析、观察及检验系统的稳定性。
二、实验原理与方法
在时域中,描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。在计算机上适合用递推法求差分方程的解,最简单的方法是采用MATLAB语言的工具箱函数filter函数.也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。
系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。
系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。或者系统的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。
实际中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数(包括零），就可以断定系统是稳定的［19]。系统的稳态输出是指当时,系统的输出。如果系统稳定，信号加入系统后,系统输出的开始一段称为暂态效应,随n的加大,幅度趋于稳定，达到稳态输出.
注意在以下实验中均假设系统的初始状态为零。
三、实验内容及步骤
（1）编制程序，包括产生输入信号、单位脉冲响应序列的子程序,用filter函数或conv函数求解系统输出响应的主程序.程序中要有绘制信号波形的功能.
（2)给定一个低通滤波器的差分方程为

输入信号

a） 分别求出系统对和的响应序列，并画出其波形。
b) 求出系统的单位冲响应，画出其波形。
（3）给定系统的单位脉冲响应为

用线性卷积法分别求系统h1（n）和h2(n)对的输出响应，并画出波形。
（4）给定一谐振器的差分方程为

令 ，谐振器的谐振频率为0。4rad。
a) 用实验方法检查系统是否稳定。输入信号为时，画出系统输出波形。
b） 给定输入信号为

求出系统的输出响应，并画出其波形。
四、思考题
（1) 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应? 如何求？
（2)如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面第一个实验结果进行分析说明。
五、实验报告要求
（1）简述在时域求系统响应的方法。
（2）简述通过实验判断系统稳定性的方法。分析上面第三个实验的稳定输出的波形。
（3）对各实验所得结果进行简单分析和解释。
(4)简要回答思考题。
(5)打印程序清单和要求的各信号波形。
六、实验程序清单
％实验2:系统响应及系统稳定性
close all；clear all
％======内容1：调用filter解差分方程，由系统对u（n)的响应判断稳定性======
A=［1，-0.9］;B=［0。05,0。05］； %系统差分方程系数向量B和A
x1n=[1 1 1 1 1 1 1 1 zeros（1,50）］; ％产生信号x1(n）=R8（n)
x2n=ones（1，128）； %产生信号x2(n）=u（n）
hn=impz（B，A,58)； %求系统单位脉冲响应h(n)
subplot(2,2,1）；y=’h（n)';tstem(hn，y）; ％调用函数tstem绘图
title（’(a） 系统单位脉冲响应h(n）')；box on
y1n=filter(B，A，x1n）； ％求系统对x1（n）的响应y1（n）
subplot（2，2,2)；y='y1(n)’;tstem（y1n，y）；
title(’（b） 系统对R8(n）的响应y1（n）'）;box on
y2n=filter（B，A,x2n）; %求系统对x2(n）的响应y2(n）
subplot(2,2,4）；y=’y2（n）'；tstem（y2n,y)；
title（'（c) 系统对u(n)的响应y2（n）'）;box on
％===内容2：调用conv函数计算卷积============================
x1n=[1 1 1 1 1 1 1 1 ］; %产生信号x1（n）=R8(n)
h1n=[ones（1,10) zeros(1，10)]；
h2n=［1 2。5 2.5 1 zeros（1，10)］；
y21n=conv(h1n,x1n）;
y22n=conv(h2n，x1n）；
figure(2）
subplot（2，2,1)；y=’h1（n）’；tstem（h1n，y）； %调用函数tstem绘图
title（’（d） 系统单位脉冲响应h1（n)’）；box on
subplot(2,2,2）;y='y21（n）’；tstem（y21n，y);
title(’(e） h1（n）与R8（n)的卷积y21(n）’）;box on
subplot(2,2，3）;y='h2（n)’；tstem（h2n，y)； %调用函数tstem绘图
title（'(f) 系统单位脉冲响应h2(n)’）；box on
subplot（2，2,4）;y='y22（n）'；tstem(y22n,y)；
title('（g) h2(n)与R8(n）的卷积y22（n)'）;box on
%=========内容3：谐振器分析========================
un=ones（1，256); %产生信号u（n)
n=0:255;
xsin=sin（0.014*n）+sin(0.4＊n）； %产生正弦信号
A=［1,—1。8237，0。9801];B=［1/100。49,0,—1/100。49]； ％系统差分方程系数向量B和A
y31n=filter(B,A,un)； %谐振器对u(n）的响应y31（n)
y32n=filter(B,A，xsin); %谐振器对u（n)的响应y31(n)
figure（3）
subplot（2，1，1）;y='y31(n）';tstem(y31n,y);
title（’（h) 谐振器对u（n）的响应y31(n）'）;box on
subplot(2,1,2);y=’y32(n)’；tstem（y32n，y)；
title('(i） 谐振器对正弦信号的响应y32(n)'）；box on
％============================= tstem 程序清单====================
function tstem（xn，yn）
％时域序列绘图函数
％ xn:信号数据序列，yn:绘图信号的纵坐标名称（字符串)
n=0：length(xn)-1;
stem（n,xn，'。’)；box on
xlabel(’n')；ylabel（yn)；
axis（[0，n（end)，min（xn),1。2*max(xn）])

七、实验程序运行结果及分析讨论
程序运行结果如图2.1所示。
实验内容（2)系统的单位冲响应、系统对和的响应序列分别如图（a）、(b)和（c）所示;
实验内容（3）系统h1(n)和h2（n)对的输出响应分别如图（e)和(g)所示;
实验内容(4）系统对和的响应序列分别如图(h）和（i)所示。由图（h)可见,系统对的响应逐渐衰减到零，所以系统稳定。由图(i)可见，系统对的稳态响应近似为正弦序列，这一结论验证了该系统的谐振频率是0。4 rad。

图2。1

实验三 时域采样与频域采样
一、实验目的
(1）掌握模拟信号采样前后频谱的变化
(2）要求掌握频率域采样会引起时域周期化的概念
(3）频率域采样定理及其对频域采样点数选择的指导作用.
二、实验原理与方法
时域采样定理的要点是：
a) 对模拟信号以间隔T进行时域等间隔理想采样，形成的采样信号的频谱是原模拟信号频谱以采样角频率（）为周期进行周期延拓。公式为：

b) 采样频率必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的频谱不产生频谱混叠。
利用计算机计算上式并不方便，下面我们导出另外一个公式，以便用计算机上进行实验。
理想采样信号和模拟信号之间的关系为:

对上式进行傅立叶变换，得到：

在上式的积分号内只有当时，才有非零值，因此：

上式中，在数值上＝,再将代入,得到:

上式的右边就是序列的傅立叶变换，即

上式说明理想采样信号的傅立叶变换可用相应的采样序列的傅立叶变换得到，只要将自变量ω用代替即可。
频域采样定理的要点是:
a) 对信号x（n）的频谱函数X(ejω）在［0，2π]上等间隔采样N点,得到

则N点IDFT[]得到的序列就是原序列x(n）以N为周期进行周期延拓后的主值区序列，公式为:

b) 由上式可知，频域采样点数N必须大于等于时域离散信号的长度M(即N≥M），才能使时域不产生混叠，则N点IDFT[］得到的序列就是原序列x(n），即=x（n）.如果N>M,比原序列尾部多N-M个零点；如果N<M，z则=IDFT［］发生了时域混叠失真，而且的长度N也比x(n)的长度M短，因此.与x(n)不相同。
在数字信号处理的应用中，只要涉及时域或者频域采样,都必须服从这两个采样理论的要点.
对比上面叙述的时域采样原理和频域采样原理,得到一个有用的结论,这两个采样理论具有对偶性:“时域采样频谱周期延拓，频域采样时域信号周期延拓”。因此放在一起进行实验。
二、实验内容及步骤
（1）时域采样理论的验证。
给定模拟信号,
式中A=444。128，=50π，=50πrad/s，它的幅频特性曲线如图3.1

图3。1 的幅频特性曲线
现用DFT(FFT)求该模拟信号的幅频特性,以验证时域采样理论。
安照的幅频特性曲线，选取三种采样频率，即=1kHz,300Hz,200Hz.观测时间选。
为使用DFT，首先用下面公式产生时域离散信号，对三种采样频率，采样序列按顺序用，,表示.

因为采样频率不同，得到的，，的长度不同， 长度(点数）用公式计算。选FFT的变换点数为M=64,序列长度不够64的尾部加零。
X(k)=FFT[x(n)］ , k=0，1,2，3，--——-，M-1
式中k代表的频率为 。
要求: 编写实验程序，计算、和的幅度特性,并绘图显示.观察分析频谱混叠失真。
（2）频域采样理论的验证。
给定信号如下：

编写程序分别对频谱函数在区间上等间隔采样32
和16点,得到:

再分别对进行32点和16点IFFT，得到：

分别画出、的幅度谱，并绘图显示x（n）、的波形，进行对比和分析，验证总结频域采样理论。
提示：频域采样用以下方法容易变程序实现。
① 直接调用MATLAB函数fft计算就得到在的32点频率域采样
② 抽取的偶数点即可得到在的16点频率域采样，即。
当然也可以按照频域采样理论,先将信号x(n)以16为周期进行周期延拓,取其主值区（16点），再对其进行16点DFT（FFT),得到的就是在的16点频率域采样。
四、思考题:
如果序列x（n）的长度为M，希望得到其频谱在上的N点等间隔采样，当N<M时， 如何用一次最少点数的DFT得到该频谱采样？
五、实验报告及要求
a）运行程序打印要求显示的图形，。
b) 分析比较实验结果，简述由实验得到的主要结论
c） 简要回答思考题
d) 附上程序清单和有关曲线。
六、实验程序清单
1 时域采样理论的验证程序清单
% 时域采样理论验证程序exp2a。m
Tp=64/1000; %观察时间Tp=64微秒
%产生M长采样序列x（n)
% Fs=1000;T=1/Fs；
Fs=1000;T=1/Fs;
M=Tp＊Fs;n=0:M—1;
A=444.128；alph=pi＊50*2^0.5;omega=pi*50＊2^0。5;
xnt=A*exp（—alph*n*T)。＊sin（omega＊n＊T）;
Xk=T＊fft（xnt,M）; ％M点FFT［xnt)]
yn=’xa（nT)'；subplot（3，2，1）;
tstem（xnt，yn)； %调用自编绘图函数tstem绘制序列图
box on;title(’(a) Fs=1000Hz’）;
k=0:M—1；fk=k/Tp；
subplot(3，2，2);plot(fk,abs（Xk）);title(’（a） T*FT[xa（nT)］，Fs=1000Hz')；
xlabel（’f（Hz)’)；ylabel（'幅度'）;axis(［0,Fs，0,1.2＊max（abs（Xk)）])
%=================================================
％ Fs=300Hz和 Fs=200Hz的程序与上面Fs=1000Hz完全相同。

2 频域采样理论的验证程序清单
%频域采样理论验证程序exp2b.m
M=27；N=32；n=0:M;
%产生M长三角波序列x(n)
xa=0：floor(M/2）; xb= ceil（M/2）—1：-1:0; xn=[xa,xb];
Xk=fft(xn，1024); ％1024点FFT[x(n)］， 用于近似序列x(n)的TF
X32k=fft（xn，32) ;％32点FFT［x（n)］
x32n=ifft（X32k）； %32点IFFT[X32(k）］得到x32(n)
X16k=X32k(1：2:N)； %隔点抽取X32k得到X16(K)
x16n=ifft（X16k,N/2）； %16点IFFT［X16（k)]得到x16(n）
subplot(3，2，2);stem（n,xn，’。')；box on
title（'（b) 三角波序列x(n)');xlabel（'n’);ylabel('x（n）’)；axis(［0,32，0，20])
k=0:1023;wk=2＊k/1024； %
subplot（3,2,1);plot(wk,abs(Xk));title（'（a）FT[x(n)]’）；
xlabel（’\omega/\pi')；ylabel(’|X(e^j^\omega）|'）；axis(［0,1，0，200]）
k=0:N/2-1；
subplot（3,2，3）；stem（k，abs(X16k）,’。'）；box on
title（'(c) 16点频域采样'）；xlabel(’k');ylabel（'|X_1_6(k）|'）；axis（[0,8,0，200])
n1=0:N/2-1；
subplot（3，2,4）；stem（n1，x16n，’。'）；box on
title('(d) 16点IDFT［X_1_6（k）］’);xlabel（’n')；ylabel（'x_1_6（n)’）；axis（［0,32,0,20]）
k=0：N-1;
subplot（3，2，5）;stem（k,abs（X32k)，’。’);box on
title（’（e) 32点频域采样’);xlabel（'k');ylabel('｜X_3_2（k)|’）;axis([0,16,0，200］)
n1=0:N—1;
subplot(3，2，6)；stem(n1，x32n,’。'）;box on
title('(f) 32点IDFT［X_3_2（k）]'）；xlabel('n'）;ylabel(’x_3_2（n)’)；axis([0，32，0,20］）
七、实验程序运行结果
1 时域采样理论的验证程序运行结果exp2a。m如图3.2所示。由图可见，采样序列的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时,在折叠频率110Hz附近频谱混叠更很严重.

图3.2

2 时域采样理论的验证程序exp2b。m运行结果如图3.3所示。

图3。3
该图验证了频域采样理论和频域采样定理。对信号x（n)的频谱函数X（ejω)在[0，2π］上等间隔采样N=16时, N点IDFT[]得到的序列正是原序列x（n)以16为周期进行周期延拓后的主值区序列:

由于N<M，所以发生了时域混叠失真,因此。与x（n)不相同，如图图3.3（c）和（d)所示。当N=32时，如图图3.3（c)和（d）所示，由于N>M，频域采样定理,所以不存在时域混叠失真，因此与x（n)相同。

实验四 用FFT对信号作频谱分析
一、实验目的
（1）学习用FFT对连续信号和时域离散信号进行谱分析的方法
（2）了解可能出现的分析误差及其原因，以便正确应用FFT。
二、实验原理
用FFT对信号作频谱分析是学习数字信号处理的重要内容.经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是，因此要求.可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱,因此N要适当选择大一些。
周期信号的频谱是离散谱,只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些.
对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号.如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。
三、实验步骤及内容
(1）对以下序列进行谱分析。

选择FFT的变换区间N为8和16 两种情况进行频谱分析.分别打印其幅频特性曲线. 并进行对比、分析和讨论.
(2）对以下周期序列进行谱分析。

选择F