数字信号处理实验教案
指导老师：李春来
指导专业：电子科学与技术
20201111年4月
实验一信号、系统及系统响应
一．实验目的
(1)熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。(2)熟悉时域离散系统的时域特性。
(3)利用卷积方法观察分析系统的时域特性。
(4)掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离
散信号及系统响应进行频域分析。
二．实验原理与方法
采样是连续信号数字处理的第一个关键环节。对一个连续信号x a (t)进行理想采样的过程可用（1）式表示。其中(t)为x a (t)的理想采样，p(t)为周期冲激脉冲，即
?()()*()
a x t p t x t =（1）
其傅里叶变换(j Ω)为
（2）
或将(1)式进行傅里叶变换，
式中的x a (nT)就是采样后得到的序列x(n)，即
x(n)的傅里叶变换为：
()()
n p t t nT δ∞
=?∞
=
?∑
^
x ^
a X ?1()[()]
a a s m X j X j m T ∞=?∞
?=???∑^
()[()()]()()()j t a a n j t a n j nT
a n X j x t t nT e dt
x t t nT e dt
x nT e δδ∞
∞
???∞
=?∞
∞
∞???∞
=?∞
∞
??=?∞
?=?=?=
∑∫
∑∫∑
()()
a x n x nT =()()j j n
n X e x n e ω
ω∞
?=?∞
=
∑
可知
在数字计算机上观察分析各种序列的频域特性，通常对X(e j ω)在［0,2π］上进行M 点采样来观察分析。对长度为N 的有限长序列x(n),有
一个时域离散线性非移变系统的输入/输出关系为
上述卷积运算也可以在频域实现
三．实验内容及步骤
(1)认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内
容，阅读本实验原理与方法。(2)编制实验用主程序及相应子程序。
①信号产生子程序，用于产生实验中要用到的下列信号序列：
x a (t)=Ae -at sin(Ω0t)u(t)
进行采样，可得到采样序列
x a (n)=x a (nT)=Ae -anT sin(Ω0nT)u(n),
0≤n<50
其中A 为幅度因子，a 为衰减因子，Ω0是模拟角频率，T 为采样间隔。这些参数都要在实验过程中由键盘输入，产生不同的x a (t)和x a (n)。
b.单位脉冲序列：x b (n)=δ(n)
c.矩形序列：x c (n)=R N (n),N=10
②系统单位脉冲响应序列产生子程序。本实验要用到两种FIR 系统。
a.h a (n)=R 10(n);
b.h b (n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)
③有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。可以直接
^
()()
j a T
X j X e ωω=??=10
()()2,0,1,,1k N j n
j k
n k X e
x m e k k M M
ωωπ
ω??===
=????∑()()()()()
m y n x n h n x m h n m ∞
=?∞
=?=
?∑
()()()
j j j Y e X e H e ωωω=
调用MATLAB 语言中的卷积函数conv 。conv 用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。调用格式如下：
y=conv (x,h)
图1实验一的主程序框图
图2x a (t)的幅频特性曲线
(3)调通并运行实验程序，完成下述实验内容：
①分析采样序列的特性。
a.取采样频率f s =1kHz,即T=1ms 。
b.改变采样频率，f s =300Hz ，观察|X(e j ω)|的变化，并做记录(打印曲线)；进一步降低采样频率，f s =200Hz ，观察频谱混叠是否明显存在，说明原因，并记录(打印)这时的|X(e j ω)|曲线。
10.8
0.60.40.200
100
200300
400
500
x a (j f )
f /H z
②时域离散信号、系统和系统响应分析。
a.观察信号x b (n)和系统h b (n)的时域和频域特性；利用线性卷积求信号x b (n)通过系统h b (n)的响应y(n)，比较所求响应y(n)和h b (n)的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。
b.观察系统h a (n)对信号x c (n)的响应特性。③卷积定理的验证。
四．思考题
(1)在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里
叶变换频谱的数字频率度量是否都相同?它们所对应的模拟频率是否相同?为什么?(2)在卷积定理验证的实验中，如果选用不同的频域采样点数M 值，例如，选M=10
和M=20，分别做序列的傅里叶变换，求得
所得结果之间有无差异?为什么?
五．实验报告要求
(1)简述实验目的及实验原理。
(2)按实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时域
和幅频特性曲线，并对所得结果进行分析和解释。(3)总结实验中的主要结论。(4)简要回答思考题。
实验二快速傅立叶变换FFT 及信号的谱分析
()()(),0,1,,1
k k k j j j a Y e X e H e k M ωωω==????
一．实验目的
1、学习时间抽取快速傅立叶变换DIT—FFT 算法。
2、掌握DIT—FFT 算法的计算机实现方法。
3、研究如何利用FFT 程序分析确定性时间连续信号的频谱。
4、用DFT/FFT 对信号进行频谱分析的误差研究
二．实验原理
1、快速傅立叶变换FFT 和快速傅立叶反变换IFFT
长度为N 的序列)(n x 的离散傅立叶变换)(k X 为：
X k x n W k N N nk n N ()(),,....,==?=?∑01
1
首先按n 的奇偶把时间序列x(n)分解为两个长为N/2点的序列
x n x r 12()()=r=0,1,...,N/2-1x n x r 221()()
=+r=0,1,...,N/2-1
则x(n)的DFT
X k ()为
X k x n W x r W x r W x r W
x r W W n N N
kn
r
N N
kr r
N N k r r
N N
kr r
N N kr
N k
()()()()()()//()
//==++=
+
=?=?=?+=?=?∑∑∑∑
∑01
21
2021
210
21
120
21
22221由于
W
e
e
W N
kr j
N
Kr j
N kr N kr 22222
2===??ππ//，故有X k x r W W
x r W X k W X k k N r
N N kr
N k r
N N kr
N k ()()()()()
,,...,/////=
+=+=?=?=?∑∑0
21
12
21
22
120121
其中X k X k 12()()和分别为x n x n 12()()和的N/2点DFT。因为X k X k 12()()和均是以N/2为周期的，且W W N
k N N k +=?/2
。因此可将N 点DFT )(k X 分解为下面的形式
X k X k W X k N k
()()()
=+12k=0,1,...,N/2-1X k N
X k W X k N k ()()()+
=?2
12k=0,1,...,N/2-1
通过上面的推导可以看出，N 点的DFT 可以分解为两个N/2点的DFT，每个N/2点的DFT 又可以分解为两个N/4点的DFT。依此类推，当N 为2的整数次幂时(M
N 2=)，由于每分
解一次降低一阶幂次，所以通过M 次的分解，最后全部成为一系列2点DFT 运算。以上就是
按时间抽取的快速傅立叶变换(FFT)算法。
序列)(k X 的离散傅立叶反变换为
x n N
X k W
n N N
nk k N ()(),
,....,=
=??=?∑1
01
1
离散傅立叶反变换与正变换的区别在于N W 变为1
?N W ，并多了一个N 1的运算。因为N W 和1?N W 对于推导按时间抽取的快速傅立叶变换算法并无实质性区别，因此可将FFT 和快速傅立叶反变换（IFFT）算法合并在同一个程序中。
2、用基2算法的F F T 程序计算任意长度的有限长序列的频谱的方法
频率分辨率可以从两个方面来定义：一是用某种算法（如谱分析、功率谱估计等）将原始信号)(n x 中两个靠得很近的谱峰仍然能保持分开的能力，二是在使用FFT 时，在频率轴上能得到的最小频率间隔f ?。
长为N 的序列f(n)的频谱可由它的付立叶变换给出:
()()n
j N n j e
n f e
F ωω
??=∑=1
其中ω是在区间0≤ω≤2π上连续取值的。
在用数字机对序列作谱分析时,不可能沿ω作连续的计算。因此只好按频率分辨力Δω的要求将频段0≤ω≤2πM 等分M≥2π／Δω,各分点的座标为:
M
K K /2πω=K=0,1,2,…,M-1
再计算各分点处的谱系数()ω
j e
F ,作为对序列f(n)的谱分析结果:
()()M
Kn j N n j e
n f e
F K
/21
πω??=∑=K=0,1,2,…,M-1
当M≥N 时,此式即为N 点序列f(n)的M 点的DFT 定义式。在采用基2算法的F F T 程序计算此式时,程序对参数的基本要求是：M=N=m 2,其中m 为正整数。为保证此要求,须:
(1)选M≥2π／Δω是最接近2π／Δω的2的整次幂
(2)当M N <时,须给)(n f 补上N M ?个零值点,即构造一个长为M 的新序列
()n f 1：
()()??
??≤≤?≤≤=1
0101M n N N n n f n f 此M 点序列f1(n)的M 点的DFT 结果就是原N 点序列f(n)的M 点的DFT 结果。而f1(n)的M 点DFT 是可直接调用基2的FFT 程序作快速计算的,这样就解决了实际点数N 不是2的整次幂时采用基2算法的困难。
下面给出做FFT 时参数选择的一般原则：
a)若已知信号的最高频率c f ，为防止混叠，应保证抽样频率s f 满足c s f f 2≥。b)根据实际需要，选定频率分辨率f ?，以确定FFT 的点数N：f f N S ?≥c)
s f 和N 确定后，可计算对于模拟信号的最小记录时间p T ：s
p f N T =d)每条谱线所代表的频率刻度k f ：2
,1,0,N k N
f k
f s
k ?
==三．实验内容：
实验中用到的信号序列：a)Gaussian 序列
b)衰减正弦序列
c)三角波序列
d)反三角波序列
上机实验内容：
1、观察高斯序列的时域和幅频特性，固定信号x a (n)中参数p=8，改变q 的值，使q 分别等于2，4，8，观察它们的时域和幅频特性，了解当q 取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p 分别等于8，13，14，观察参数p 变化对信号序列的时域及幅频特性的影响，观察p 等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。
2、观察衰减正弦序列x b (n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f 分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。
3、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT 分析信号序列x c (n)和x d (n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。
在x c (n)和x d (n)末尾补零，用N=16点FFT 分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT 频谱还有相同之处吗？这些变化说明了什么？
4、编写基2DIT—FFT 子程序，代替MATLAB 中命令fft,观察两者异同。
四．思考题
1、根据实验中各)(n x 的X(k)值以及频谱图，说明参数的变化对信号频谱产生哪些影
响？
2、基2FFT相对于DFT在运算速度上有什么改进？
五．实验报告要求
1、简述实验目的及实验原理。
2、整理好经过运行并证明是正确的程序，并且加上详细的注释。
3、N＝1024点DIT—FFT分析所用时间与直接计算DFT所用时间对比，理解“快”
字。
实验三IIR 滤波器的设计
一.实验目的
1、掌握用脉冲响应不变法，双线性变换法设计IIR DF 的原理和具体的设计方法。
2、掌握数字滤波器的计算机仿真。
3、了解模拟Butterworth，Chebyshev,Cauer 滤波器的计算机实现。
二.实验原理和方法
IIR 数字滤波器的系统函数为z ?1的有理分式
1
1
1
01)(?=?=∑∑+=
z a z b z H k N
k k N
k
设计IIR 滤波器的系统函数，就是要确定H(z)的阶数N 以及分子分母多项式的系数a k 和
b k ，使其ω
ωj e z j z H e H ==)
()(满足指定的频率特性
由于模拟滤波器的设计有许多简单而严谨的设计公式和大量的图表可以利用，因此IIR 滤波器设计的方法之一是，先设计一个合适的模拟滤波器，然后将模拟滤波器通过适当的变换转换成满足给定指标的数字滤波器。本次实验设计的变换方法是双线性变换法。
双线性变换法是从频域响应出发，直接使数字滤波器的频域响应H e j ()ω逼近模拟滤波器的频域响应|H j ()?|,进而求出H (z )。
双线性变换法通过两次映射，避免了脉冲响应不变法的频谱混叠：首先将S 平面压缩到某一中介的S 1平面的一条横带域（-
πT ，π
T
）内，然后再通过标准的变换将此横带域映射到整个Z 平面上去，这样就使S 平面与Z 平面建立了一一对应的关系，消除了多值性。
??=t g T
(
)12：实现S 平面的jΩ轴压缩到S 1平面jΩ1轴上的(-j πT ，j πT
)段上。T s e z 1=：将Ω1映射到Z 平面的单位圆上。
得到S 平面到Z 平面的映射关系：
1
1112??+?=
z z T s 或
s
T s T z )2/(1)2/(1?+=
模拟频率与数字频率的关系为
?=t g (
)ω2
这是一种非线性的关系，这种非线性关系使得模拟滤波器与数字滤波器的响应与对应的频率关系上发生了畸变，也造成了相位的非线性变化。为保证各边界频率点为预先指定的频率，例如低通滤波器的截止频率ωc 与阻带边频ωs ，在确定模拟低通滤波器系统函数之前必须按下式进行频率预畸变：
?
c
c
t g =(
)ω2
用双线性变换法设计IIR DF 的步骤如下：
1、确定DF 的通带频率、阻带频率，通带最大衰减和阻带最小衰减。
2、计算频率预畸变后对应的模拟低通滤波器的频率。
3、确定模拟低通滤波器的阶数N 和3dB 截止频率c ?。
4、确定模拟低通滤波器的系统函数H(s)。
5、由H(s)经过反归一化、双线性变换确定数字低通滤波器的系统函数H(s)。
6、设计其它形式的滤波器时，由模拟低通到所需类型滤波器的频率域变换直接得到。
三.实验内容及步骤：
1、编程实现模拟Butterworth、Chebyshev 滤波器。
用脉冲响应不变法和双线性变换法设计下列满足下列指标Butterworth 模拟低通滤
波器:fsam=1Hz ；fp=0.2Hz;fs=0.3Hz Ap=1dB;As=25dB
可能用到的函数：
[N,wc]=buttord(wp,ws,Ap,As,'s')[numa,dena]=butter(N,wc,'s');
[numd,dend]=impinvar(numa,dena,fsam);[numd,dend]=bilinear(numa,dena,fsam);2、用双线性变换法设计下列IIR 滤波器:
1）、设计一Chebyshv1型高通IIR 滤波器：
fsam=1Hz ；fp=0.2Hz;fs=0.3Hz ；Ap=0.8dB;As=20dB 2）、设计一Butterworth 型带通IIR 滤波器：
fpl=20Hz;fph=30Hz;fsl=15Hz;fsh=35HzAp=1dB;As=40dB;fsam=100Hz 可能用到的函数：
[N,wc]=cheb1ord(wp,ws,Ap,As,'s');[num,den]=cheby1(N,Ap,wc,'s');[numt,dent]=lp2hp(num,den,1);
[numt,dent]=lp2bp(numa,dena,w0,B);[numd,dend]=bilinear(numt,dent,fsam)
3、设计一椭圆IIR 低通滤波器，求出其传输函数，并画出其增益响应曲线：
fp=80Hz;fs=1000Hz;fsam=4000Hz ；Ap=0.5dB;As=40dB
可能用到的函数：
[N,wc]=ellipord(wp,ws,Ap,As,'s');[num,den]=ellip(N,Ap,As,wc,'s');
[numd,dend]=impinvar(num,den,fsam);[numd,dend]=bilinear(num,den,fsam)
4、若信号由5Hz,15Hz,30Hz 三个正弦频率成分构成，设计一个椭圆滤波器滤除5Hz,30Hz 频率成分
提示:N=4阶，Rp=0.1dB,Rs=40dB
四.思考题
用双线性变换法设计数字滤波器过程中，变换公式
1
1112??+?=
z z T s 中T 的取值对设计结果有无影响？为什么？
五.实验报告要求
1、简述实验目的和原理。
2、整理好经过运行并证明是正确的程序并且加上详细的注释。
3、给出程序运行结果，并对结果加以分析。
4、简述双线性变换法和脉冲响应不变法的异同。
5、简要回答问题。
实验四FIR 滤波器的设计
一.实验目的
1、掌握用窗函数法设计FIR DF 的原理和方法。
2、熟悉线性相位FIR DF 的特性。
3、了解各种窗函数对滤波特性的影响。
4、掌握等波纹线性相位FIR DF 设计方法
二.实验原理和方法
1、线性相位FIR DF 的特性
FIR 数字滤波器相对于IIR 数字滤波器的最大优点是能够做到严格线性相位。当FIR 滤波器具有严格的线性相位时，系统的单位脉冲响应h(n)满足下述条件：
)
1()(n N h n h ??±=针对单位脉冲响应h(n)满足奇对称或偶对称，FIR 滤波器的h(n)长度为奇数或偶数，线性相位FIR 滤波器的幅频特性共分成四种情况：
a、h n h N n ()()=??1，N 取奇数：适合于设计任何在ωππ=02,,是偶对称频率特性的滤波器。
b、h n h N n ()()=??1，N 取偶数:不适合作在ωπ=处不等于零的滤波器，例如高通、带阻滤波器：
c、h n h N n ()()=???1，N 为奇数：不适合作在ωππ=02,,处为偶对称的滤波器，例如低通和高通滤波器
d、h n h N n ()()=???1，N 为偶数：不适合作在ωπ=02,处为偶对称的滤波器，例如低通滤波器
2、窗函数法设计FIR 滤波器
设欲设计滤波器的理想频率响应为H e d j ()ω，单位脉冲响应为h n d ()，h n d ()与H e d j ()ω是一对富氏变换，因此有
H e
h n e d j n
d j n ()()ω
ω
=
=?∞
∞
?
∑h n H e e d d d j j n ()()=
?
∫1
2π
ω
π
π
ωω根据给定的H e d j ()ω求得的h n d ()一般是无限长的且是非因果的。为了得到一个因果的有限长的滤波器h n ()，最直接的方法是截断h n d ()，或者说用一个窗口函数w n ()对h n d ()进行加窗处理
h n ()
=h n d ()w n ()
h n ()成为实际设计FIR 滤波器的
单位脉冲响应,其频率响应为)(ω
j e
H 为：
∑?==10
)()(N n j j e n h e H ω
ω
其中N 为窗口w n ()的长度。窗口函数的形状和窗口长度N 决定了窗函数法设计出的FIR 滤波器的性能。设计时，要根据阻带的最小衰减和过渡带宽度来选择恰当的窗函数类型和窗口长度N。各种窗函数能够达到的阻带的最小衰减和过渡带宽可参见课本的相关章节，以编写程序时使用。1、设窗函数法设计FIR 滤波器的步骤：（1）给定频响函数H e d j ()ω；（2）求出单位抽样响应h n ()；
（3）根据过渡带宽度和阻带最小衰减，确定窗函数类型和窗口长度N；（4）求实际设计FIR 滤波器的h n ()=h n d ()w n ()以及)(ω
j e H ；
（5）验证)(ωj e H 是否满足要求；
此外，如果要求线性相位，h(n)还应该满足对称条件：
)
1()(n N h n h ??±=根据设计滤波器的滤波特性，正确选择一种类型，例如，要设计线性相位低通滤波器时，可选择偶对称一类h n h N n ()()=??1，而不是奇对称一类h n h N n ()()=???1。
三.实验内容：
1复习窗函数法设计FIR 滤波器的内容，阅读实验原理，掌握设计步骤。2编写程序：
a)编写矩形窗、三角窗、海宁窗、汉宁窗和凯塞窗的窗函数子程序。b)编写主程序实现FIR DF。
3N=23,4πω=c ，用上述五种窗口设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察3dB 和20dB 带宽以及阻带最小衰减，比较五种窗口对滤波器性能的影响。
4确定一个等波纹线性相位FIR LPF 的传输函数，并画出其增益响应，指标如下：通带截止频率Fp=800Hz,阻带截止频率Fs=1000Hz,通带波纹Rp=0.5dB,最小阻带衰减Rs=40dB，抽样频率Ft=4000Hz.
5若信号由5Hz,15Hz,30Hz 三个正弦频率成分构成，设计一个等波纹线性相位FIR 滤波器滤除5Hz,30Hz 频率成分
四.思考题
1、窗口的长度对滤波器的性能有什么影响？
2、吉布斯效应是怎样产生的？增加窗口长度能否消除吉布斯效应，为什么？
五.实验报告要求
1、简述实验目的和原理。
2、整理好经过运行并证明是正确的程序并且加上详细的注释。
3、给出程序运行结果，并对结果加以分析。
4、分析各种窗口类型对滤波器性能的影响。
5、简要回答问题。
实验参考源程序：
实验一、信号、系统及系统响应
%高斯序列
clear;
n=0:50;
A=444.128;a=50*sqrt(2)*pi;
T=0.001;
%T=1/300(1/200);
w0=50*sqrt(2)*pi;
x=A*exp(-a*n*T).*sin(w0*n*T); close all;
subplot(3,1,1);stem(x);
title('理想采样序列');
k=-50:50;
w=(pi/25)*k;
X=x*(exp(-j*pi/25)).^(n'*k); magX=abs(X);
subplot(3,1,2);stem(magX);
title('理想采样序列的幅度谱'); angX=angle(X);
subplot(3,1,3);stem(angX);
title('理想采样序列的相位谱');
%单位脉冲序列;
clear;
n=1:50;
x=zeros(1,50);
x（1）=1;close all;
subplot(3,1,1);stem(x);
title('单位脉冲序列');
k=-50:50;
X=x*(exp(-j*pi/25)).^(n'*k); magX=abs(X);
subplot(3,1,2);stem(magX);
title('单位脉冲序列的幅度谱'); angX=angle(X);
subplot(3,1,3);stem(angX);
title('单位脉冲序列相位谱');
%矩形序列
clear;n=1:5;
x=sign(sign(10-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);
subplot(3,1,2);stem(magX);
title('矩形序列的幅度谱');
angX=angle(X);
subplot(3,1,3);stem(angX);
title('矩形序列相位谱');
%冲击序列
close;
n=1:50;
x=zeros(1,50);
x(1)=1;x(2)=2.5;x(3)=2.5;x(4)=1;
close all;
subplot(3,1,1);stem(x);title('单位冲击信号序列')
k=-50:50;X=x*(exp(-j*pi/25)).^(n'*k); magX=abs(X);
subplot(3,1,2);stem(magX);
title('单位冲击信号的幅度谱')
angX=angle(X);
subplot(3,1,3);stem(3,1,3);
title('单位冲击信号的相位谱');
%卷积计算及其定律验算
clear;
close all;
n=1:50;
hb=zeros(1,50);
hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1; figure(1);
subplot(3,1,1);stem(hb);
title('系统响应hb(n)');
m=1:50;
T=0.001;
A=444.128;
a=50*sqrt(2.0)*pi;
w0=50*sqrt(2.0)*pi;
x=A*exp(-a*m*T).*sin(w0*m*T); subplot(3,1,2);stem(x);
title('输入信号x(n)');
y=conv(x,hb);
subplot(3,1,3);stem(y);
title('输出信号y(n)');
k=-50:50;
X=x*(exp(-j*pi/25)).^(n'*k); figure(2);
magX=abs(X);
subplot(3,2,1);stem(magX);
title('输入信号的幅度谱'); angX=angle(X);
subplot(3,2,2);stem(angX);
title('输入信号的相位谱');
Hb=hb*(exp(-j*pi/25)).^(n'*k); magHb=abs(Hb);subplot(3,2,3);stem(magHb);
title('系统响应的幅度谱');
angHb=angle(Hb);
subplot(3,2,4);stem(angHb);
title('系统响应的相位谱');
n=1:99;k=1:99;
Y=y*(exp(-j*pi/25)).^(n'*k);
magY=abs(Y);
subplot(3,2,5);stem(magY);
title('输出信号的幅度谱');
angY=angle(Y);
subplot(3,2,6);stem(angY);
title('输出信号的相位谱');
XHb=X.*Hb;figure(3);
subplot(2,1,1);stem(abs(XHb));
title('输入信号的幅度谱与系统响应的幅度谱相乘');
subplot(2,1,2);stem(abs(Y));
title('y(n)的幅度谱');
axis([0,60,0,8000]);
实验二快速傅立叶变换FFT及信号的谱分析
function x=revorder(xn)%%时域序列抽取
N=length(xn);
LH=N/2;
j1=LH;
N2=N-2;
for I=1:N2
if(I t=xn(I+1); xn(I+1)=xn(j1+1); xn(j1+1)=t; end k=LH; while(j1>=k) j1=j1-k; k=k/2; end j1=j1+k end x=xn function ffts xi=input('xi=');%%[1,5,3] Ni=length(xi); N=input('N=');%%8 M=log2(N); xn=[xi,zeros(1,N-Ni)]; %x=revorder(xn); x=bitrevorder(xn); for L=1:M B=pow2(L-1); for J=0:B-1 P=pow2(J,M-L); for k=J+1:(2^L):N X(k)=x(k)+x(k+B)*(exp(-j*2*pi/N*P)); X(k+B)=x(k)-x(k+B)*(exp(-j*2*pi/N*P)); x(k)=X(k); x(k+B)=X(k+B); end end end figure; subplot(311);stem(xi);title('xi'); subplot(312);stem(abs(x));title('abs(x)'); subplot(313);stem(abs(fft(xn)));title('abs(fft(xn))'); %高斯序列 clear; n=0:15; p=8;q=2;x=exp(-1*(n-p).^2/q); %close all; subplot(321);stem(abs(fft(x)));%%%%%%%%fft(x,N)%%fft(x,64) subplot(322);stem(x);title('p=8,q=2'); p=8;q=4;x=exp(-1*(n-p).^2/q); subplot(323);stem(abs(fft(x))); subplot(324);stem(x);title('p=8,q=4'); p=8;q=8;x=exp(-1*(n-p).^2/q); subplot(325);stem(abs(fft(x))); subplot(326);stem(x);title('p=8,q=8'); clear; n=0:15; p=8;q=8;x=exp(-1*(n-p).^2/q); figure(2); subplot(321);stem(abs(fft(x)));%%%%%%%%fft(x,N)%%fft(x,64) subplot(322);stem(x);title('p=8,q=8'); p=13;q=8;x=exp(-1*(n-p).^2/q); subplot(323);stem(abs(fft(x))); subplot(324);stem(x);title('p=13,q=8'); p=14;q=8;x=exp(-1*(n-p).^2/q); subplot(325);stem(abs(fft(x))); subplot(326);stem(x);title('p=14,q=8'); %正弦衰减序列 clear; n=0:15; a=0.1;f=0.0625;x=sin(2*pi*f*n); figure(3); subplot(321);stem(abs(fft(x)));%%%%%%%%fft(x,N)%%fft(x,64) subplot(322);stem(x);title('a=0.1,f=0.0625'); a=0.1;f=0.4375;x=exp(-a*n).*sin(2*pi*f*n); subplot(323);stem(abs(fft(x)));%%%%%%%%fft(x,N)%%fft(x,64) subplot(324);stem(x);title('a=0.1,f=0.4375'); a=0.1;f=0.5625;x=exp(-a*n).*sin(2*pi*f*n); subplot(325);stem(abs(fft(x)));%%%%%%%%fft(x,N)%%fft(x,64) subplot(326);stem(x);title('a=0.1,f=0.5625'); %三角波和反三角波 clear for i=1:4 x(i)=i; end for i=5:8 x(i)=9-i; end figure(4);subplot(2,2,1);stem(x); subplot(222);stem(abs(fft(x,16)));%%%%%%%%fft(x,N)%%fft(x,64) for i=1:4 x(i)=5-i; end for i=5:8 x(i)=i-4; end subplot(223);stem(x); subplot(224);stem(abs(fft(x,16)));%%%%%%%%fft(x,N)%%fft(x,64)
相关主题
文本预览