数字信号处理实验报告1

对时域进行了采样， 所以相当于对模拟域进行了以 f s 1 / T 1600Hz 为周期的周期延拓， 并以 f s 2 变换到数字域。又因为仅取了 [0, 2 ) 这一段的值，所以频域显示出在 2 / 32 处有值。 作为验证，利用 MATLAB 绘制了 DTFT 变换的时频域图：
e) 信号频率 F＝50Hz，采样点数 N=64 ，采样间隔 T=0.000625s ；
相较于(a)问，本问中采样点数变为 64 点，采样频率 f s 1600Hz 。在数字域， sin (100 t )信号为频 谱在 / 16 上的冲激,并以 2 为周期做周期 延拓。矩形窗的长度变为 64，频谱中 sinc 函数主瓣宽度为 / 32 ，在频域这二者的卷积，相当于把 sinc 函数搬移到 2k 的位置上，再在 2 的周期内对频域进行 64 点采样，得到 X (k ) 。 本问题中，对正弦波采了 2 个周期，每个周期采 32 个点。即矩形窗对周期信号恰周期截断，所得到的 频谱并没有发生泄漏，很好的反映了正弦波的单频频谱特性。
可以看出推导出的频率 6.25 Hz 与真实频率 250Hz 相差很大。 通过分析可以发现一系列问题。 采样频率 f s 200Hz ， 信号频率为 250Hz .则在数字域，sin(500 t ) 信号为频谱在 5 / 2 上的冲激,并以 2 为周期做周期延拓。 矩形窗的长度变为了 32， 频谱中 sinc 函数主 瓣宽度为 /16 ，在频域这二者的卷积，相当于把 sinc 函数搬移到 (5 / 2) 2 k 的位置上，再在 2 的周期内对频域进行 32 点采样，得到 X ( k ) 。 本问题中，
k
Hale Waihona Puke Baidu
X (k ) X 1 (k ) WN k X 2 (k ), k X ( N / 2 k ) X 1 (k ) WN X 2 ( k ),
用蝶形图表示为：
k 0,1,..., N/ 2 1 k 0,1,..., N/ 2 1
只要 N 是 2 的整数次幂，这种分解就可一直进行下去，将 DFT 运算转换为.. 级运算，每一级运 算都由 N / 2 个碟形运算组成。利用长度为 N 的数组存放输入数据和每一级计算的结果。
f) 信号频率 F＝250Hz，采样点数 N=32 ，采样间隔 T=0.005s ；
本问中，看似所得到的频谱并没有破坏正弦波的单频特性，但是根据(a)问中，由频域推导时域信号的 频率：


T

1 1 1 1 2 2 =2 6.25 TN 0.005 32
g) 将 c)中信号后补 32 个 0 ，做 64 点 FFT ，并与直接采样 64 个点做 FFT 的 结果进行对比。
本问中实际上是对正弦信号进行了 32 点采样之后，再补零至 64 点，与长度为 64 的矩形窗卷积所 得。补零实际上，系统的 z 变换并没有发生改变，相当于在 z 平面单位圆周上采了更多的点，使得其更光 滑，更接近系统本身的 DTFT，减轻栅栏效应。 因此，本小问实际上是在 c)中卷积所得到的频谱上（ ）内借去了更多的采样点。 在数字域对应的是 的冲激， 并不是 sinc 函数主瓣宽度 的整数倍， 在 的位置有最大值。 此位置恰好是 的 15 倍， 因此在 K=15,49 上有峰值，其余奇数点上由于主瓣宽度的影响为 0，而偶数点上相当于 c)中每隔的采样点，值不变。
(2).原位计算 对 N 2M 点的 FFT 共进行 M 级运算，每级由 N / 2 个蝶形运算组成。在同一级中，每个蝶形的
输入数据只对本蝶有用， 且输出节点与输入节点在同一水平线上， 这就意味着每算完一个蝶形块后， 所得数据可立即存入原输入数据所占用的数组元素(存储单元)，这种原位(址)计算的方法可节省大 量内存。
由 / f s Ts ，可以由频谱图得到原始正弦信号的周期为：

