1. 项目概述从时域到频域的工程视角信号处理工程师拿到一段时域波形第一反应往往不是看它的起伏而是问“它的频谱长什么样” 这就像医生看心电图不仅要看心跳的节奏时域更要分析其频率成分是否健康。Matlab 中的 FFT快速傅里叶变换就是我们完成这个“频谱体检”最核心、最常用的“听诊器”。这个项目标题“Matlab的信号频谱分析——FFT变换”看似简单但它背后串联的是从理论到实践从理想模型到工程落地的完整链条。我处理过从音频分析、振动诊断到通信解调的无数案例深知仅仅调用fft()函数是远远不够的如何正确设置参数、理解结果物理意义、避开频谱泄露和栅栏效应这些坑才是从“会用”到“精通”的关键。这篇文章我就以一个老工程师的视角拆解用 Matlab 做 FFT 频谱分析时那些必须掌握的细节和容易踩的坑让你拿到的频谱图不仅好看更要准确、可信。2. 核心原理与工程意义为什么是 FFT在深入代码之前我们必须搞清楚 FFT 到底是什么以及它为什么在工程领域不可替代。傅里叶变换的本质是告诉我们一个复杂的信号可以分解成许多不同频率、不同幅度、不同相位的正弦波的叠加。而 FFT是离散傅里叶变换DFT的一种快速算法它将计算复杂度从 O(N²) 降低到 O(N log N)这才使得在普通计算机上实时分析海量数据成为可能。2.1 FFT 与 DFT 的关系从“算盘”到“计算器”你可以把 DFT 理解成一种定义对于 N 个离散的时域数据点通过一套固定的复数乘加公式计算出 N 个离散的频域点。这套公式是精确的但计算量巨大就像用算盘做大规模乘法。FFT 则是一套巧妙的“心算”技巧它利用正弦函数的对称性和周期性将一个大 DFT 分解成多个小 DFT 的组合从而极大地减少了重复计算。在 Matlab 中我们用的fft(x)函数就是实现了这种高效算法的黑箱。但作为工程师我们不能只满足于黑箱至少要明白它的输入输出特性它假设你给它的 N 点数据是某个无限长周期信号的一个周期然后帮你算出这个假设周期信号的频谱。2.2 频谱的物理意义幅度谱、相位谱与功率谱FFT 直接输出的结果是复数包含实部和虚部。这对数学是完备的但对工程师不够直观。我们通常从中提取三种更有工程意义的谱幅度谱每个频率分量的强度大小。abs(fft(x))就能得到。这是最常用的用于查找信号中的主频、谐波。相位谱每个频率分量的初始相位。angle(fft(x))得到。在图像处理、通信系统同步和某些故障诊断中至关重要。功率谱幅度谱的平方或归一化的幅度平方。它反映信号功率在频域的分布在能量分析中常用单位是 dB。一个关键概念是单边谱和双边谱。由于实数信号的频谱具有共轭对称性FFT 输出的前半部分正频率和后半部分负频率是对称的。为了更直观我们通常只显示正频率部分并将除直流分量0Hz和奈奎斯特频率点外的幅度乘以2以补偿被折叠掉的负频率能量这就是单边幅度谱。这是工程展示的标准做法。3. 完整实操流程从数据导入到频谱图解读理论说再多不如动手做一遍。下面我以一个包含 50Hz 和 120Hz 正弦波的混合信号为例展示从零开始的完整分析流程并穿插解释每一个参数的选择依据。3.1 环境准备与信号合成首先我们合成一个干净的测试信号这有助于验证我们流程的正确性。% 1. 基本参数设置 Fs 1000; % 采样频率单位 Hz。根据香农采样定理要大于信号最高频率的2倍。 T 1/Fs; % 采样间隔 L 1500; % 信号长度点数。选择1500而非1024是为了演示非2的幂次时FFT依然高效。 t (0:L-1)*T; % 时间向量 % 2. 合成信号一个包含直流、基波和高频噪声的典型工程信号 f1 50; % 基频 50Hz f2 120; % 二次谐波 120Hz S 0.7*sin(2*pi*f1*t) sin(2*pi*f2*t); % 两个正弦波叠加 % 模拟实际情况加入直流偏移和高斯白噪声 X S 2.5*cos(2*pi*300*t); % 加入一个强干扰的高频分量300Hz X X 0.5*randn(size(t)); % 加入高斯白噪声模拟传感器噪声这里有几个经验点采样频率 Fs必须事先知道或设定。它决定了你能分析的最高频率Fs/2即奈奎斯特频率。如果信号本身有高于 Fs/2 的频率会发生混叠频谱将完全失真且不可逆。在数据采集阶段就必须用抗混叠滤波器。信号长度 L不一定非要选 2 的幂如1024。早期FFT算法对2的幂有优化但现代Matlab的fft函数对任意长度的输入都采用了高效算法。选择 L 更多考虑的是频率分辨率。3.2 执行 FFT 与频谱计算% 3. 执行FFT Y fft(X); % Y是复数长度也为L % 4. 计算双边谱和单边谱 P2 abs(Y/L); % 计算双边幅度谱并除以L进行归一化解释见下文 P1 P2(1:L/21); % 取前半部分包含直流和奈奎斯特频率点 P1(2:end-1) 2*P1(2:end-1); % 除直流和奈奎斯特点外幅度乘2得到单边谱 % 5. 构建频率轴 f Fs*(0:(L/2))/L; % 单边谱对应的频率向量从0Hz到Fs/2关键解释abs(Y/L)为什么除以 L这是为了幅度归一化。使得时域中幅度为 A 的正弦波在频域对应谱线的幅度也是 A。如果不除谱线幅度将是 A*L/2失去了物理意义。这是新手最常忽略的一步。P1(2:end-1) 2*P1(2:end-1)这就是生成单边谱的核心操作。因为能量平均分布在正负频率只显示正频率时需要将幅度加倍功率则需要乘4。直流分量0Hz和奈奎斯特频率Fs/2是特例它们没有对应的负频率成分所以不加倍。3.3 结果可视化与解读% 6. 绘制时域图和频谱图 figure(Position, [100, 100, 1200, 500]) % 子图1时域信号 subplot(1,2,1) plot(t(1:200), X(1:200)) % 只画前200个点便于观察波形 title(时域信号 (前200点)) xlabel(时间 (s)) ylabel(幅度) grid on % 子图2单边幅度频谱 subplot(1,2,2) plot(f, P1) title(单边幅度频谱) xlabel(频率 (Hz)) ylabel(|幅度|) xlim([0, Fs/2]) % 通常只显示到奈奎斯特频率 grid on % 标记出主要频率峰值 [~, locs] findpeaks(P1, MinPeakHeight, max(P1)*0.1); % 简单找峰 hold on plot(f(locs), P1(locs), ro, MarkerFaceColor, r) text(f(locs)5, P1(locs), cellstr(num2str(round(f(locs)))), VerticalAlignment,bottom) hold off运行这段代码你会得到清晰的时域/频域对比图。在频谱图上你应该能清晰地看到 50Hz 和 120Hz 处尖锐的谱线幅度大约为 0.7 和 1。300Hz 处也有一个明显的峰。而噪声则表现为整个频带底部的“毛刺”基底。注意findpeaks是信号处理工具箱里的一个非常实用的函数用于自动寻找局部极大值。在实际分析中我们经常需要用它来提取主频、谐波频率及其幅度。4. 高级议题与工程陷阱规避如果只是分析理想合成信号那太简单了。工程中的信号是“脏”的数据是有限的下面这些才是真正考验功力的地方。4.1 频谱泄露与加窗函数我们的 FFT 隐含了一个假设截取的那段数据长度L是周期信号的整数个周期。如果不是就会发生频谱泄露——能量从主频点“泄露”到旁边的频点导致谱线变宽、幅值不准旁边还会出现虚假的旁瓣。解决方案就是加窗。窗函数在时域上对信号两端进行平滑衰减减少截断带来的突变。Matlab 提供了hamming,hann,blackman,flattopwin等函数。% 加窗处理示例 win hann(L); % 生成汉宁窗转置成行向量 X_windowed X .* win; % 时域点乘窗函数 % 对加窗后的信号做FFT并修正窗函数带来的幅度损失 Y_win fft(X_windowed); P2_win abs(Y_win / (sum(win)/L)); % 关键归一化因子变为窗函数的平均高度 P1_win P2_win(1:L/21); P1_win(2:end-1) 2*P1_win(2:end-1); figure; plot(f, P1, b, LineWidth, 1.5); hold on; plot(f, P1_win, r--, LineWidth, 1.2); legend(无窗, 汉宁窗); title(加窗对频谱的影响); xlabel(频率 (Hz)); ylabel(|幅度|); grid on;关键点加窗后归一化因子不再是信号长度 L而是窗函数的和sum(win)或者更常用的是窗函数的平均高度sum(win)/L。这是因为窗函数削弱了信号两端的能量直接用 L 除会低估幅度。汉宁窗能有效抑制旁瓣但主瓣会略微变宽频率分辨率下降这是抑制泄露必须付出的代价。选择窗函数本质是在频谱分辨率主瓣宽度和频谱泄露旁瓣高度之间做权衡。4.2 频率分辨率与栅栏效应频率分辨率Δf Fs / L。它表示频谱图上相邻两个频点之间的间隔。如果你的信号中有两个频率分量 f1 和 f2只有当 |f1 - f2| Δf 时才能在频谱图上被区分开。增加 L采集更长时间的数据或降低 Fs在满足采样定理的前提下可以提高分辨率。栅栏效应是离散采样的固有特性FFT 只计算频率为 k*Δf (k0,1,2,...) 这些离散点上的频谱。如果信号的真实频率正好落在两个“栅栏”之间那么它的能量就会分散到相邻的频点上即使没有泄露幅值测量也会不准。解决方案增加数据长度 L这是最根本的方法直接提高分辨率。使用高分辨率谱估计方法如 Chirp-Z 变换可以对特定频段进行“细看”。补零在信号末尾补零后再做 FFTfft(x, N)其中 N length(x)。这不能提高真实的频率分辨率但可以通过对频谱进行插值让曲线更光滑有助于更精确地通过曲线拟合来定位峰值频率是一种“视觉增强”手段。% 演示栅栏效应和补零 L_short 128; % 短数据分辨率低 t_short (0:L_short-1)*T; X_short 0.7*sin(2*pi*50*t_short) sin(2*pi*55*t_short); % 两个很近的频率 % 不补零 Y1 fft(X_short); P1_short 2*abs(Y1(1:L_short/21)/L_short); f_short Fs*(0:(L_short/2))/L_short; % 补零到1024点 N_fft 1024; Y2 fft(X_short, N_fft); P1_zpad 2*abs(Y2(1:N_fft/21)/L_short); % 归一化仍用原数据长度 f_zpad Fs*(0:(N_fft/2))/N_fft; figure; subplot(2,1,1); stem(f_short, P1_short, b, LineWidth, 1.5); % 用 stem 更显离散性 title([短数据 (L, num2str(L_short), )分辨率低栅栏效应明显]); xlabel(频率 (Hz)); ylabel(|幅度|); grid on; xlim([40, 70]); subplot(2,1,2); plot(f_zpad, P1_zpad, r-, LineWidth, 1.2); title([补零到, num2str(N_fft), 点频谱插值更光滑但分辨率未变]); xlabel(频率 (Hz)); ylabel(|幅度|); grid on; xlim([40, 70]);从图中可以清晰看到短数据时两个频率峰混叠在一起补零后谱线变密能更好地描绘出包的形状但两个峰依然无法分开证明物理分辨率未变。4.3 平均与平滑从瞬时谱到统计谱对于平稳随机信号如噪声或为了抑制分析中的随机波动我们通常不直接分析一段数据的频谱而是计算平均功率谱密度。Welch 方法这是工程上的标准方法。它将长数据分段可重叠对每一段加窗并计算周期图单个段的功率谱最后对所有段的周期图求平均。Matlab 中的pwelch函数实现了它。[pxx, f_welch] pwelch(X, hann(256), 128, 1024, Fs); % 窗长256重叠128FFT点数1024 figure; plot(f_welch, 10*log10(pxx)); % 以dB为单位绘制 title(Welch方法估计的平均功率谱密度 (PSD)); xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); grid on;Welch 方法通过平均显著平滑了频谱降低了方差更能反映信号的统计特性特别适合分析噪声和随机振动信号。滑动平均平滑对于已经计算出的幅度谱也可以进行频域平滑比如使用移动平均滤波器。但这是一种后处理会损失频率分辨率。P1_smooth movmean(P1, 5); % 5点移动平均5. 实战案例电机振动信号分析让我们用一个更接近实战的场景来串联以上所有知识。假设我们采集了一段电机轴承的振动加速度信号采样频率 Fs 10 kHz数据长度 N 20000 点即2秒数据。我们怀疑轴承存在故障其故障特征频率约为 120 Hz。% 模拟电机振动信号 (包含故障频率和宽带噪声) Fs_motor 10000; t_motor (0:19999)/Fs_motor; f_fault 120; % 假设的故障特征频率 vibration 1.5 * sin(2*pi*f_fault*t_motor) ... % 故障特征 0.3 * sin(2*pi*2*f_fault*t_motor) ... % 二次谐波 0.1 * sin(2*pi*3*f_fault*t_motor) ... % 三次谐波 randn(size(t_motor)); % 强烈的随机振动噪声 % 1. 直接FFT分析效果可能不佳 L_motor length(vibration); Y_raw fft(vibration); P1_raw 2*abs(Y_raw(1:L_motor/21)/L_motor); f_motor Fs_motor*(0:(L_motor/2))/L_motor; % 2. 使用Welch方法估计PSD突出周期性成分 [pxx, f_welch] pwelch(vibration, hann(2048), 1024, 4096, Fs_motor); % 绘图对比 figure(Position, [50, 50, 1400, 600]); subplot(1,2,1); plot(f_motor, 10*log10(P1_raw.^2)); % 将幅度谱转为功率谱粗略对比 title(直接FFT得到的频谱 (dB)); xlabel(频率 (Hz)); ylabel(功率 (dB)); grid on; xlim([0, 500]); % 标记故障频率及其谐波 hold on; plot([f_fault, 2*f_fault, 3*f_fault], [-20, -30, -35], rv, MarkerFaceColor, r); hold off; subplot(1,2,2); plot(f_welch, 10*log10(pxx)); title(Welch方法估计的功率谱密度 (PSD)); xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); grid on; xlim([0, 500]); hold on; plot([f_fault, 2*f_fault, 3*f_fault], [max(10*log10(pxx))-10, -45, -50], rv, MarkerFaceColor, r); hold off; legend(PSD, 故障频率点);在这个案例中直接 FFT 的频谱被强大的背景噪声淹没故障频率的峰值并不明显。而经过 Welch 平均后的 PSD 图噪声基底变得平坦120Hz、240Hz、360Hz 处的故障特征频率及其谐波清晰地凸现出来这对于故障诊断具有决定性意义。6. 常见问题与调试技巧实录在实际操作中你一定会遇到各种奇怪的现象。下面是我总结的一些典型问题及排查思路。问题现象可能原因排查与解决方法频谱幅值不对远大于或小于预期未进行幅度归一化。FFT结果未除以信号长度L。计算幅度谱时务必使用abs(fft(x)/N)其中N是参与FFT运算的数据点数。频谱图在中间频率出现对称的“镜像”峰错误地绘制了双边谱。对于实数信号这是正常现象但通常我们只看单边谱。确保你只取了FFT结果的前半部分N/21点并对幅度进行了乘2处理直流和奈奎斯特点除外。单一频率的正弦波频谱却是一个很宽的“包”频谱泄露。信号截取长度不是信号周期的整数倍。加窗处理。使用hann或hamming窗。同时尽量采集更长时间的数据使截取长度接近周期的整数倍。已知信号频率为f0但频谱峰值在f0旁边栅栏效应。f0 不在频率分辨率的整数倍上。1.增加数据长度L以提高分辨率。2. 在数据后补零并进行FFT插值然后通过抛物线插值或寻找最大值点来更精确地估计真实频率。高频部分出现不应该有的低频成分频谱混叠。信号中包含高于奈奎斯特频率Fs/2的成分。这是硬件问题无法通过软件完全修复。必须在ADC采样前使用模拟抗混叠滤波器将高于Fs/2的频率成分滤除。检查你的采样率是否足够。功率谱密度PSD的计算结果单位不对对PSD的定义和归一化理解有误。pwelch默认返回的是单边PSD单位是x^2/Hz。理解pwelch的输出。如果要转换为 dB使用10*log10(pxx)。确保你的窗函数参数设置合理pwelch内部已经考虑了窗函数的能量归一化。findpeaks找不到正确的峰值或找到太多杂峰阈值设置不当。噪声基底过高。使用findpeaks(P1, MinPeakHeight, threshold)设置绝对阈值或MinPeakProminence设置最小峰凸起度。通常可以先估算噪声水平将阈值设为噪声水平的3-5倍。处理大量数据时FFT速度慢数据长度过长或循环中多次调用FFT。1. 确保数据长度是许多小素数的乘积Matlab的FFT对此优化最好。2. 考虑使用分段处理或实时频谱分析技术。3. 对于固定长度的FFT可以预先计算旋转因子。几个私房调试技巧从简单信号开始验证任何新的频谱分析流程先用一个幅度、频率已知的纯净正弦波测试。确保频谱图上在正确位置出现一根干净、幅度正确的谱线。这是检验你流程是否正确的“试金石”。关注频率轴很多错误源于频率轴构建错误。反复检查f Fs*(0:(N/2))/N这个公式。确保你的频率向量长度与单边谱数据长度完全一致。理解fft函数的第二个参数 NY fft(X, N)指定了进行 N 点 FFT。如果 N length(X)会自动补零如果 N length(X)会截断 X。这个特性在需要固定长度FFT或进行补零时非常有用。使用fftshift可视化fftshift可以将FFT输出的零频分量移动到频谱中心这对于观察以零频对称的信号如基带信号很方便。但注意这之后频率轴也需要相应调整。保存中间结果在编写复杂的分析脚本时将关键的中间变量如加窗后的信号、原始的FFT复数结果、频率向量保存下来或单独绘图检查能帮你快速定位问题出在哪一步。最后记住频谱分析是一门“艺术”需要在分辨率、精度、速度和平滑度之间根据具体应用做取舍。没有一种设置能通吃所有场景。多动手多对比用已知信号去验证你的流程你的“频谱直觉”就会慢慢建立起来。当你拿到一段陌生的信号能迅速在脑海中勾勒出它大致的频谱模样并知道用什么工具和方法去验证时你就真正掌握了这门技能。
1. 项目概述从时域到频域的工程视角信号处理工程师拿到一段时域波形第一反应往往不是看它的起伏而是问“它的频谱长什么样” 这就像医生看心电图不仅要看心跳的节奏时域更要分析其频率成分是否健康。Matlab 中的 FFT快速傅里叶变换就是我们完成这个“频谱体检”最核心、最常用的“听诊器”。这个项目标题“Matlab的信号频谱分析——FFT变换”看似简单但它背后串联的是从理论到实践从理想模型到工程落地的完整链条。我处理过从音频分析、振动诊断到通信解调的无数案例深知仅仅调用fft()函数是远远不够的如何正确设置参数、理解结果物理意义、避开频谱泄露和栅栏效应这些坑才是从“会用”到“精通”的关键。这篇文章我就以一个老工程师的视角拆解用 Matlab 做 FFT 频谱分析时那些必须掌握的细节和容易踩的坑让你拿到的频谱图不仅好看更要准确、可信。2. 核心原理与工程意义为什么是 FFT在深入代码之前我们必须搞清楚 FFT 到底是什么以及它为什么在工程领域不可替代。傅里叶变换的本质是告诉我们一个复杂的信号可以分解成许多不同频率、不同幅度、不同相位的正弦波的叠加。而 FFT是离散傅里叶变换DFT的一种快速算法它将计算复杂度从 O(N²) 降低到 O(N log N)这才使得在普通计算机上实时分析海量数据成为可能。2.1 FFT 与 DFT 的关系从“算盘”到“计算器”你可以把 DFT 理解成一种定义对于 N 个离散的时域数据点通过一套固定的复数乘加公式计算出 N 个离散的频域点。这套公式是精确的但计算量巨大就像用算盘做大规模乘法。FFT 则是一套巧妙的“心算”技巧它利用正弦函数的对称性和周期性将一个大 DFT 分解成多个小 DFT 的组合从而极大地减少了重复计算。在 Matlab 中我们用的fft(x)函数就是实现了这种高效算法的黑箱。但作为工程师我们不能只满足于黑箱至少要明白它的输入输出特性它假设你给它的 N 点数据是某个无限长周期信号的一个周期然后帮你算出这个假设周期信号的频谱。2.2 频谱的物理意义幅度谱、相位谱与功率谱FFT 直接输出的结果是复数包含实部和虚部。这对数学是完备的但对工程师不够直观。我们通常从中提取三种更有工程意义的谱幅度谱每个频率分量的强度大小。abs(fft(x))就能得到。这是最常用的用于查找信号中的主频、谐波。相位谱每个频率分量的初始相位。angle(fft(x))得到。在图像处理、通信系统同步和某些故障诊断中至关重要。功率谱幅度谱的平方或归一化的幅度平方。它反映信号功率在频域的分布在能量分析中常用单位是 dB。一个关键概念是单边谱和双边谱。由于实数信号的频谱具有共轭对称性FFT 输出的前半部分正频率和后半部分负频率是对称的。为了更直观我们通常只显示正频率部分并将除直流分量0Hz和奈奎斯特频率点外的幅度乘以2以补偿被折叠掉的负频率能量这就是单边幅度谱。这是工程展示的标准做法。3. 完整实操流程从数据导入到频谱图解读理论说再多不如动手做一遍。下面我以一个包含 50Hz 和 120Hz 正弦波的混合信号为例展示从零开始的完整分析流程并穿插解释每一个参数的选择依据。3.1 环境准备与信号合成首先我们合成一个干净的测试信号这有助于验证我们流程的正确性。% 1. 基本参数设置 Fs 1000; % 采样频率单位 Hz。根据香农采样定理要大于信号最高频率的2倍。 T 1/Fs; % 采样间隔 L 1500; % 信号长度点数。选择1500而非1024是为了演示非2的幂次时FFT依然高效。 t (0:L-1)*T; % 时间向量 % 2. 合成信号一个包含直流、基波和高频噪声的典型工程信号 f1 50; % 基频 50Hz f2 120; % 二次谐波 120Hz S 0.7*sin(2*pi*f1*t) sin(2*pi*f2*t); % 两个正弦波叠加 % 模拟实际情况加入直流偏移和高斯白噪声 X S 2.5*cos(2*pi*300*t); % 加入一个强干扰的高频分量300Hz X X 0.5*randn(size(t)); % 加入高斯白噪声模拟传感器噪声这里有几个经验点采样频率 Fs必须事先知道或设定。它决定了你能分析的最高频率Fs/2即奈奎斯特频率。如果信号本身有高于 Fs/2 的频率会发生混叠频谱将完全失真且不可逆。在数据采集阶段就必须用抗混叠滤波器。信号长度 L不一定非要选 2 的幂如1024。早期FFT算法对2的幂有优化但现代Matlab的fft函数对任意长度的输入都采用了高效算法。选择 L 更多考虑的是频率分辨率。3.2 执行 FFT 与频谱计算% 3. 执行FFT Y fft(X); % Y是复数长度也为L % 4. 计算双边谱和单边谱 P2 abs(Y/L); % 计算双边幅度谱并除以L进行归一化解释见下文 P1 P2(1:L/21); % 取前半部分包含直流和奈奎斯特频率点 P1(2:end-1) 2*P1(2:end-1); % 除直流和奈奎斯特点外幅度乘2得到单边谱 % 5. 构建频率轴 f Fs*(0:(L/2))/L; % 单边谱对应的频率向量从0Hz到Fs/2关键解释abs(Y/L)为什么除以 L这是为了幅度归一化。使得时域中幅度为 A 的正弦波在频域对应谱线的幅度也是 A。如果不除谱线幅度将是 A*L/2失去了物理意义。这是新手最常忽略的一步。P1(2:end-1) 2*P1(2:end-1)这就是生成单边谱的核心操作。因为能量平均分布在正负频率只显示正频率时需要将幅度加倍功率则需要乘4。直流分量0Hz和奈奎斯特频率Fs/2是特例它们没有对应的负频率成分所以不加倍。3.3 结果可视化与解读% 6. 绘制时域图和频谱图 figure(Position, [100, 100, 1200, 500]) % 子图1时域信号 subplot(1,2,1) plot(t(1:200), X(1:200)) % 只画前200个点便于观察波形 title(时域信号 (前200点)) xlabel(时间 (s)) ylabel(幅度) grid on % 子图2单边幅度频谱 subplot(1,2,2) plot(f, P1) title(单边幅度频谱) xlabel(频率 (Hz)) ylabel(|幅度|) xlim([0, Fs/2]) % 通常只显示到奈奎斯特频率 grid on % 标记出主要频率峰值 [~, locs] findpeaks(P1, MinPeakHeight, max(P1)*0.1); % 简单找峰 hold on plot(f(locs), P1(locs), ro, MarkerFaceColor, r) text(f(locs)5, P1(locs), cellstr(num2str(round(f(locs)))), VerticalAlignment,bottom) hold off运行这段代码你会得到清晰的时域/频域对比图。在频谱图上你应该能清晰地看到 50Hz 和 120Hz 处尖锐的谱线幅度大约为 0.7 和 1。300Hz 处也有一个明显的峰。而噪声则表现为整个频带底部的“毛刺”基底。注意findpeaks是信号处理工具箱里的一个非常实用的函数用于自动寻找局部极大值。在实际分析中我们经常需要用它来提取主频、谐波频率及其幅度。4. 高级议题与工程陷阱规避如果只是分析理想合成信号那太简单了。工程中的信号是“脏”的数据是有限的下面这些才是真正考验功力的地方。4.1 频谱泄露与加窗函数我们的 FFT 隐含了一个假设截取的那段数据长度L是周期信号的整数个周期。如果不是就会发生频谱泄露——能量从主频点“泄露”到旁边的频点导致谱线变宽、幅值不准旁边还会出现虚假的旁瓣。解决方案就是加窗。窗函数在时域上对信号两端进行平滑衰减减少截断带来的突变。Matlab 提供了hamming,hann,blackman,flattopwin等函数。% 加窗处理示例 win hann(L); % 生成汉宁窗转置成行向量 X_windowed X .* win; % 时域点乘窗函数 % 对加窗后的信号做FFT并修正窗函数带来的幅度损失 Y_win fft(X_windowed); P2_win abs(Y_win / (sum(win)/L)); % 关键归一化因子变为窗函数的平均高度 P1_win P2_win(1:L/21); P1_win(2:end-1) 2*P1_win(2:end-1); figure; plot(f, P1, b, LineWidth, 1.5); hold on; plot(f, P1_win, r--, LineWidth, 1.2); legend(无窗, 汉宁窗); title(加窗对频谱的影响); xlabel(频率 (Hz)); ylabel(|幅度|); grid on;关键点加窗后归一化因子不再是信号长度 L而是窗函数的和sum(win)或者更常用的是窗函数的平均高度sum(win)/L。这是因为窗函数削弱了信号两端的能量直接用 L 除会低估幅度。汉宁窗能有效抑制旁瓣但主瓣会略微变宽频率分辨率下降这是抑制泄露必须付出的代价。选择窗函数本质是在频谱分辨率主瓣宽度和频谱泄露旁瓣高度之间做权衡。4.2 频率分辨率与栅栏效应频率分辨率Δf Fs / L。它表示频谱图上相邻两个频点之间的间隔。如果你的信号中有两个频率分量 f1 和 f2只有当 |f1 - f2| Δf 时才能在频谱图上被区分开。增加 L采集更长时间的数据或降低 Fs在满足采样定理的前提下可以提高分辨率。栅栏效应是离散采样的固有特性FFT 只计算频率为 k*Δf (k0,1,2,...) 这些离散点上的频谱。如果信号的真实频率正好落在两个“栅栏”之间那么它的能量就会分散到相邻的频点上即使没有泄露幅值测量也会不准。解决方案增加数据长度 L这是最根本的方法直接提高分辨率。使用高分辨率谱估计方法如 Chirp-Z 变换可以对特定频段进行“细看”。补零在信号末尾补零后再做 FFTfft(x, N)其中 N length(x)。这不能提高真实的频率分辨率但可以通过对频谱进行插值让曲线更光滑有助于更精确地通过曲线拟合来定位峰值频率是一种“视觉增强”手段。% 演示栅栏效应和补零 L_short 128; % 短数据分辨率低 t_short (0:L_short-1)*T; X_short 0.7*sin(2*pi*50*t_short) sin(2*pi*55*t_short); % 两个很近的频率 % 不补零 Y1 fft(X_short); P1_short 2*abs(Y1(1:L_short/21)/L_short); f_short Fs*(0:(L_short/2))/L_short; % 补零到1024点 N_fft 1024; Y2 fft(X_short, N_fft); P1_zpad 2*abs(Y2(1:N_fft/21)/L_short); % 归一化仍用原数据长度 f_zpad Fs*(0:(N_fft/2))/N_fft; figure; subplot(2,1,1); stem(f_short, P1_short, b, LineWidth, 1.5); % 用 stem 更显离散性 title([短数据 (L, num2str(L_short), )分辨率低栅栏效应明显]); xlabel(频率 (Hz)); ylabel(|幅度|); grid on; xlim([40, 70]); subplot(2,1,2); plot(f_zpad, P1_zpad, r-, LineWidth, 1.2); title([补零到, num2str(N_fft), 点频谱插值更光滑但分辨率未变]); xlabel(频率 (Hz)); ylabel(|幅度|); grid on; xlim([40, 70]);从图中可以清晰看到短数据时两个频率峰混叠在一起补零后谱线变密能更好地描绘出包的形状但两个峰依然无法分开证明物理分辨率未变。4.3 平均与平滑从瞬时谱到统计谱对于平稳随机信号如噪声或为了抑制分析中的随机波动我们通常不直接分析一段数据的频谱而是计算平均功率谱密度。Welch 方法这是工程上的标准方法。它将长数据分段可重叠对每一段加窗并计算周期图单个段的功率谱最后对所有段的周期图求平均。Matlab 中的pwelch函数实现了它。[pxx, f_welch] pwelch(X, hann(256), 128, 1024, Fs); % 窗长256重叠128FFT点数1024 figure; plot(f_welch, 10*log10(pxx)); % 以dB为单位绘制 title(Welch方法估计的平均功率谱密度 (PSD)); xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); grid on;Welch 方法通过平均显著平滑了频谱降低了方差更能反映信号的统计特性特别适合分析噪声和随机振动信号。滑动平均平滑对于已经计算出的幅度谱也可以进行频域平滑比如使用移动平均滤波器。但这是一种后处理会损失频率分辨率。P1_smooth movmean(P1, 5); % 5点移动平均5. 实战案例电机振动信号分析让我们用一个更接近实战的场景来串联以上所有知识。假设我们采集了一段电机轴承的振动加速度信号采样频率 Fs 10 kHz数据长度 N 20000 点即2秒数据。我们怀疑轴承存在故障其故障特征频率约为 120 Hz。% 模拟电机振动信号 (包含故障频率和宽带噪声) Fs_motor 10000; t_motor (0:19999)/Fs_motor; f_fault 120; % 假设的故障特征频率 vibration 1.5 * sin(2*pi*f_fault*t_motor) ... % 故障特征 0.3 * sin(2*pi*2*f_fault*t_motor) ... % 二次谐波 0.1 * sin(2*pi*3*f_fault*t_motor) ... % 三次谐波 randn(size(t_motor)); % 强烈的随机振动噪声 % 1. 直接FFT分析效果可能不佳 L_motor length(vibration); Y_raw fft(vibration); P1_raw 2*abs(Y_raw(1:L_motor/21)/L_motor); f_motor Fs_motor*(0:(L_motor/2))/L_motor; % 2. 使用Welch方法估计PSD突出周期性成分 [pxx, f_welch] pwelch(vibration, hann(2048), 1024, 4096, Fs_motor); % 绘图对比 figure(Position, [50, 50, 1400, 600]); subplot(1,2,1); plot(f_motor, 10*log10(P1_raw.^2)); % 将幅度谱转为功率谱粗略对比 title(直接FFT得到的频谱 (dB)); xlabel(频率 (Hz)); ylabel(功率 (dB)); grid on; xlim([0, 500]); % 标记故障频率及其谐波 hold on; plot([f_fault, 2*f_fault, 3*f_fault], [-20, -30, -35], rv, MarkerFaceColor, r); hold off; subplot(1,2,2); plot(f_welch, 10*log10(pxx)); title(Welch方法估计的功率谱密度 (PSD)); xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); grid on; xlim([0, 500]); hold on; plot([f_fault, 2*f_fault, 3*f_fault], [max(10*log10(pxx))-10, -45, -50], rv, MarkerFaceColor, r); hold off; legend(PSD, 故障频率点);在这个案例中直接 FFT 的频谱被强大的背景噪声淹没故障频率的峰值并不明显。而经过 Welch 平均后的 PSD 图噪声基底变得平坦120Hz、240Hz、360Hz 处的故障特征频率及其谐波清晰地凸现出来这对于故障诊断具有决定性意义。6. 常见问题与调试技巧实录在实际操作中你一定会遇到各种奇怪的现象。下面是我总结的一些典型问题及排查思路。问题现象可能原因排查与解决方法频谱幅值不对远大于或小于预期未进行幅度归一化。FFT结果未除以信号长度L。计算幅度谱时务必使用abs(fft(x)/N)其中N是参与FFT运算的数据点数。频谱图在中间频率出现对称的“镜像”峰错误地绘制了双边谱。对于实数信号这是正常现象但通常我们只看单边谱。确保你只取了FFT结果的前半部分N/21点并对幅度进行了乘2处理直流和奈奎斯特点除外。单一频率的正弦波频谱却是一个很宽的“包”频谱泄露。信号截取长度不是信号周期的整数倍。加窗处理。使用hann或hamming窗。同时尽量采集更长时间的数据使截取长度接近周期的整数倍。已知信号频率为f0但频谱峰值在f0旁边栅栏效应。f0 不在频率分辨率的整数倍上。1.增加数据长度L以提高分辨率。2. 在数据后补零并进行FFT插值然后通过抛物线插值或寻找最大值点来更精确地估计真实频率。高频部分出现不应该有的低频成分频谱混叠。信号中包含高于奈奎斯特频率Fs/2的成分。这是硬件问题无法通过软件完全修复。必须在ADC采样前使用模拟抗混叠滤波器将高于Fs/2的频率成分滤除。检查你的采样率是否足够。功率谱密度PSD的计算结果单位不对对PSD的定义和归一化理解有误。pwelch默认返回的是单边PSD单位是x^2/Hz。理解pwelch的输出。如果要转换为 dB使用10*log10(pxx)。确保你的窗函数参数设置合理pwelch内部已经考虑了窗函数的能量归一化。findpeaks找不到正确的峰值或找到太多杂峰阈值设置不当。噪声基底过高。使用findpeaks(P1, MinPeakHeight, threshold)设置绝对阈值或MinPeakProminence设置最小峰凸起度。通常可以先估算噪声水平将阈值设为噪声水平的3-5倍。处理大量数据时FFT速度慢数据长度过长或循环中多次调用FFT。1. 确保数据长度是许多小素数的乘积Matlab的FFT对此优化最好。2. 考虑使用分段处理或实时频谱分析技术。3. 对于固定长度的FFT可以预先计算旋转因子。几个私房调试技巧从简单信号开始验证任何新的频谱分析流程先用一个幅度、频率已知的纯净正弦波测试。确保频谱图上在正确位置出现一根干净、幅度正确的谱线。这是检验你流程是否正确的“试金石”。关注频率轴很多错误源于频率轴构建错误。反复检查f Fs*(0:(N/2))/N这个公式。确保你的频率向量长度与单边谱数据长度完全一致。理解fft函数的第二个参数 NY fft(X, N)指定了进行 N 点 FFT。如果 N length(X)会自动补零如果 N length(X)会截断 X。这个特性在需要固定长度FFT或进行补零时非常有用。使用fftshift可视化fftshift可以将FFT输出的零频分量移动到频谱中心这对于观察以零频对称的信号如基带信号很方便。但注意这之后频率轴也需要相应调整。保存中间结果在编写复杂的分析脚本时将关键的中间变量如加窗后的信号、原始的FFT复数结果、频率向量保存下来或单独绘图检查能帮你快速定位问题出在哪一步。最后记住频谱分析是一门“艺术”需要在分辨率、精度、速度和平滑度之间根据具体应用做取舍。没有一种设置能通吃所有场景。多动手多对比用已知信号去验证你的流程你的“频谱直觉”就会慢慢建立起来。当你拿到一段陌生的信号能迅速在脑海中勾勒出它大致的频谱模样并知道用什么工具和方法去验证时你就真正掌握了这门技能。
)
