1. 从时域到频域为什么我们需要频谱分析作为一名在信号处理领域摸爬滚打了十多年的工程师我处理过从电机振动到无线通信的各种信号。我敢说至少有80%的信号问题在时域里看就是一团乱麻但切换到频域视角一切就豁然开朗了。这就像你听一首交响乐在时域里你听到的是所有乐器声音混合在一起、随时间变化的波形复杂且难以分辨而在频域里你可以清晰地看到小提琴、大提琴、定音鼓各自在哪个音高频率上演奏以及它们的音量幅度有多大。Matlab就是我们进行这种“音乐分解”的绝佳工具。很多刚接触信号处理的工程师朋友一上来就想用Matlab的fft函数但往往被横坐标、纵坐标的单位和幅值搞得一头雾水画出来的图怎么看都不对劲。这很正常因为fft函数本身只是一个“计算器”它只负责完成数学变换而如何呈现一个物理意义上正确的频谱图需要我们自己去构建频率轴、处理幅值。今天我就结合一个最经典的正弦波例子把Matlab频域分析从原理到实操的“坑”都填平让你不仅能画出正确的频谱图更能理解每一个步骤背后的“为什么”。2. 核心概念与准备工作采样、频率与频谱在动手写代码之前我们必须先理清几个核心概念。这是避免后续所有迷惑的基础。2.1 采样连接模拟与数字世界的桥梁我们处理的信号无论是麦克风采集的声音还是传感器读取的电压最初都是连续的模拟信号。计算机无法处理连续的数据必须对其进行采样将其变成离散的数字序列。这里有两个至关重要的参数采样频率Fs每秒采集多少个数据点单位是赫兹Hz。根据奈奎斯特采样定理为了无失真地还原信号Fs必须大于信号最高频率成分的两倍。比如你想分析一个最高频率为1000Hz的信号你的Fs至少需要2000Hz。采样点数N你总共采集了多少个数据点。它决定了频谱的频率分辨率。N越大频率分辨率越高能区分开更近的两个频率但计算量也越大。注意Fs决定了你能分析的最高频率Fs/2即奈奎斯特频率而N决定了你分析得有多“精细”。在实际项目中你需要根据信号特性和分析需求在这两者间权衡。2.2 频谱信号在频率上的“身份证”对一个信号做傅里叶变换FFT是其快速算法目的就是得到它的频谱。频谱告诉我们这个信号是由哪些不同频率的正弦波组成的以及这些正弦波的幅度和相位分别是什么。一个理想的单一频率正弦波sin(2πf0*t)其频谱应该是在f0和-f0处各有一根谱线。为什么有负频率这是数学上复数表示带来的在物理上f0和-f0共同代表了一个实际存在的频率为f0的振荡。对于实信号我们遇到的大多数信号其频谱总是关于0Hz共轭对称的这意味着负频率部分的信息是冗余的我们通常只关心正频率部分。2.3 Matlab环境与代码管理我强烈建议你为这个教程新建一个文件夹比如FFT_Tutorial。将你写的所有.m脚本文件和后面要定义的自定义函数文件都放在这里并在Matlab中将当前文件夹切换到此处。这能避免因路径问题导致的“函数未定义”错误也是工程化思维的良好起点。3. 从零开始生成一个时域信号并观察理论说得再多不如一行代码。我们从一个最简单的信号开始一个频率为4Hz幅度为2的正弦波。% 参数定义 fo 4; % 信号频率单位Hz A 2; % 信号幅度 Fs 100; % 采样频率单位Hz。远大于2*fo满足奈奎斯特定理。 Ts 1/Fs; % 采样间隔单位秒 duration 1; % 信号持续时间单位秒 % 生成时间向量t和信号y t 0:Ts: (duration - Ts); % 从0开始到 (1-Ts) 结束总共 duration*Fs 个点 n length(t); % 采样点数 N y A * sin(2*pi*fo*t); % 生成正弦信号 % 绘制时域波形 figure(Name, 时域正弦波, Position, [100, 100, 800, 400]); plot(t, y, b-, LineWidth, 1.5); xlabel(时间 (秒)); ylabel(幅度 y(t)); title(sprintf(时域正弦波: 频率 %.1f Hz, 幅度 %.1f, fo, A)); grid on; xlim([0, 1]); % 显示1秒的长度足以看到多个周期运行这段代码你会看到一个清晰的正弦波形。横轴是时间纵轴是幅度非常直观。我们的目标就是将这个时域图转换成下面这样的频域图——在4Hz的位置有一根高度为2的谱线注意由于FFT计算方式的差异实际峰值高度为A/2我们后面会详细解释并修正。4. 初探FFT为什么直接使用fft()会让人困惑现在我们对信号y进行FFT。这是新手最常踩的坑我们一起来看看直接使用会出什么问题。% 直接使用Matlab内置的fft函数 Y_raw fft(y); % 对信号y做FFT % 绘制FFT结果的幅度谱取绝对值 figure(Name, 直接使用FFT的结果, Position, [100, 100, 800, 400]); stem(abs(Y_raw), filled, MarkerSize, 3); xlabel(谱线索引 (k)); ylabel(|Y(k)|); title(直接使用 fft() 函数得到的幅度谱); grid on; xlim([0, n]); % 显示所有谱线你会得到一张令人困惑的图。横坐标是0到N-1的索引号k而不是我们期望的频率Hz。幅度值也非常大在索引4和索引96附近有两个很高的尖峰但它们的值大约是100而不是我们信号的幅度2。为什么会这样横坐标问题fft函数默认输出的横坐标是数字频率索引k范围是0到N-1。它对应的物理频率f需要通过公式计算f k * (Fs / N)。当k N/2时对应的是负频率部分f (k - N) * (Fs / N)。纵坐标问题fft计算的结果没有进行归一化。对于一个长度为N的序列其FFT结果的幅度是原始信号幅度的N倍。对于实信号能量会平均分布在正负频率上所以单个频率点的幅度会是A * N / 2。在我们的例子中A2,N100所以峰值约为2 * 100 / 2 100。所以这张图在数学上没错但它对工程师不友好。我们需要对它进行“后处理”才能得到物理意义清晰的频谱图。5. 构建专业频谱图自定义函数实战为了解决上述问题我习惯编写两个自定义函数。它们封装了频率轴生成、幅值归一化和频谱搬移的步骤让你一键生成想要的频谱。5.1 双边频谱函数centeredFFT双边频谱展示了完整的正负频率信息对于理解FFT的对称性和进行某些复杂运算如滤波有帮助。将以下代码保存为centeredFFT.m文件放在你的工作目录下。function [X, freq] centeredFFT(x, Fs) % CENTEREDFFT 计算并返回中心化双边频谱 % 输入 % x - 时域输入信号向量 % Fs - 采样频率 (Hz) % 输出 % X - 中心化后的复数频谱 % freq - 对应的频率轴向量 (Hz)范围从 -Fs/2 到 Fs/2近似 N length(x); % 获取信号长度 % 生成对称的频率轴向量 % 如果N是偶数频率点从 -N/2 到 N/2-1 % 如果N是奇数频率点从 -(N-1)/2 到 (N-1)/2 if mod(N, 2) 0 k -N/2 : (N/2 - 1); else k -(N-1)/2 : (N-1)/2; end T N / Fs; % 信号总时长 freq k / T; % 生成物理频率轴单位Hz % 计算FFT并进行处理 X fft(x) / N; % 关键一步除以N进行归一化使幅度具有物理意义 X fftshift(X); % 关键一步将零频率分量移动到频谱中心 end关键点解析X fft(x) / N这是幅度归一化。除以采样点数N后频谱的幅度就对应了原始信号中该频率正弦分量的实际幅度对于单频信号正负频率点幅度各为 A/2加起来为A。fftshift这个函数将FFT输出的前半部分正频率和后半部分负频率交换使得频率轴从[-Fs/2, Fs/2]的顺序排列零频率在中心。这对于观察对称的频谱非常直观。频率轴生成我们根据k和总时长T来生成频率轴freq。freq k / T是数字频率到物理频率的转换公式。现在使用这个函数% 使用自定义的centeredFFT函数 [Y_centered, freq_axis] centeredFFT(y, Fs); % 绘制双边频谱 figure(Name, 双边频谱 (centeredFFT), Position, [100, 100, 900, 400]); subplot(1,2,1); stem(freq_axis, abs(Y_centered), filled, MarkerSize, 4, LineWidth, 1); xlabel(频率 (Hz)); ylabel(幅度 |Y(f)|); title(双边幅度谱); grid on; xlim([-10, 10]); % 聚焦在-10Hz到10Hz范围 ylim([0, 1.2]); % 幅度应该在1附近 (A/2 2/2 1) % 为了更清晰绘制相位谱可选 subplot(1,2,2); stem(freq_axis, angle(Y_centered)*180/pi, filled, MarkerSize, 4); % 相位转换为度 xlabel(频率 (Hz)); ylabel(相位 (度)); title(双边相位谱); grid on; xlim([-10, 10]);运行后你会看到完美的双边频谱图。幅度谱在4Hz和-4Hz处各有一个峰值幅度为1即 A/2。相位谱在4Hz处约为 -90度-π/2-4Hz处约为 90度π/2这正是纯正弦波sin(2πft)的相位特性余弦波cos(2πft)的相位则在0度。5.2 单边频谱函数positiveFFT在绝大多数工程应用中我们只关心正频率部分因为实信号的频谱是共轭对称的负频率部分不提供新的信息。单边频谱更简洁也是仪器如频谱分析仪通常显示的格式。将以下代码保存为positiveFFT.m文件。function [X, freq] positiveFFT(x, Fs) % POSITIVEFFT 计算并返回单边频谱 % 输入 % x - 时域输入信号向量 % Fs - 采样频率 (Hz) % 输出 % X - 单边复数频谱仅正频率部分0Hz和奈奎斯特频率点特殊处理 % freq - 对应的正频率轴向量 (Hz)范围从 0 到 Fs/2 N length(x); k 0:N-1; % 索引从0开始 T N / Fs; freq k / T; % 此时freq范围是 0 到 Fs*(N-1)/N ≈ Fs % 计算FFT并归一化 X fft(x) / N; % 取前半部分正频率部分包括0Hz和可能的奈奎斯特频率点 cutOff floor(N/2) 1; % 例如 N100则取前51个点 (0到50) X X(1:cutOff); freq freq(1:cutOff); % 关键调整除0Hz和奈奎斯特频率点外其他点的幅度乘以2 % 这是因为能量分布在正负频率上单边谱需要将负频率的能量加回来。 X(2:end-1) 2 * X(2:end-1); % 索引2到end-1避开了0Hz和最后一个点如果Fs/2存在 % 如果N是偶数最后一个点对应Fs/2奈奎斯特频率这个点没有对应的负频率所以不乘2。 end关键点解析取前半部分cutOff floor(N/2) 1确保了无论N是奇数还是偶数我们都能正确取到从0Hz到奈奎斯特频率Fs/2附近的所有正频率点。幅度乘以2这是单边谱的核心修正。在双边谱中一个频率分量的能量幅度被平均分到了f和-f两个点上。当我们只显示正频率部分时需要将f点的幅度乘以2以代表该频率分量的总幅度。但有两个例外直流分量0Hz它没有对应的负频率所以不乘2。奈奎斯特频率点Fs/2当N为偶数时这个点存在且是实数它也没有对应的负频率或者说它自己就是正负频率的边界所以也不乘2。现在使用这个函数% 使用自定义的positiveFFT函数 [Y_positive, freq_positive] positiveFFT(y, Fs); % 绘制单边频谱 figure(Name, 单边频谱 (positiveFFT), Position, [100, 100, 800, 400]); stem(freq_positive, abs(Y_positive), filled, MarkerSize, 5, LineWidth, 1.2, Color, [0.85, 0.33, 0.10]); xlabel(频率 (Hz)); ylabel(幅度 |Y(f)|); title(单边幅度谱 (幅度已校正)); grid on; xlim([0, 20]); % 显示到20Hz足以 ylim([0, 2.2]); % 幅度应该在2附近 (A 2) % 在峰值处添加标注 hold on; [peak_mag, peak_idx] max(abs(Y_positive)); plot(freq_positive(peak_idx), peak_mag, ro, MarkerSize, 10, LineWidth, 2); text(freq_positive(peak_idx)0.5, peak_mag, sprintf(峰值: %.2f Hz, 幅度: %.2f, freq_positive(peak_idx), peak_mag), ... VerticalAlignment, bottom, FontWeight, bold); hold off;这次你将得到一张非常干净的频谱图。只在4Hz处有一个峰值并且其幅度为2这正是我们原始正弦信号的幅度A。这张图才是一个工程师期望看到的、具有直接物理意义的频谱。6. 应对复杂信号多频信号与噪声的频谱分析现实世界的信号很少是单一频率的。它们通常是多个频率的叠加并且混杂着噪声。让我们用自定义的函数来分析一个更复杂的信号。% 生成一个包含多个频率分量和噪声的复杂信号 Fs 1000; % 采样率 1kHz t 0:1/Fs:1-1/Fs; % 1秒时长 N length(t); % 信号成分10Hz幅度350Hz幅度1.5120Hz幅度1 signal_clean 3*sin(2*pi*10*t) 1.5*cos(2*pi*50*t) 1*sin(2*pi*120*t); % 添加随机噪声 noise 0.5 * randn(size(t)); y_complex signal_clean noise; % 绘制时域信号 figure(Name, 复杂时域信号, Position, [50, 50, 1200, 400]); subplot(1,2,1); plot(t, y_complex, b-, LineWidth, 0.5); xlabel(时间 (秒)); ylabel(幅度); title(含噪声的复杂时域信号); grid on; xlim([0, 0.5]); % 只看前0.5秒否则太密 % 使用positiveFFT进行频域分析 [Y_comp, freq_comp] positiveFFT(y_complex, Fs); % 绘制单边频谱 subplot(1,2,2); stem(freq_comp, abs(Y_comp), filled, MarkerSize, 3, LineWidth, 0.8); xlabel(频率 (Hz)); ylabel(幅度); title(复杂信号的单边幅度谱); grid on; xlim([0, 200]); % 显示到200Hz覆盖所有信号频率 ylim([0, 3.5]); % 标记出主要的频率成分 hold on; freq_of_interest [10, 50, 120]; for f freq_of_interest % 找到最接近的频点索引 [~, idx] min(abs(freq_comp - f)); mag abs(Y_comp(idx)); plot(freq_comp(idx), mag, ro, MarkerSize, 8, LineWidth, 2); text(freq_comp(idx)5, mag, sprintf(%d Hz\n%.2f, f, mag), ... VerticalAlignment, bottom, FontWeight, bold); end hold off;在这张频谱图上你可以清晰地看到在10Hz、50Hz和120Hz处有三个突出的峰值它们的幅度分别接近3、1.5和1。尽管时域波形因为噪声的加入显得杂乱无章但频域分析成功地剥离出了信号的核心频率成分。背景中起伏的“毛刺”就是宽带噪声的频谱体现。7. 高级话题与实战避坑指南掌握了基础操作后我们来看看在实际工程中会遇到哪些问题以及如何解决。7.1 频谱泄露与加窗在上面的例子中我们的信号频率4Hz恰好是频率分辨率Fs/N 1 Hz的整数倍所以能量完美地集中在了一个频点上。但如果不是整数倍呢% 频谱泄露演示 Fs 100; t 0:1/Fs:1-1/Fs; % 1秒N100 fo 4.5; % 频率不是1Hz的整数倍 y_leak sin(2*pi*fo*t); [Y_leak, freq_leak] positiveFFT(y_leak, Fs); figure(Name, 频谱泄露现象); stem(freq_leak, abs(Y_leak), filled); xlabel(频率 (Hz)); ylabel(幅度); title(sprintf(频谱泄露: 信号频率 %.1f Hz, fo)); grid on; xlim([0, 10]);你会发现峰值不再尖锐能量“泄露”到了旁边的频点上主瓣变宽旁瓣出现。这是因为我们截取了一段有限长的信号1秒相当于用一个矩形窗去乘无限长的正弦波在频域上引入了卷积效应。解决方案加窗。在FFT前对时域信号乘以一个窗函数如汉宁窗、汉明窗可以抑制频谱泄露降低旁瓣电平代价是主瓣会稍微变宽。% 加窗处理 window hanning(length(y_leak)); % 生成汉宁窗 y_windowed y_leak .* window; % 加窗 [Y_win, freq_win] positiveFFT(y_windowed, Fs); figure(Name, 加窗效果对比); subplot(2,1,1); stem(freq_leak, abs(Y_leak), b); hold on; stem(freq_win, abs(Y_win), r, LineWidth, 1.5); hold off; xlabel(频率 (Hz)); ylabel(幅度); legend(无窗矩形窗, 汉宁窗); title(幅度谱对比); grid on; xlim([0, 10]); subplot(2,1,2); plot(t, y_leak, b-, LineWidth, 1); hold on; plot(t, y_windowed, r-, LineWidth, 1.5); hold off; xlabel(时间 (秒)); ylabel(幅度); legend(原始信号, 加汉宁窗后信号); title(时域波形对比); grid on;加窗后频谱的旁瓣两侧的“小鼓包”被显著抑制频谱看起来更“干净”但主峰确实变宽了。在需要精确测量频率和幅度的场合如振动分析、音频分析选择合适的窗函数是必须的。7.2 频率分辨率与补零频率分辨率Δf Fs / N。要提高分辨率让谱线更密要么降低Fs不能低于奈奎斯特率要么增加N采集更长时间的数据。如果数据长度固定可以通过补零来增加FFT点数使频谱图看起来更平滑但这不能提高真正的频率分辨率只是对已有频谱的插值。% 补零演示 Fs 100; t 0:1/Fs:0.5-1/Fs; % 只采集0.5秒N50 y_short sin(2*pi*4*t); % 不补零 N_original length(y_short); [Y_orig, f_orig] positiveFFT(y_short, Fs); % 补零到256点 N_zeropad 256; y_zeropadded [y_short, zeros(1, N_zeropad - N_original)]; [Y_zp, f_zp] positiveFFT(y_zeropadded, Fs); figure(Name, 补零效果); subplot(2,1,1); stem(f_orig, abs(Y_orig), filled); title(sprintf(原始数据N%dΔf%.2f Hz, N_original, Fs/N_original)); xlabel(频率 (Hz)); ylabel(幅度); grid on; xlim([0, 20]); subplot(2,1,2); plot(f_zp, abs(Y_zp), r-, LineWidth, 1); title(sprintf(补零后N%d谱线更密插值, N_zeropad)); xlabel(频率 (Hz)); ylabel(幅度); grid on; xlim([0, 20]);补零后的频谱曲线更平滑但两个频率非常接近的信号在补零前后能否被分辨开的能力是一样的。真正的分辨率只由原始数据长度决定。7.3 幅度校准与功率谱我们之前的positiveFFT函数进行了幅度校准使得正弦波的峰值等于其幅度。但在有些场合我们关心的是功率幅度平方。Matlab也提供了pwelch、periodogram等函数来直接估计功率谱密度PSD。了解幅度谱和功率谱的关系很重要幅度谱abs(Y)单位与原始信号相同如伏特V。功率谱(abs(Y).^2)单位是原始信号的平方如V²。功率谱密度PSD(abs(Y).^2) / (Fs * 窗函数能量)单位是V²/Hz用于描述功率在频率上的分布密度在分析噪声时特别有用。8. 常见问题与排查技巧实录在我多年的调试经历中下面这些问题几乎每个人都遇到过。我把它们整理成表方便你快速排查。问题现象可能原因解决方案与检查点频谱峰值幅度不对不是A或A/21. 未做归一化 (/N)。2. 单边谱未对非直流分量乘2。3. 信号中含有直流偏移。1. 检查FFT后是否除以了采样点数N。2. 检查单边谱函数是否正确处理了乘2的逻辑避开0Hz和Fs/2。3. 绘制时域信号检查均值是否为零可使用y y - mean(y)去除直流。频谱图看起来是镜像的左右对称但中心不在0Hz使用了fftshift但频率轴没跟着变或者反之。确保fftshift操作和频率轴freq的生成是匹配的。centeredFFT函数中先计算freq再fftshift(X)顺序不能错。频率轴显示不对峰值位置不是预期的Hz值频率轴计算公式错误。确认公式freq k * (Fs / N)。在双边谱中k的范围是-N/2 : N/2-1在单边谱中是0 : floor(N/2)。频谱图非常杂乱看不到明显峰值1. 信号信噪比太低噪声淹没信号。2. 频谱泄露严重。3. FFT点数太少频率分辨率太低。1. 尝试对信号进行平均如果信号是周期性的或使用带通滤波。2. 尝试加窗如汉宁窗。3. 增加数据采集时间增大N以提高分辨率。在Fs/2附近出现奇怪的峰值1. 信号中含有高于Fs/2的频率成分混叠。2. 可能是频谱泄露的旁瓣。1.这是严重问题检查信号最高频率是否满足奈奎斯特定理。在采样前必须使用抗混叠滤波器。2. 尝试加窗观察峰值是否消失或减弱。单边谱在0Hz或Fs/2处幅度异常高单边谱幅度修正时错误地对直流或奈奎斯特频率点乘了2。检查positiveFFT函数中乘2的逻辑确保X(1)直流和X(end)当N为偶数时的Fs/2点没有被乘以2。最后分享一个我调试时的小习惯永远先从一个已知的、干净的信号开始。就像我们教程里用的4Hz正弦波。先确保对这个简单信号的频谱分析是完全正确的峰值位置、幅度、相位然后再把复杂的真实信号代入。这能帮你快速定位问题是出在信号本身还是你的分析代码上。频域分析是信号处理工程师的“眼睛”把这套流程练熟无论是调试电路、分析振动还是处理通信数据你都能游刃有余。