
简介实验三用FFT对信号作频谱分析实验报告是一份面向数字信号处理课程学习者与相关实验人员的实验报告资源系统讲解如何利用FFT对连续信号和时域离散信号进行谱分析并探讨频谱分辨率与分析误差的来源。报告包含完整的实验原理、步骤与结果分析覆盖有限长序列、周期序列和模拟周期信号三类典型对象分别在不同FFT变换区间N下绘制幅频特性曲线并讨论差异。压缩包为doc格式文档共1个文件大小187KB便于直接阅读与打印。资源目前已有1818人学习。通过这份报告读者可以掌握FFT变换区间选择对频率分辨率和谱分析结果的影响理解离散谱与连续谱逼近的误差原因并能参考附录中的MATLAB实验代码自行复现适合作为实验预习、课程报告撰写或期末复习的参考资料。1. 用 FFT 做频谱分析结果差多少取决于变换区间 N 怎么选这个实验的核心结论其实一句话就能说清FFT 本身没有参数可调真正决定频谱“像不像”的是变换区间 N。同一个 x6(n)cos(8πt)cos(16πt)cos(20πt)采样频率 Fs64HzN 取 16 时第三条谱线 10Hz 几乎分辨不出来N 取 32 时三条线整整齐齐排开。差别不在算法而在 FFT 的频率网格间隔 2π/N 以及截断长度是否满足整数周期条件。实验用三类信号——有限长序列 R4(n)、三角波序列、周期序列和模拟周期信号——分别做 N8/16 和 16/32/64 的谱分析把频谱泄漏、栅栏效应、频谱分辨率这几个最容易被忽略的概念一次性暴露出来。适合正在做数字信号处理课程实验的学生也适合日常用 FFT 判断信号成分、却被“多出来的谱线”困扰的工程师。2. 频谱分辨率 D、栅栏效应与变换区间 N 的数学关系2.1 为什么 FFT 的频率分辨率是 2π/NN 点 FFT 等价于把序列的 DTFT 在单位圆上均匀取 N 个样本X(k) Σ_{n0}^{N-1} x(n) e^(-j2πkn/N)k 0, 1, …, N-1对应的频率点是 ω_k 2πk/N相邻两个 k 的频率间隔固定为 2π/N。这个间隔就是 FFT 能提供的频率网格密度也就是实验报告里说的“FFT 能够实现的频率分辨率是 2π/N”。注意这里说的是网格间距不是谱峰的宽度更不是信号真实频率的估计误差。网格间距决定了两根频率差小于它的谱线不可能被分开。比如 N8 时频率网格间隔是 π/4而 x5 里 cos(πn/8) 的频率是 π/8正好落在两个网格点中间这就是 8 点 FFT 看不出第二根谱线的直接原因。如果换算成模拟频率网格间隔是 ΔfFs/N。Fs64HzN16 时 Δf4HzN32 时 Δf2HzN64 时 Δf1Hz。要把 8Hz 和 10Hz 这两根线分开频率间隔至少要小于 2Hz于是有 2π/N ≤ D 或者 N ≥ Fs/Δf。实验中的 N 一般取 2 的幂但这不是 FFT 的要求只是因为 radix-2 算法在这类长度上效率最高如果分辨率要求算出 N24直接取 24 也能算只是慢一点。2.2 D 的实际换算法从角频率到模拟频率实验报告里的公式 2π/N ≤ D 用的是归一化角频率D 是被分析信号中两个频率分量的最小间隔。实际处理模拟信号时归一化角频率要除以采样周期才能变成物理频率。一个频率为 f0 的余弦采样后数字角频率是 ω2πf0/FsFFT 谱线位置 kf0·N/Fs。于是 2π/N ≤ D 等价于 Fs/N ≤ Δf也就是 N ≥ Fs/Δf。这是一个非常实用的换算x6 里 8Hz 和 10Hz 的间隔 Δf2HzFs64Hz算出来 N 至少是 32实验设计正好卡着这个临界值故意让 N16 出错再让 N32 恢复正常对比效果非常明显。我在复现这个实验时习惯先把 kf0·N/Fs 对每个频率分量算一遍能整除说明谱线落在网格上不能整除说明这一项必然泄漏。这个方法比盯着幅频曲线猜要快得多也适合在嵌入式平台上做 FFT 前预判分辨率是否够用。2.3 栅栏效应谱线刚好错过真实频率DFT 只在离散点 ω_k 上取值相当于透过栅栏看连续谱必然漏掉栅栏缝里的细节。如果信号频率恰好落在两个 bin 之间真实谱峰就永远不会被任何一条离散谱线采到只能看到两侧相对较低的样本这就是栅栏效应。最直接的表现是峰值幅度偏低、频率定位有偏差。对 10Hz 在 N16、Fs64 的情况k2.5实际谱峰在 k2 和 k3 之间两条线的幅度都不是真正的峰值。栅栏效应和频谱泄漏经常同时出现但成因不同栅栏效应是采样网格太稀泄漏是截断长度非整数周期。缓解栅栏效应的常规做法是补零延长 FFT 点数让谱线变密。但要强调补零只是让网格变密并没有增加观察时间所以不能提高分辨率也不能把真实频率很近的两个峰分开。分辨率只取决于采样点数 N观察时长这一点在实验第三部分会有非常直观的验证。2.4 频谱泄漏非整数周期截断的后果周期信号在理论上只有离散谱线。但对 cos(πn/8) 做 8 点 FFT 时信号周期是 16 个样本8 个点只覆盖半个周期截断造成了首尾不连续。数学上看这相当于把无限长的余弦乘上一个 8 点矩形窗频域上就是余弦谱线与矩形窗频谱的卷积能量被摊到整个频带产生一排旁瓣。这就是 N8 时 x5 的幅频图“多出一堆小谱线”的原因。误差来源成因表现缓解方式栅栏效应DFT 只在离散频率点取值峰值幅度偏低、频率定位偏差补零加密谱线或增大 N频谱泄漏截断长度不是信号周期的整数倍谱线展宽、旁瓣出现、幅度下降整数周期截断、加窗、增大 N实验要求比较 N8 和 N16、N16 和 N32本质就是在看这两类误差如何随变换区间 N 变化。注意这两个误差在 N 较小时会叠加泄漏让能量散开栅栏效应让采样点落在降低后的包络上两者共同造成幅度失真。这也是为什么实验结论里会强调“N 要适当选择大一些”——但大多少靠的正是 2π/N ≤ D 这个定量关系。3. 序列 x1、x2、x3 的 FFT 实验矩形窗、三角波与幅度谱不变性3.1 序列构造与 8 点、16 点 DFT 复现原文给出的三个序列中x1(n)R4(n) 是最简单的有限长矩形窗x2(n) 在 0≤n≤3 时为 n1在 4≤n≤7 时为 8-n展开是 [1,2,3,4,4,3,2,1]x3(n) 在 0≤n≤3 时为 4-n在 4≤n≤7 时为 n-3展开是 [4,3,2,1,1,2,3,4]。这里有个值得注意的细节x3 并不是 x2 的时间反转而是 x2 的圆周移位——把 x2 向右循环移位 4 位就得到 x3。这个关系后面会用来核对结果。用 MATLAB 复现时x2 和 x3 用两段向量拼接比手写 8 个元素更不容易出错x1n ones(1,4); % R4(n)n0..3 取 1其余为 0 X1k8 fft(x1n, 8); % fft 第二参数指定变换区间 N8 X1k16 fft(x1n, 16); % N16序列不足部分自动补零 xa 1:4; xb 4:-1:1; % xa[1 2 3 4]xb[4 3 2 1] x2n [xa, xb]; % [1 2 3 4 4 3 2 1] x3n [xb, xa]; % [4 3 2 1 1 2 3 4] X2k8 fft(x2n, 8); X2k16 fft(x2n, 16); X3k8 fft(x3n, 8); X3k16 fft(x3n, 16); N 8; f 2/N*(0:N-1); % 频率轴单位是 ω/π stem(f, abs(X1k8), .);fft(x, N) 的第二个参数是关键序列长度小于 N 时自动补零大于 N 时强制截断。对 R4(n) 做 8 点和 16 点 FFT本质是在它原有的 4 点 DTFT 连续谱上分别按 2π/8 和 2π/16 的间隔重新采样补零不会改变谱的包络形状只会改变采样密度。绘制幅频特性必须用 stem 而不是 plotDFT 输出是离散谱用连续曲线连接会让人误以为中间也有值。频率轴 f2/N*(0:N-1) 的单位是 ω/πk1 对应 ω2π/N换算成 f 就是 2/N所以横坐标出现 0.25、0.5 是正常的它不是 Hz。3.2 读幅频曲线主瓣、零点与旁瓣R4(n) 的 8 点幅度谱很容易手算验证k0 时是 4k1 时约 2.613k2 时是 0k3 时约 1.082k4 时是 0后半段关于 k4 对称。第一次做这个实验的人看到 k2 和 k4 的幅值是 0会以为是程序写错了实际上这是矩形窗 DTFT 的过零点4 点矩形窗的 DTFT 在 ω2π/4π/2 处第一次过零8 点 FFT 的频率网格步长是 2π/8π/4k2 正好落在 ωπ/2 上所以幅值为零。N 换成 16 后网格步长变小过零点位置对应的 k 变成 4同样为零。这说明零点位置由信号本身决定bin 位置由 N 决定两者吻合时就会出现零值谱线。x2(n) 是两端不为 0 的三角窗它的频谱不再有零点主瓣比矩形窗宽旁瓣电平明显低。从卷积角度看三角序列可以看成两个矩形序列的卷积频域是矩形窗频谱的平方旁瓣衰减速度比矩形窗快一倍。实际跑出来的幅频图里x2 的旁瓣贴着主瓣快速衰减而 R4 的旁瓣会出现周期性的零点。x2 的直流分量是 20所以 k0 的谱线幅值 20这也是核对序列构造是否正确的一个标志。3.3 用圆周移位性质核对 x2、x3 的幅度谱x3(n)x2((n4) mod 8)即向右循环移位 4 位。DFT 有一条性质时域圆周移位只改变相位不改变幅度谱。因此无论 N8 还是 N16|X3(k)| 与 |X2(k)| 都应当完全重合。这是一个零成本的代码自检如果两张幅频图画出来不重合先检查序列构造是否把 [1,2,3,4,4,3,2,1] 和 [4,3,2,1,1,2,3,4] 的顺序写错了。isequal(x3n, circshift(x2n, 4)) % 结果为 1 说明构造正确 isequal(x3n, fliplr(x2n)) % 结果为 0说明不是简单左右翻转常见误用是把 x3 当时间反转去验证因为序列形状看起来像反转。实际上 x2 本身左右对称反转后还是自己而 x3 是移出来的两者的相位谱完全不同。做谱分析时如果只关心幅度容易忽略这类差异但后续做相位分析或 FIR 滤波器设计时时移导致的线性相位会直接影响结果。这里花一分钟确认性质比事后排查快得多。4. 周期序列与模拟周期信号泄漏、分辨率与 N 的临界值4.1 x4 与 x5为什么 N8 时 x5 的第二根谱线消失x4(n)cos(πn/4) 的周期是 8 个样本x5(n)cos(πn/4)cos(πn/8) 中第二项的周期是 16 个样本。N8 时x4 恰好截取一个完整周期谱线落在 k1 和 k7幅度各为 4x5 的第一项同样落在 k1但第二项的真实频率 π/8 在两个 bin 之间能量泄漏到全部 8 条谱线所以 N8 的幅频图里只能看到一根主谱线加一圈毛刺。N16 时两个分量都满足整数周期条件x5 的谱线分别落在 k1 和 k2幅度各为 8。这里有个容易被误读的细节DFT 对幅度为 1 的余弦谱峰幅度是 N/2N 越大峰值越高这不是信号变强而是 DFT 没有 1/N 归一化。比较不同 N 下的幅频曲线时不要拿绝对幅度直接对比要看谱线位置和是否存在多余的非零谱线。原报告里“图形与理论分析相符”指的就是这个x4 无论 N 取 8 还是 16都只有两根谱线而 x5 在 N8 时出现拖尾。如果想验证 x5 在 N16 时确实没有泄漏可以直接查非零谱线下标n 0:15; X5k16 fft(cos(pi*n/4) cos(pi*n/8), 16); find(abs(X5k16) 1e-6) % 返回 [2 3 15 16]只有四根谱线索引 2、3 对应 k1、2索引 15、16 对应 k14、15也就是正负频率两侧各两根幅度都是 8。注意这里的 n 必须从 0 开始因为 DFT 公式要求第一个采样点在 n0如果从 1 开始相当于信号循环移位一个样本幅度谱不变但相位全变初学时容易忽略这个问题。4.2 模拟信号 x6 的采样与 16/32/64 点 FFTx6(n)cos(8πt)cos(16πt)cos(20πt) 名义上是连续时间函数实验必须先采样。tnTn/FsFs64Hz三个分量变成 cos(8πn/64)、cos(16πn/64)、cos(20πn/64)模拟频率分别是 4Hz、8Hz、10Hz数字角频率分别是 π/8、π/4、5π/16。实验报告正文写的是 x6(n)附录代码里变量名却叫 x8n复现时统一用 x6n 就行不用管命名不一致。import numpy as np Fs 64 for N in (16, 32, 64): n np.arange(N) # 0..N-1 的采样序号 x (np.cos(8 * np.pi * n / Fs) # 4 Hz 分量 np.cos(16 * np.pi * n / Fs) # 8 Hz 分量 np.cos(20 * np.pi * n / Fs)) # 10 Hz 分量 X np.fft.fft(x) # 默认按输入长度做 FFT mag np.abs(X) print(N , N, 峰值索引 , np.argsort(mag)[-3:][::-1])np.fft.fft 不指定第二个参数时按序列长度计算等价于 MATLAB 的 fft(x)。np.argsort(mag)[-3:][::-1] 取出幅度最大的三个谱线位置结果能直接看到 N 变化带来的差异N16 时最大峰聚集在 k1、k2 附近10Hz 的分量混在 8Hz 的峰里N32 时索引变成 [2,4,5]N64 时变成 [4,8,10]三个频率精确落在对应 bin 上。这组索引就是判断实验是否做对的硬指标。4.3 N 与栅栏效应的对照表把三个频率在不同 N 下对应的 k 值全部列出来问题一目了然NΔfFs/N4Hz 的 k8Hz 的 k10Hz 的 k谱分析结果164 Hz122.5不在网格上8Hz 与 10Hz 无法分辨10Hz 泄漏322 Hz245三条谱线清晰分离641 Hz4810全部整数 bin峰值精确等于 N/2kf0·N/Fs 必须是整数截断才是整数周期。10Hz 在 N16 时算出 2.5这就是分辨不出的根源N32 变成 5一切正常。换个角度看采样时长N/64 是观察时间N16 对应 0.25 秒10Hz 信号在这段时间里有 2.5 个周期N32 对应 0.5 秒10Hz 正好 5 个周期。所以实验比较 N16、32、64 不是看 N 越大越好而是看 N 是否越过整数周期这个临界点。注意N16 的谱线间隔 4Hz 本身已大于 8Hz 与 10Hz 的间距 2Hz信息论意义上就不可能分辨。这不是 FFT 实现的问题是观察时间不够长换任何算法都救不回来。实验报告要求讨论“三种情况的区别”核心论点就两条N32 恰好满足 2π/N ≤ D 的临界条件N64 额外保证了每个分量都是整数周期因此峰值幅度无泄漏、精确为 N/2。实际工程中如果只关心频率位置N 取临界值就够如果还要测幅度建议留出冗余让所有感兴趣频率都落在整数 bin 上。5. 三个核对技巧不重跑实验也能发现 FFT 结果对不对5.1 由峰值索引反推频率并核对 N/2 幅度对幅度为 A 的余弦整数周期截断时 DFT 谱峰幅度是 A·N/2。x4 在 N16 时峰值应为 8x6 在 N64 时三个峰都应为 32。检查时用一条命令X4k16 fft(cos(pi*(0:15)/4), 16); find(abs(X4k16) 1e-6) % 期望 [3, 15]对应 k2 与 k14 abs(X4k16(3)) % 期望 8索引 3 对应 k2因为 MATLAB 数组从 1 开始。如果 find 返回一串连续非零索引说明截断不是整数周期先检查采样长度不要急着改窗函数。5.2 补零与加长 N 的区分fft(x, 256) 与 fft(x, 64) 都合法但含义完全不同。补零到 256 点只是在 64 点网格之间插值加密谱峰位置不移动包络形状不变要真正把 8Hz 和 10Hz 分开必须让实际采样点数满足 N ≥ Fs/Δf。判断方法很简单补零后相邻谱线的包络连接起来跟原来一样增大 N 后包络本身会变化峰变得更尖、旁瓣更低。在 FPGA 上做这个实验时Vivado 的 FFT IP 核会把变换点数固定死在配置里点数定小了只能重综合所以上板之前先用这段代码把 N 的临界值算清楚。5.3 Parseval 能量核对不经归一化的 DFT 满足 Σ|x(n)|² (1/N)Σ|X(k)|²。x6 在 64 点时三个幅度为 1 的余弦各贡献 32 点能量总能量应为 96x6n cos(8*pi*(0:63)/64) cos(16*pi*(0:63)/64) cos(20*pi*(0:63)/64); sum(abs(x6n).^2) % 期望 96 sum(abs(fft(x6n)).^2) / 64 % 期望同样为 96如果数值对不上先检查三个分量是否同时满足整数周期条件再检查采样序号是否从 0 开始。这个核对不依赖任何外部工具箱是离开实验室也能做的最后一道防线。本文还有配套的精品资源点击获取