实验目的
加深对FFT 算法原理的理解（因为FFT 只是DFT
的一种快速算法，所以FFT的运算结果必然满足DFT 的基本性质）。
熟悉FFT 子程序的应用。
用FFT
对连续信号和时域离散信号（如语音和图像信号）进行谱分析的实现方法，能理解可能出现的分析误差，能分析误差的原因，以便在实际应用中能正确使用FFT。
对同一信号，可分别在时域和频域中进行表示或描述，以便从不同域，也就是不同的角度对信号进行特性分析，或进行有效处理，因此要建立从多种角度看待问题的哲学思想和科学思维。
实验原理
DFT的定义
离散傅里叶变换（Discrete
FourierTransform，简称DFT）是一种将有限长度的时域信号（离散信号）转换为有限长度的频域信号的数学方法。DFT是傅里叶变换在离散信号上的应用，它允许我们分析数字信号的频率成分。对于长度为N的序列，其DFT定义为：
其中， 是序列 的DFT， 是频率索引， 是虚数单位，是自然对数的底数。
DFT的性质
线性性质
如果 和 是两个有限长序列，长度分别为 和，且 其中，、
为常数，取，则 的
点 DFT 为 其中 和 分别为 和 的 点 DFT。
循环移位性质
1.序列循环移位
设 为有限长序列，长度为
， ，则 的循环移位定义为
公式表明，将 以 为周期进行周期延拓得到，再将 左移 得到，最后取
的主值序列则得到有限长序列
的循环移位序列 。
2.时域循环移位定理
设 是长度为 的有限长序列， 为 的循环移位，即
则
3.频域循环移位定理
如果
则
用MATLAB计算序列的DFT
MATLAB提供了用快速傅里叶变换算法FFT（算法见第
4章介绍）计算DFT的函数fft,其调用格式如下:
调用参数
为被变换的时域序列向量, 是
DFT变换区间长度,当 大于 的长度时, fft函数自动在 后面补零。函数返回 的 点 DFT变换结果向量 ,这时,
。当 小于 的长度时, fft函数计算 的前面 个元素构成的 长序列的 点 DFT,忽略 后面的元素。
用DFT信号进行谱分析
所谓信号的谱分析，就是计算信号的傅里叶变换。连续信号与系统的傅里叶分析显然不便于直接用计算机进行计算，使其应用受到限制。而DFT是一种时域和频域均离散化的变换，适合数值运算，成为用计算机分析离散信号和系统的有力工具。对连续信号和系统，可以通过时域采样，应用DFT进行近似谱分析。下面分别介绍用DFT对连续信号和离散信号(序列)进行谱分析的基本原理和方法。
用DFT对连续信号进行谱分析
工程实际中，经常遇到连续信号 ，其频谱函数 也是连续函数。为了利用DFT对
进行频谱分析，先对 进行时域采样，得到 ，再对 进行 DFT，得到的X(k)则是
x(n)的傅里叶变换 在频率区间 上的 N点等间隔采样。这里
x(n)和X(k)均为有限长序列。然而，由傅里叶变换理论知道，若信号持续时间有限长，则其频谱无限宽；若信号的频谱有限宽，则其持续时间必然为无限长。所以严格地讲，持续时间有限的带限信号是不存在的。
因此，按采样定理采样时，上述两种情况下的采样序列
均应为无限长，不满足DFT的变换条件。实际上对频谱很宽的信号，为防止时域采样后产生频谱混叠失真，可用预滤波器滤除幅度较小的高频成分，使连续信号的带宽小于折叠频率。对于持续时间很长的信号，采样点数太多，以致无法存储和计算，只好截取有限点进行DFT。
由上述可见，用DFT对连续信号进行频谱分析必然是近似的，其近似程度与信号带宽、采样频率和截取长度有关。实际上从工程角度看，滤除幅度很小的高频成分和截去幅度很小的部分时间信号是允许的。因此，在下面分析中，假设
是经过预滤波和截取处理的有限长带限信号。
设连续信号 持续时间为
，最高频率为 。的傅里叶变换为 ，对 进行时域采样得到， 的傅里叶变换为。由假设条件可知
的长度为
式中，T为采样间隔，
为采样频率。用 X(k)表示 x(n)的 N点DFT，下面推导出X(k)与
的关系，最后由此关系归纳出用X(k)表示 的方法，即用
DFT对连续信号进行谱分析的方法。
由(2.4.3)式知道，
的傅里叶变换 与 的傅里叶变换 满足如下关系：
将
代入上式，得到：
式中
表示模拟信号频谱
的周期延拓函数。由 x(n)的 N点DFT的定义有
将上式代入下式，得到：
上式说明了 与 的关系。为了符合一般的频谱描述习惯，以频率
f为自变量，整理上式。令
则上式变为
由此可得
式中，F表示对模拟信号频谱的采样间隔，所以称之为频率分辨率， 为截断时间长度。
在对连续信号进行谱分析时，主要关心两个问题，这就是谱分析范围和频率分辨率。谱分析范围为，直接受采样频率 的限制。为了不产生频谱混叠失真，通常要求信号的最高频率。频率分辨率用频率采样间隔F描述，F表示谱分析中能够分辨的两个频率分量的最小间隔。显然，F越小，谱分析就越接近，所以
F较小时，我们称频率分辨率较高。下面讨论用DFT对连续信号谱分析的参数选择原则。
在已知信号的最高频率（即谱分析范围）时，为了避免频率混叠现象，要求采样速率
满足下式：
按照上式，谱分辨率 ，如果保持采样点数N不变，要提高频谱分辨率（减小F），就必须降低采样频率，采样频率的降低会引起谱分析范围变窄和频谱混叠失真。如维持
不变，为提高频率分辨率可以增加采样点数 N，因为，只有增加对信号的观察时间，才能增加 和
N可以按照下面两式进行选择：
用DFT对序列进行谱分析
我们知道单位圆上的 Z 变换就是序列的傅里叶变换，即
是
的连续周期函数。如果对序列 进行
N 点 DFT 得到 X(k)，则 X(k) 是在区间 上对 的 N
点等间隔采样，频谱分辨率就是采样间隔。因此序列的傅里叶变换可利用
DFT（即 FFT）来计算。
对周期为 N 的周期序列 ，由(2.3.10)式知道，其频谱函数为
其中
由于 以 N
为周期，因而 也是以 为周期的离散谱，每个周期有 N
条谱线，第 k 条谱线位于 处，代表 的
k次谐波分量。而且，谱线的相对大小与 成正比。由此可见，周期序列的频谱结构可用其离散傅里叶级数系数 表示。由 DFT
的隐含周期性知道，截取 的主值序列 ，并进行 N 点
DFT，得到：
所以可用 X(k)表示
的频谱结构。
如果截取长度 M 等于
的整数个周期，即 ，m为正整数，即
令 ，则
因为
所以
由此可见， 也能表示 的频谱结构，只是在 时，，表示 的 i次谐波谱线，其幅度扩大 m
倍。而其他 k值时，，当然， 与
对应点频率是相等的。所以，只要截取
的整数个周期进行DFT，就可得到它的频谱结构，达到谱分析的目的。
如果
的周期预先不知道，可先通过计算 的自相关函数 估算 的周期，然后按上述方法对周期信号进行谱分析。
用DFT进行谱分析的误差问题
DFT（实际中用FFT计算）可用来对连续信号和数字信号进行谱分析。在实际分析过程中，要对连续信号采样和截断，有些非时限数据序列也要截断，由此可能引起分析误差。下面分别对可能产生误差的三种现象进行讨论。
(1)混叠现象。对连续信号进行谱分析时，首先要对其采样，变成时域离散信号后才能用DFT(FFT)进行谱分析。采样速率
必须满足采样定理，否则会在 （对应模拟频率）附近发生频谱混叠现象。这时用DFT分析的结果必然在
附近产生较大误差。因此，理论上必须满足 （为连续信号的最高频率）。对 确定的情况，一般在采样前进行预滤波，滤除高于折叠频率
的频率成分，以免发生频谱混叠现象。
(2)栅栏效应。我们知道，N点 DFT是在频率区间 上对时域离散信号的频谱进行
N点等间隔采样，而采样点之间的频谱是看不到的。这就好像从
N个栅栏缝隙中观看信号的频谱情况，仅得到
N个缝隙中看到的频谱函数值。因此称这种现象为栅栏效应。由于栅栏效应，有可能漏掉（挡住）大的频谱分量。为了把原来被"栅栏"挡住的频谱分量检测出来，就必须提高频率分辨率。对有限长序列，可以在原序列尾部补零；对无限长序列，可以增大截取长度及DFT变换区间长度，从而使频域采样间隔变小，增加频域采样点数和采样点位置，使原来漏掉的某些频谱分量被检测出来。对连续信号的谱分析，只要采样速率足够高，且采样点数满足频率分辨率要求（见(3.4.14)式和(3.4.15)式），就可以认为DFT
后所得离散谱的包络近似代表原信号的频谱。
(3)截断效应。实际中遇到的序列 可能是无限长的，用
DFT对其进行谱分析时，必须将其截短，形成有限长序列y(n) = x(n) w(n)w(n)$ 称为窗函数，长度为N。，称为矩形窗函数。根据傅里叶变换的频域卷积定理，有
其中
幅度谱
曲线如图1所示（ 以 为周期，只画低频部分）。图中，的部分称为主瓣，其余部分称为旁瓣。
矩形窗的频谱图
例如，，其频谱为
x(n)的频谱 如图2(a)所示。将x(n)截断后， 的幅频曲线如图2(b)所示。
由上述可见，截断后序列的频谱
与原序列频谱必然有差别，这种差别对谱分析的影响主要表现在如下两个方面：
(1)泄露。由图 2(b)可知，原来序列
x(n)的频谱是离散谱线，经截断的离散谱线向附近展宽，通常称这种展宽为泄露。显然，泄露使频谱变模糊，使谱分辨率降低。从图2
可以看出，频谱泄露程度与窗函数幅度谱的主瓣宽度直接相关，在第7章将证明，在所有的窗函数中，矩形窗的主瓣是最窄的，但其旁瓣的幅度也最大。所以，在窗函数长度N
相同时，用矩形窗截取，产生的泄露最小。
(2)谱间干扰。在主谱线两边形成很多旁瓣，引起不同频率分量间的干扰（简称谱间干扰），特别是强信号谱的旁瓣可能湮没弱信号的主谱线，或者把强信号谱的旁瓣误认为是另一频率的信号的谱线，从而造成假信号，这样就会使谱分析产生较大偏差。由于矩形窗的旁瓣最大，所以，用矩形窗截取时，产生的谱间干扰最大。
由于上述两种影响是由对信号截断引起的，因此称之为截断效应。由图1可以看出，增加N可使 的主瓣变窄，减小泄露，提高频率分辨率，但旁瓣的相对幅度并不减小。为了减小谱间干扰，应用其它形状的窗函数代替矩形窗（窗函数将在FIR数字滤波器设计中介绍）。但在N一定时，旁瓣幅度越小的窗函数，其主瓣就越宽。所以，在DFT变换区间（即截取长度）N一定时，只能以降低谱分析分辨率为代价，换取谱间干扰的减小。通过进一步学习数字信号处理的功率谱估计等现代谱估计内容可知，减小截断效应的最好方法是用近代谱估计的方法。但谱估计只适用于不需要相位信息的谱分析场合。
最后要说明的是，栅栏效应与频率分辨率是不同的两个概念。如果截取长度为N的一段数据序列，则可以在其后面补N个零，再进行2N点DFT，使栅栏宽度减半，从而减轻了栅栏效应。但是这种截短后补零的方法不能提高频率分辨率。因为截短已经使频谱变模糊，补零后仅使采样间隔变小，但得到的频谱采样的包络仍是已经变模糊的频谱，所以频率分辨率没有提高。因此，要提高频率分辨率，就必须对原始信号截取的长度加长（对模拟信号，就是增加采样时间
的长度）。
实验步骤
复习DFT 的定义、性质和用DFT 作谱分析的有关内容。
复习FFT 算法原理与编程思想，并对照DIT-FFT
运算流图和程序框图，读懂
本实验提供的FFT 子程序。
- 编制信号产生子程序，产生以下典型信号供谱分析用：
应当注意，如果给出的是连续信号 ，则首先要根据其最高频率确定采
样速率
以及由频率分辨率选择采样点数N，然后对其进行软件采样（即计算
)，产生对应序列 。对信号 ，频率分辨率的
选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周
期的整数倍进行谱分析，否则有可能产生较大的分析误差。请实验者根据DFT
的隐含周期性思考这个问题。
实验原理
【实例4-1】对x3 做8 点FFT，绘制出离散幅度谱。
1
2
3
4
5
6
| N=8;
x3=[4:-1:1 1:4];
xk=fft(x3,N);
figure;
subplot(211); stem(0:length(x3)-1,x3,'.'); title('x3 的波形');
subplot(212); stem(0:N-1,abs(xk),'.'); title('x3 的8 点离散幅度谱');
|
1．编写 Matlab
M 文件对信号x1(n)做 8 点和 16 点的 FFT
1
2
3
4
5
6
7
8
9
10
11
| close;clear;
N1=8;
N2=16;
x1=[1 1 1 1 0 0 0 0];
xk1=fft(x1,N1);
xk2=fft(x1,N2);
figure;
subplot(221); stem(0:length(x1)-1,x1,'.'); title('x1 的波形');
subplot(222); stem(0:N1-1,abs(xk1),'.'); title('x1 的8 点离散幅度谱');
subplot(223); stem(0:length(x1)-1,x1,'.'); title('x1 的波形');
subplot(224); stem(0:N2-1,abs(xk2),'.'); title('x1 的16 点离散幅度谱');
|
2、编写Matlab
M 文件对信号 x2(n)做 8 点和 16 点的 FFT
1
2
3
4
5
6
7
8
9
10
11
| close;clear;
N1=8;
N2=16;
x2=[1:4 4:-1:1];
xk1=fft(x2,N1);
xk2=fft(x2,N2);
figure;
subplot(221); stem(0:length(x2)-1,x2,'.'); title('x2 的波形');
subplot(222); stem(0:N1-1,abs(xk1),'.'); title('x2 的8 点离散幅度谱');
subplot(223); stem(0:length(x2)-1,x2,'.'); title('x2 的波形');
subplot(224); stem(0:N2-1,abs(xk2),'.'); title('x2 的16 点离散幅度谱');
|
3、编写Matlab
M 文件对信号 x4(n) 做 8 点和 16 点的 FFT
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
| close;
clear;
N1=8;
N2=16;
n1=0:1:8;
n2=0:1:16;
x4_1=cos(pi*n1/4);
x4_2=cos(pi*n2/4);
xk1=fft(x4_1,N1);
xk2=fft(x4_2,N2);
figure;
subplot(221); stem(0:length(x4_1)-1,x4_1,'.'); title('x4 的波形');
subplot(222); stem(0:N1-1,abs(xk1),'.'); title('x4 的8 点离散幅度谱');
subplot(223); stem(0:length(x4_2)-1,x4_2,'.'); title('x4 的波形');
subplot(224); stem(0:N2-1,abs(xk2),'.'); title('x4 的16 点离散幅度谱');
|
4、编写Matlab
M 文件对信号 x6(t) 以 fs=64（Hz）采样后做 N=16、32、64 点的FFT
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| clear;close;
fs=64;
N=16;
N1=32;
N2=64;
t = (0:N-1)/fs;
t1 = (0:N1-1)/fs;
t2 = (0:N2-1)/fs;
x6 = cos(8*pi*t) + cos(16*pi*t) + cos(20*pi*t);
x61 = cos(8*pi*t1) + cos(16*pi*t1) + cos(20*pi*t1);
x62 = cos(8*pi*t2) + cos(16*pi*t2) + cos(20*pi*t2);
xk=fft(x6,N);
xk1=fft(x61,N1);
xk2=fft(x62,N2);
figure;
subplot(411); stem(0:length(x6)-1,x6,'.'); title('x6 的波形');
subplot(412); stem(0:N-1,abs(xk),'.'); title('x6 的16 点离散幅度谱');
subplot(413); stem(0:N1-1,abs(xk1),'.'); title('x6 的32 点离散幅度谱');
subplot(414); stem(0:N2-1,abs(xk2),'.'); title('x6 的64 点离散幅度谱');
|
5、编写
Matlab M 文件读取 motherland.wav 音频 数据，分析第 8 000 至 8199
共200个采样点的频谱（提示 这里的频谱指的 是 信号的 傅里叶变换）。
实现方法为
对这 200 个点数据做 N 512 的 DFT （采用 F FT 实现）。要求：画出其在[
0,2 π）
的连续幅度谱和相位谱图。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
| close;clear;
[audio,fs] = audioread(...
".\2024-2025（1）《信号处理实验》资料（请拷贝到u盘中，每次试验带到实验室）\实验中用到的一些函数和音频图像文件\motherland.wav"...
);
audioSegment = audio(8000:8199, :);
N=512;
t = (0:N-1)/fs;
xk=fft(audioSegment,N);
f = linspace(0, 2, N);
figure;
subplot(2, 1, 1);
plot(f, abs(xk));
title('幅度谱');
xlabel('\omega/\pi ');
ylabel('幅度');
subplot(2, 1, 2);
plot(f, angle(xk));
title('相位谱');
xlabel('\omega/\pi');
ylabel('相位 (rad)');
|
1
2
3
4
5
6
7
| close;clear;
A = imread('.\2024-2025（1）《信号处理实验》资料（请拷贝到u盘中，每次试验带到实验室）\实验中用到的一些函数和音频图像文件\lena.bmp');
figure(1);
imshow(A);
title('原图');
|
1
2
3
4
5
6
7
8
9
10
11
12
13
|
fftI = fft2(double(A),512,512);
A1 = abs(fftI);
B1 = (A1-min(min(A1)))/(max(max(A1))-min(min(A1)))*255;
figure(2);
imshow(B1);
title('二维幅度谱图');
|
1
2
3
4
5
|
B2 = fftshift(B1);
figure(3);
imshow(B2);
title('移到中心位置的二维频谱图');
|
思考题
1.
在 N=8 和 N=16 两种情况下， x2(n)和 x3(n)
的幅频特性会相同吗？为什么？
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| close;clear;
N1=8;
N2=16;
x2=[1:4 4:-1:1];
xk1=fft(x2,N1);
xk2=fft(x2,N2);
x3=[4:-1:1 1:4];
xk3=fft(x3,N1);
xk4=fft(x3,N2);
figure;
subplot(321); stem(0:length(x2)-1,x2,'.'); title('x2 的波形');
subplot(322); stem(0:length(x3)-1,x3,'.'); title('x2 的波形');
subplot(323); stem(0:N1-1,abs(xk1),'.'); title('x2 的8 点离散幅度谱');
subplot(324); stem(0:N2-1,abs(xk2),'.'); title('x2 的16 点离散幅度谱');
subplot(325); stem(0:N1-1,abs(xk3),'.'); title('x3 的8 点离散幅度谱');
subplot(326); stem(0:N2-1,abs(xk4),'.'); title('x3 的16 点离散幅度谱');
|
由 MATLAB 绘图可以发现，N=8时，x2(n) 和 x3(n)
的幅频特性是相同的，因为x3(n)=x2((n+4))R8(n)，为循环移位关系，所以 x3(n)
与 x2(n) 的DFT的幅频特性相同，如图所示
但是，当 N=16 时，x3(n) 与 x2(n)
就不满足循环移位关系了，所以幅频特性不同
2.
如果周期信号的周期预先不知道，如何用FFT 进行分析？
基于MATLAB的数字信号处理(3)
用FFT对信号作频谱分析-云社区-华为云 (huaweicloud.com)
周期信号的周期预先不知道时，可先截取 M
点进行DFT，再将截取长度扩大一倍截取，比较结果，如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱，如果不满足误差要求就继续将截取长度加倍，重复比较，直到结果满足要求。
3. 已知序列x=[1,1,2,2,3,3,2,2,1,1]。
（1）对x 进行2 选1
的抽取，得到序列x1=[1,2,3,2,1]；
1
2
3
| clear;
x=[1,1,2,2,3,3,2,2,1,1];
x1=x(1:2:end)
|
（2）
对x 在两个序列值之间插入一个0 值点进行内插，
得到序列x2=[1,0,1,0,2,0,2,0,3,0,3,0,2,0,2,0,1,0,1]。
1
2
3
4
| x2 = zeros(1, length(x) * 2 - 1);
x2(1:2:end) = x;
x2(2:2:end) = 0;
x2
|
1
2
3
| x2 = 1x19
1 0 1 0 2 0 2 0 3 0 3 0 2 0 2 0 1 0 1
|
（3）使用函数fft
分别画出x、x1 和x2 在[0,2π)的连续幅度谱图（提示是序列FT）；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
N = 16*length(x);
X = fft(x, N);
X1 = fft(x1, N);
X2 = fft(x2, N);
f = linspace(0, 2, N);
figure;
subplot(3, 1, 1);
plot(f, abs(X));
title('序列 x 的幅度谱');
xlabel('\omega/\pi');
ylabel('幅度');
subplot(3, 1, 2);
plot(f, abs(X1));
title('序列 x1 的幅度谱');
xlabel('\omega/\pi');
ylabel('幅度');
subplot(3, 1, 3);
plot(f, abs(X2));
title('序列 x2 的幅度谱');
xlabel('\omega/\pi');
ylabel('幅度');
|
（4）分别写出x1、x2
与x 频谱关系的数学表达式。根据（3）显示出的幅度频谱图，解释x1、x2 与x
幅度频谱变化的原因。
序列 与 的频谱关系:
抽取操作会导致频谱的周期性重复。抽取因子为 2，所以 的频谱是 的频谱以 2
为周期的周期性重复。
这意味着 的频谱是 的频谱在频率轴上每隔 弧度压缩一次。
- 序列 与 的频谱关系:
内插操作会导致频谱的周期性扩展。内插因子为 2，所以 的频谱是 的频谱以 2 为周期的周期性扩展。
这意味着 的频谱是 的频谱在频率轴上每隔 弧度扩展一次。
结论：
抽取操作导致频谱的周期性重复，可能导致频谱的混叠。内插操作导致频谱的周期性扩展，可能导致频谱的平滑。