Matlab频谱分析实战：从FFT原理到工程应用避坑指南

1. 从时域到频域：为什么我们需要频谱分析

作为一名在信号处理领域摸爬滚打了十多年的工程师，我处理过从电机振动到无线通信的各种信号。我敢说，至少有80%的信号问题，在时域里看就是一团乱麻，但切换到频域视角，一切就豁然开朗了。这就像你听一首交响乐，在时域里你听到的是所有乐器声音混合在一起、随时间变化的波形，复杂且难以分辨；而在频域里，你可以清晰地看到小提琴、大提琴、定音鼓各自在哪个音高（频率）上演奏，以及它们的音量（幅度）有多大。Matlab，就是我们进行这种“音乐分解”的绝佳工具。

很多刚接触信号处理的工程师朋友，一上来就想用Matlab的fft函数，但往往被横坐标、纵坐标的单位和幅值搞得一头雾水，画出来的图怎么看都不对劲。这很正常，因为fft函数本身只是一个“计算器”，它只负责完成数学变换，而如何呈现一个物理意义上正确的频谱图，需要我们自己去构建频率轴、处理幅值。今天，我就结合一个最经典的正弦波例子，把Matlab频域分析从原理到实操的“坑”都填平，让你不仅能画出正确的频谱图，更能理解每一个步骤背后的“为什么”。

2. 核心概念与准备工作：采样、频率与频谱

在动手写代码之前，我们必须先理清几个核心概念。这是避免后续所有迷惑的基础。

2.1 采样：连接模拟与数字世界的桥梁

我们处理的信号，无论是麦克风采集的声音，还是传感器读取的电压，最初都是连续的模拟信号。计算机无法处理连续的数据，必须对其进行采样，将其变成离散的数字序列。

这里有两个至关重要的参数：

1. 采样频率Fs：每秒采集多少个数据点，单位是赫兹（Hz）。根据奈奎斯特采样定理，为了无失真地还原信号，Fs必须大于信号最高频率成分的两倍。比如，你想分析一个最高频率为1000Hz的信号，你的Fs至少需要2000Hz。
2. 采样点数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。

为什么会这样？

1. 横坐标问题：fft函数默认输出的横坐标是数字频率索引k，范围是0到N-1。它对应的物理频率f需要通过公式计算：f = k * (Fs / N)。当k > N/2时，对应的是负频率部分（f = (k - N) * (Fs / N)）。
2. 纵坐标问题：fft计算的结果没有进行归一化。对于一个长度为N的序列，其FFT结果的幅度是原始信号幅度的N倍。对于实信号，能量会平均分布在正负频率上，所以单个频率点的幅度会是A * N / 2。在我们的例子中，A=2,N=100，所以峰值约为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; % 例如 N=100，则取前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，以代表该频率分量的总幅度。但有两个例外：
  1. 直流分量（0Hz）：它没有对应的负频率，所以不乘2。
  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（幅度3），50Hz（幅度1.5），120Hz（幅度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秒，N=100 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秒，N=50 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. 常见问题与排查技巧实录

在我多年的调试经历中，下面这些问题几乎每个人都遇到过。我把它们整理成表，方便你快速排查。

最后分享一个我调试时的小习惯：永远先从一个已知的、干净的信号开始。就像我们教程里用的4Hz正弦波。先确保对这个简单信号的频谱分析是完全正确的（峰值位置、幅度、相位），然后再把复杂的真实信号代入。这能帮你快速定位问题是出在信号本身，还是你的分析代码上。频域分析是信号处理工程师的“眼睛”，把这套流程练熟，无论是调试电路、分析振动还是处理通信数据，你都能游刃有余。