同时，

T

1 1 1 1 2 2 =2 50 TN 0.000625 32
N T 32 0.000625 1 ，相当于对正弦波采样了一个周期，恰好形成周期截断，频谱并没有 1/ f s 1/ 50
（2M L 1) 2 L R N 2 L R
本级中各蝶的第二个节点与第一个节点都相距 B 点。
应用原位计算，蝶形运算可表示成如下形式：
P AL (k ) AL 1 (k ) WN AL -1 (k B) P A L (k B) AL 1 (k ) WN AL 1 (k B)
32 0.004 =6.4 1/ 50
sin(100 t ) 0 t NT x(t ) 0 t NT
图1
x(t )
的频谱图
的频谱图， 如图 1 所示，可以看出是此题所求频谱的包络图。 又令 N 35 ,则此时 期截断，没有频谱泄漏。
35 0.004 =7 ，是整数周期，可以看出与问题(a),(b)具有同样的形式，形成周 1/ 50
(3). 蝶形运算
蝶形运算是分级进行的； 每级的蝶形运算可以按旋转因子的指数大小排序进行； 如果指数大小一样则 可从上往下依次计算。对 N 2M 点的 FFT 共有
L 1, 2, ... M。 第 L 级共有 B 2
L 1
M 级运算，用 L 表示从左到右的运 算 级 数
个不同指数的旋转因子， 用 R 表示这些不同指数旋转因子从上到
可以看出，采样点除了在 1 时采样到了最高点外，其它点都恰好采样在了零值点上。 另一方面，信号相当于一个连续正弦信号乘以一个长度为 32 的矩形窗，则频谱中 sinc 函数的主瓣宽 度为 2
（2 / 32） +2k 处。 / 32 ， 在频域对两者进行卷积后相当于把 sinc 函数搬移到 再对频域进行采样。
N T 32 0.005 8 为整数，所以恰好周期截断，频谱未发生泄漏，但一个周期内采样的 1 / fs 1 / 50
密集度低于 a)中所示，较好的反映了正弦信号的单频特性。
c) 信号频率 F＝50Hz，采样点数 N=32 ，采样间隔 T=0.0046875s ；
本问中，由得到的 X ( k ) 就可以明显看出频谱的采样发生了严重的泄漏，并没能反映出正弦信号的单频 特性。 这是因为由于频域的采样是时域的周期延拓，从时域图中我们可以看出，当对这 32 个点进行周期延拓 时，信号的包络已经不是正弦波形了，而在(d)中做了一些工作探究这个现象。
数字信号处理实验
快速傅里叶变换与信号频谱分析
一． 实验目的
1. 在理论学习的基础上，通过本实验加深对离散傅里叶变换的理解。 2. 熟悉并掌握按时间抽取编写快速傅里叶变换（FFT）算法的程序。 3. 了解应用 FFT 进行信号频谱分析过程中可能出现的问题，例如频谱混淆
二．实验原理，内容与分析
1. 仔细分析教材第六章“时间抽取法 FFT 的 FORTRAN 程序” ，编写出相应的使用 FFT 进行信 号频谱分析的 Matlab 程序。
(1) 时间抽取 FFT 算法
有限长序列 x(n)的 N 点 DFT 定义为：
kn X (k ) x(n)WN n 0
2 N
N 1
式中 WN e
j
。
若对 x(n) 做一次 DFT 运算，时间复杂度为 N 2 级，其运算效率低下。而利用基于时间（或频率） 抽取 FFT 方法，时间复杂度可降至 N log N 级，大大提高了效率。 基于时间抽取的 FFT 方法将序列 x(n) 分解为两组：偶数和技术，利用 WN 的周期性和对称性， 可将 X ( k ) 表达为前后两部分：
32 0.0046875 另一方面， 这是由于加了时窗之后， 即相当于对 7.5 个周期的正弦波进行了截断， 7.5 ， 1 / 50
每周期 4.266 个采样点。非整数周期的截断导致了频谱的泄漏，此时所做的 FFT 运算实际是没有工程意义 的，因为已经破坏了原信号的频谱特性。 本问中， sin 100 t 在数字域对应的是 0.468 的冲激，并不是 sinc 函数主瓣宽度 / 16 的整数倍， 在 0.468 2k 的位置有最大值。 X (k ) 是对频谱的采样，由此可知：当时窗非整数周期截断时会造成原 信号频谱的泄漏。
d) 信号频率 F＝50Hz，采样点数 N=32 ，采样间隔 T=0.004s ；
本问中，频谱仍然发生了泄漏。原因同上一个问题相同， 即相当于对 6.4 个周期的正弦波进行了截断，每周期 5 个采 样点。 同样是非整数周期的截断导致了频谱的泄漏， 此时所做 的 FFT 运算因为已经破坏了原信号的频谱特性实际是没有工 程意义的。 但是与上一问比较发现，周期内采样点数的增加 使得频谱的变换更为平缓，分辨率适当提高。 同样在本问中， sin(100 t ) 在数字域对应的是 0.4 的 冲激，并不是 sinc 函数主瓣宽度的整数倍，在 0.4 2k 的位置有最大值。 X ( k ) 是对频谱的采样。 如上一问题(c)题中所说，本题中，为了验证以上结论，使 用 MATLAB 画出了模拟信号：
2.
用 FFT 程序分析正弦信号，分别在以下情况进行分析，并讨论所得的结果：
a) 信号频率 F＝50Hz，采样点数 N=32 ，采样间隔 T=0.000625s ；
1 j 2 ft j 2 ft e ) ， 本题中，sin(2 ft ) 的频谱是离散的且由于 sin 2 ft 2 j (e 仅在 2 50 处有值。 由于
ii
c)与 d) 比较 c)与 d)，虽然同为非整数周期的截断，但是由于 d)中每个周期采样点数的增加，使得所得到的
实验结果对比分析
i
a)与 c)
这两个图所得到频谱的不一致，仅仅在于采样频率的不一致。在 a)问中，对正弦波恰好周期截断，共采了 一个周期，采了 32 个点。而在 c)中，采了 7.5 个周期，每个周期 4.266 个点。因此我们发现，矩形窗对 正弦波非整数周期的截断，即非整数点非整数周期的采样严重破坏了正弦波自身的单频特性，将会导致频 谱发生泄漏。
32 0.005 40 ，即相当于对正弦波采了 40 个周期，每个周期采 0.8 个点。由奈奎斯采 1/ 250
样定理，采样频率 f s 2 f ，而在本问中，已经不满足奈奎斯定理的要求，频谱已经发生了混淆。混淆导 致在 / 2,3 / 2 处也产生了冲激信号，实际上参与卷积的是这两个信号，并不是正弦信号的频谱信号，而 是混淆之后的。虽然此时矩形窗对周期信号恰周期截断，所得到的频谱并没有发生泄漏，但是已发生了混 淆，得到的频谱也是没有价值的。
发生泄漏，很好的反映了正弦波的频谱特性。
b) 信号频率 F＝50Hz，采样点数 N=32 ，采样间隔 T=0.005s ；
同上问题类似，采样频率. f s 1/ T 200 Hz .，频谱在 / 2 有值，并以 2 为周期做周期延拓。宽 度为 32 的矩形窗，对应的频谱中 sinc 函数主瓣宽度为 / 16 .。 并且由于
下的顺序 R 0,1,..., B 1 。第 R 个旋转因子的指数 P 2M L R ，旋转因子指数为 P 的第一个蝶的第一 节点标号 k 从 R 开始， 由于本级中旋转因子指数相同的蝶共有 2M L 个， 且这些蝶的相邻间距为 2L ， 故旋转因子指数为 P 的最后一个蝶的第一节点标号 k 为：
实现方框图如下所示：
(4). 变址运算 为了保证运算输出的 X k 按顺序排列，要求 序列
x n 倒序输入，即在运算前要先对输入的序列进行位
序颠倒。一般实践中采用的是将十进制序列先转换为二 进制，将得到的二进制进行倒序，再转换为十进制序列 输出。实现方框图如右图所示： 程序见附录所示。