尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

MATLAB FFT频谱仿真:从DFT原理到参数设置与窗函数选择

MATLAB FFT频谱仿真:从DFT原理到参数设置与窗函数选择 简介面向MATLAB频谱分析学习场景的轻量例程包以快速傅里叶变换FFT的编程实现为核心适合刚接触信号处理、希望快速建立信号频域直观认识的初学者。压缩包整体仅951B虽然总共只有2个文件却覆盖了从时域信号读入、FFT运算到频谱图绘制的完整流程其中一个MATLAB脚本承担核心计算另一个则是自动保存的备份文件不影响主程序运行。已有310人学习该资源其最大价值在于提供一个可直接运行的仿真示例帮助读者通过实际代码理解fft(y)返回复数频谱、abs(fft(y))提取幅度谱、angle(fft(y))获取相位信息等关键操作并在最终绘制的频率-幅度图中观察信号的特征频率、噪声水平与谐波成分。例程同时适合通信、音频处理、图像处理等领域的频谱基础分析通过运行与修改脚本读者可以直观对比周期信号与非周期信号的频谱结构为后续功率谱密度估计、滤波器设计等进阶应用打下基础。1. 为什么 MATLAB 的 FFT 频谱仿真要先搞懂 DFT 点数与频率轴很多人拿到一个 fft.m 例程点击运行看到一条尖峰就以为完事了。实际上我见过不少人在这一步栽跟头横轴画的是“点数”而不是“频率Hz”峰值位置对不上信号频率甚至因为脚本起名 fft.m 导致调用内置 fft 函数直接报错。这个 rar 包里给出的正是一个最简复现流程读入或生成信号、调用 fft、用 plot 画出幅度谱。它适合刚接触频谱分析的学生也适合用 MATLAB 做振动、音频、通信数据预研的工程师。要想把这条流程真正用于仿真先要搞清楚 FFT 的点数 N、采样率 fs 和频率分辨率 fs/N 三者之间的关系否则后面加窗、滤波、分段频谱都无从谈起。2. FFT 的数学底子与 MATLAB fft 函数的关键参数2.1 从 DFT 到 FFT计算量为什么从 N^2 降到 N log N离散傅里叶变换DFT把长度为 N 的时域序列 x(n) 转换为频域序列 X(k)定义式是X(k) sum_{n0}^{N-1} x(n) * exp(-j2pikn/N), k 0, 1, ..., N-1直接用这个式子做完整变换每个频点需要 N 次复数乘加N 个频点就是 N^2 次。当 N65536 时N^2 大约是 4.3e9这个量级在普通电脑上即使预分配矩阵也需要几百毫秒到几秒。FFT 能利用旋转因子 exp(-j2pi/N) 的周期性与对称性把一次大 DFT 递归地拆成多个小 DFT整体运算量降为 O(N log N)。N65536 时约 1e6 次量级两者相差两个数量级。所以频谱仿真里没有理由直接写 DFT 循环matlab 内置的 fft 就是首选。MATLAB 的 fft 实现会根据 N 的因子选择不同策略当 N 是 2 的幂时走基-2 FFTN 含有其他小素因子时走混合基算法完全分解不了时退化为 Bluestein 算法。对实时性要求高的时候尽量让 N 成为 2 的幂这也是后面仿真的一个惯用手段。实际代码里经常出现 fft(x, 2^nextpow2(N))就是为了把 N 凑成 2 的幂。但要注意补零只能让频谱曲线更密不能提高物理分辨率真实分辨率由实际信号长度决定。2.2 fft(y)、abs、angle 与单边频谱的换算matlab 里最基础的调用是 Y fft(x)它返回一个和 x 等长的复数数组下标从 0 到 N-1 对应数字角频率 0 到 2pi(N-1)/N。如果你直接 plot(abs(Y))会看到左右两个对称峰那是正负频率各自贡献的幅度。绝对值是幅度谱angle(Y) 是相位谱单位是弧度。要画出工程上常用的单边幅度谱一般要这么写fs 1000; N 1000; t (0:N-1)/fs; x 0.7*sin(2*pi*50*t) 0.3*sin(2*pi*120*t); Y fft(x); % 双边复数谱 P2 abs(Y/N); % 归一化到实际幅度 P1 P2(1:N/21); % 保留 0 ~ fs/2 P1(2:end-1) 2*P1(2:end-1); % 合并正负频率能量 f fs*(0:(N/2))/N; % 频率轴单位 Hz plot(f, P1, LineWidth, 1.2); xlabel(Frequency (Hz)); ylabel(Amplitude); grid on;代码里最容易错的两行是 P2 abs(Y/N) 和 P1(2:end-1) 2*P1(...)。abs(Y/N) 把 FFT 的输出按点数平均否则 50 Hz 正弦的幅度会显示成 350 而不是 0.7。P1 后半段的前半部分实际上是负频率镜像对于实信号幅度只占一半所以除了 DC 和 fs/2 两个点外都要乘 2。关于 fft 的可选参数fft(x, n) 会补零到 n 点fft(x, n, dim) 可以沿指定维度做批量 FFT适合一次性处理多条振动测点。具体差异可以看下面这张表调用形式返回内容适用场景Y fft(x)与 x 等长的双边复数谱默认按第一维变换Y fft(x, n)n 点 FFT长度不足补零加速到 2 的幂Y fft(x, n, dim)沿 dim 维变换多测点批量处理处理多测点时dim1 表示按列变换dim2 按行变换这个参数在批量仿真中非常有用。如果信号是实数还可以用 rfft 的思路只算一半但 MATLAB 没有内置 rfft实际中直接 fft 再截断更常见。3. 拆解 fft.m 例程数据读取、仿真跑通与 .asv 文件恢复3.1 主脚本流程与关键行功能常见的 fft.m 例程结构非常紧凑核心流程如下clear; close all; clc; Fs 1000; % 采样率单位 Hz L 1500; % 序列长度 t (0:L-1)/Fs; % 时间轴 S 0.8*sin(2*pi*50*t) sin(2*pi*120*t); % 原始信号 Y fft(S); P2 abs(Y/L); P1 P2(1:L/21); P1(2:end-1) 2*P1(2:end-1); f Fs*(0:(L/2))/L; plot(f, P1); xlabel(Frequency (Hz)); ylabel(|P1(f)|);这段代码的每一行都有讲究。Fs1000 意味着奈奎斯特频率是 500 Hz所以上述两个正弦峰分别落在 50 和 120。L1500 不是 2 的幂但 fft 底层会自动按混合基执行速度也不慢。fft(S) 返回 1500 个复数点abs(Y/L) 是归一化幅度谱取 1:(L/21) 去掉负频率镜像最后 f 把 bin 索引换算成物理频率。plot 的结果就是频谱仿真最常看到的单边幅度谱。有个和包名直接相关的坑这个脚本命名为 fft.m和内置函数 fft 重名。在 MATLAB 里当前工作目录中的脚本优先级高于内置函数所以脚本一旦运行内部的 fft(S) 会被解析成调用自身导致报错或奇怪的行为。我拿到这类压缩包的第一件事就是把 fft.m 改名为 myfft_demo.m然后重新运行。这也是为什么网上很多频谱仿真代码会特意避免用内置函数名当脚本名。3.2 采样率、时长与频率分辨率怎么设FFT 频谱仿真里参数设置往往决定结果是可信还是不可用。先把关键参数表放出来参数典型值作用错误后果Fs 采样率1000 Hz决定可分析最高频率 Fs/2小于两倍最高信号频率时产生混叠L 信号长度1500 点决定频率分辨率 df Fs/L太短时邻近谱峰无法分离NFFT 计算点数2048补零插值使曲线平滑被人误认为提高了物理分辨率window 窗函数hann(256)控制泄漏和主瓣宽度窗与信号不匹配时旁瓣污染实际设置时先按信号最高频率 hf 确定 Fs 2*hf再按你最关心的最小频率间隔 df_min 确定 L Fs/df_min。比如要分辨 1 Hz 以内的振动边带Fs1000 时 L 至少 1000 点对应 1 秒采样时长。反之如果采样率给高了L 不变df 会被拉大频率就更粗。这是很多人都卡过的 trade-off。3.3 snrtry.asv 是什么以及如何从自动保存恢复rar 里的 snrtry.asv 是 MATLAB 编辑器自动生成的备份文件命名规则是原文件名加 .asv。当源 .m 文件丢失或 MATLAB 异常崩溃可以把 snrtry.asv 重命名为 snrtry.m 打开继续编辑。恢复后注意检查文件修改时间因为自动保存通常有 5 分钟左右的滞后最后一次改动可能不在里面。也可以用 MATLAB 的 compare 功能把 .asv 和现有 .m 做对比确认差异。如果仿真要换成真实数据常见做法是读 CSVdata readmatrix(vibration.csv); t data(:, 1); x data(:, 2); fs 1 / (t(2) - t(1)); % 根据时间间隔估算采样率readmatrix 会自动跳过表头并处理数值类型比老旧的 csvread 更省心。如果数据中间有缺失值先补一下x fillmissing(x, linear);补值后再做 FFT否则 NaN 会让整段频谱变成空洞。到这里fft.m 例程的替换逻辑就清楚了把 t 和 x 代入主流程L 改为 length(x)就能对真实信号做频谱仿真。4. 频谱仿真实战未知采样率、加窗与谐波/噪声识别4.1 不知道采样率时如何求频率轴实际工程里从采集器导出的 CSV 经常没有采样率字段。这时有两个办法。第一如果时间列有完整时间戳直接用 fs1/mean(diff(t)) 估算第二只知道时间间隔就取中位数避免个别坏点干扰。完全没有任何时间信息时只能在归一化频率上分析横轴范围是 0 到 0.5对应直流到奈奎斯特频率。代码如下Y fft(x); P2 abs(Y/length(x)); P1 P2(1:floor(length(x)/2)1); P1(2:end-1) 2*P1(2:end-1); f_norm (0:length(P1)-1) / length(x); % 归一化频率0~0.5 plot(f_norm, P1);这时的 f_norm 单位是 cycles/sample也就是每个采样间隔内重复的周期数。若后来搞清楚真实采样率是 fs把横轴乘以 fs 就变成 Hz。还有一种反推办法如果信号里有已知特征频率比如 50 Hz 工频先读归一化峰值位置 p则 fs 50 / p。这个方法虽然粗糙但在没有文档的旧数据上很管用。4.2 频谱泄漏与窗函数选择直接对截断的有限长数据做 FFT相当于在矩形窗下分析这会在真实频率附近产生主瓣和旁瓣。主瓣宽度限制分辨率旁瓣会把附近大信号的尾巴铺到其他频点这就是频谱泄漏。减小泄漏的标准手段是加窗。下表列出常用窗的选择倾向窗函数主瓣宽度旁瓣衰减适用场景矩形窗2*df-13 dB瞬态信号、频率精确已知汉宁窗4*df-31 dB随机振动、连续谱海明窗4*df-43 dB语音信号布莱克曼窗6*df-58 dB需要强抑制旁瓣加窗在 MATLAB 里一句话就能做N length(x); w hann(N, periodic); xw x(:) .* w; Y fft(xw); P2 abs(Y / sum(w)); % 除以窗面积而不是 N幅度才接近原值注意 hann(N,periodic) 适合频谱分析因为它满足周期延拓条件。如果没有 Signal Processing Toolbox可以改用 hanning(N)效果几乎一致。除以 sum(w) 是为了补偿窗函数带来的幅度损失汉宁窗的 sum(w) 约等于 N/2所以也有人直接乘 2。加了窗之后主瓣变宽两个很近的谱峰可能粘在一起这是泄漏和分辨率的天然矛盾。4.3 振动频谱图怎么分析谐波、底噪与毛刺拿到频谱图先别急着看最高峰。优先顺序应该是先找基频再看基频整数倍位置是否有等间隔谐波然后评估底噪平台最后检查异常单根谱线。下面的代码用 findpeaks 提取主要谱峰并判断每个峰频率与基频的比例[pks, locs] findpeaks(P1, MinPeakHeight, max(P1)*0.1); f_peaks f(locs); ratio f_peaks / f_peaks(1); harmonic_mask abs(ratio - round(ratio)) 0.02;如果 ratio 接近 1, 2, 3这些峰基本可以定为基频的谐波。而振动问题里常见的 2 倍频、0.5 倍频分别对应不对中、松动等特征。底噪表现为整段频谱抬高的平坦区通常由机械冲击或流噪声引起。如果出现一根非常窄、幅度远超周围的谱线但又不是任何谐波就要怀疑是电干扰或数据采集毛刺先把原始时域波形拉出来看确认是否有单点异常。真正干净的信号频谱包络应当是平滑的一旦出现很多高度参差的尖峰编程上先检查时间轴是否均匀再决定是否需要重采样。5. 分段频谱图验证非平稳信号spectrogram 的窗口与重叠参数5.1 用 spectrogram 观察时变频率与验证加窗效果一次 FFT 只能反映整段数据的平均频谱。齿轮箱启停、语音、变频器输出这类非平稳信号频率随时间漂移必须用分段频谱图短时傅里叶变换。MATLAB 命令是fs 1000; t 0:1/fs:5; x sin(2*pi*(5020*t).*t); % 频率从 50 线性扫到 150 Hz spectrogram(x, hann(256), 192, 512, fs, yaxis); colorbar;这里 window 长度 256 个采样noverlap 192 表示相邻窗口重叠 75%nfft 512 用于细化频点。窗长越短时间分辨率越高但短窗像一个滑动平均频率轴会被拉粗noverlap 越大图像越平滑代价是计算量增大。一般我会从 50% 重叠起步观察图形后上调到 75%很少用 90% 以上因为相邻窗口高度相关视觉上只是浪费时间。分段频谱图同时也是验证整段 FFT 是否正确的好工具。具体做法是取 spectrogram 在某一时刻返回的切片与普通 fft 的结果对比。例如[S, F, T] spectrogram(x, hann(256), 192, 512, fs); middle_slice S(:, ceil(length(T)/2)); plot(F, abs(middle_slice));若切片主峰位置与直接 fft 的主峰一致说明加窗、重叠和归一化的参数没有破坏频域幅度关系。如果最大频率超过奈奎斯特频率spectrogram 会在频带上出现折叠峰此时回到时域检查时间间隔是否均匀。对于离线分析我还会把 spectrogram 结果导出成矩阵再画等高线方便定位谐振点随转速的变化。对照一下窗口长度和频率分辨率就能确定是改窗长还是改重叠率。本文还有配套的精品资源点击获取
返回列表