简介：快速傅里叶变换（FFT）是信号处理中将时域信号转换为频域的关键技术，广泛应用于频谱分析和功率谱分析。本文以MATLAB为工具，深入讲解如何利用 fft 函数进行信号的频域变换，涵盖从数据读取、加窗处理到频谱与功率谱计算的完整流程。通过 fft.m 脚本实例，介绍采样率、频率分辨率、窗口函数应用及频率轴校正等核心概念，并结合 plot 、 specgram 等函数实现可视化，帮助用户掌握噪声分析、滤波设计和系统识别中的关键技术。本项目适用于信号处理初学者及工程实践者，具备较强的可操作性与实用价值。
傅里叶变换的工程艺术：从理论到MATLAB实战
你有没有遇到过这样的情况？采集了一段振动信号，兴冲冲地调用 fft() 函数，结果出来的频谱图却像一团模糊的毛线球——主峰歪斜、旁瓣乱飞、还冒出几个莫名其妙的“幽灵频率”😱。别急，这可不是你的代码写错了，而是我们忽略了FFT背后的那些 隐藏规则 。
今天咱们就来一场深度拆解之旅，不讲教科书式的定义堆砌，而是以一个工程师的视角，带你真正搞懂：
- 为什么补零不能提高分辨率？
- 加窗到底是救星还是副作用制造者？
- 如何让频谱图既准确又好看？
准备好了吗？Let’s go！🚀
🔧 FFT的本质：不是魔法，是数学巧劲
先说个真相： fft() 函数本身并不创造信息 ，它只是把已有的数据用另一种方式“重新排列组合”。就像把一堆杂乱的衣服按颜色分类，衣服总量没变，但你看得更清楚了。
DFT → FFT：从 $ O(N^2) $ 到 $ O(N \log N) $ 的飞跃
设想你要分析一段1024点的信号。如果老老实实按DFT公式算：
$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi kn/N}
$$
每个频率点都要做1024次复数乘法和加法，总共要算 $ 1024^2 = 1,048,576 $ 次运算。这在实时系统里简直是灾难 ❌。
而FFT呢？利用旋转因子 $ W_N = e^{-j2\pi/N} $ 的周期性和对称性，把大问题拆成小问题递归处理。最经典的 库利-图基算法（Cooley-Tukey） 就是把序列按奇偶分组，分别做FFT，再通过“蝶形运算”合并：
$$
\begin{aligned}
X[k] &= E[k] + W_N^k O[k] \
X[k+N/2] &= E[k] - W_N^k O[k]
\end{aligned}
$$
其中：
- $ E[k] $：偶数项的DFT
- $ O[k] $：奇数项的DFT
这个过程可以用下面这张流程图直观展示：
graph TD
A[N点输入序列x[n]] --> B{N == 1?}
B -- Yes --> C[返回x[0]]
B -- No --> D[分离奇偶索引]
D --> E[对偶数项做FFT]
D --> F[对奇数项做FFT]
E --> G[计算E[k]]
F --> H[计算O[k]]
G & H --> I[执行蝶形运算]
I --> J[输出X[k]]
看到了吧？这就是典型的“分而治之”策略。MATLAB里的 fft() 内部正是基于这类高效算法实现的，而且聪明得很——即使你的数据长度不是2的幂，它也能自动选择最优路径（比如混合基FFT或Bluestein算法），确保性能最大化 💪。
📌 小贴士 ：虽然速度快了上百倍，但FFT的数学本质仍是DFT。所有关于采样定理、频谱泄漏、分辨率限制的问题，一个都不会少！
⚙️ MATLAB中 fft() 的真实打开方式
你以为 Y = fft(X) 很简单？其实里面藏着不少“潜规则”。
输入输出那些事
| 参数 | 含义 |
|---|---|
X | 可以是向量、矩阵或多维数组 |
n | 指定FFT点数；若未指定则取 length(X) |
dim | 沿哪个维度变换，默认为第一个非单例维 |
举个例子：
X_matrix = [x; x]; % 2×1000 矩阵
Y_matrix = fft(X_matrix); % 默认对每列独立做FFT
这对多通道信号（如三轴加速度计、EEG脑电）特别有用，不需要循环就能并行处理。
复数输出的秘密
执行完 fft() 后得到的是复数数组，形式为 $ a + bj $，其中：
- 幅度： abs(Y) → 表示各频率成分的能量大小
- 相位： angle(Y) → 表示初始相位角
对于实信号，还有一个重要性质： 共轭对称性
$$
Y[k] = Y^*[N-k], \quad k=1,\dots,N-1
$$
这意味着负频率部分完全由正频率镜像而来。所以通常只看前半段就够了，省时又省心 ✅。
🎯 频率分辨率陷阱：补零真的能看清细节吗？
这是新手最容易掉进去的坑之一！很多人以为：“我把FFT点数从1000拉到10000，分辨率不就提升了10倍？” —— 错！大错特错！🚫
先搞清两个概念
| 名称 | 公式 | 决定因素 |
|---|---|---|
| 频率步长（Δf） | $ f_s / N $ | 采样率和FFT点数 |
| 真实分辨率 | $ 1/T $ | 信号持续时间 |
只有当 $ N = f_s \cdot T $ 时，两者才一致。否则你就只是在“插值”，相当于把一张低清图片放大，像素多了，但依然模糊 😒。
实验说话：三种情况对比
我们构造一个50Hz正弦波，分别测试以下三种情形：
fs = 1000; % 采样率
t_short = 0:1/fs:0.1-1/fs; % 仅0.1秒
x_short = sin(2*pi*50*t_short);
% 情况1：直接100点FFT
X1 = fft(x_short, 100);
f1 = (0:99)*fs/100;
% 情况2：补零至1000点
x_padded = [x_short, zeros(1, 900)];
X2 = fft(x_padded);
f2 = (0:999)*fs/1000;
% 情况3：真实采集1秒信号
t_long = 0:1/fs:1-1/fs;
x_long = sin(2*pi*50*t_long);
X3 = fft(x_long, 1000);
f3 = f2;
画出来一看就知道区别有多大：
figure;
subplot(3,1,1); plot(f1, abs(X1)); title('100点原始FFT');
subplot(3,1,2); plot(f2, abs(X2)); title('补零至1000点');
subplot(3,1,3); plot(f3, abs(X3)); title('真实1000点FFT');
👀 观察重点：
- 子图1：步长10Hz，勉强能看到峰值；
- 子图2：步长变1Hz，谱线密了，但主瓣宽度没变 → 还是分不清50Hz和51Hz；
- 子图3：因观测时间延长到1秒，主瓣明显收窄，这才具备真正分辨能力！
🎯 结论 ：
想提升分辨率？请耐心等信号多跑一会儿，而不是靠补零“作弊”。
| 场景 | 推荐做法 |
|---|---|
| 提升频率分辨率 | 延长采集时间 $ T $ |
| 改善频谱平滑度 | 适当补零（如2~4倍原长） |
| 实时系统资源紧张 | 使用重叠FFT或Welch方法 |
| 区分密集频率成分 | 结合MUSIC、ESPRIT等超分辨算法 |
🛠 手动还原DFT：理解比调包更重要
虽然 fft() 很快，但我们不妨自己动手实现一次DFT，看看底层发生了什么。
function X = my_dft(x)
N = length(x);
X = zeros(1, N);
for k = 0:N-1
for n = 0:N-1
X(k+1) = X(k+1) + x(n+1) * exp(-1i*2*pi*k*n/N);
end
end
end
这段代码逻辑清晰，但也暴露了问题：双重循环导致复杂度高达 $ O(N^2) $。我们来测一下性能：
| 方法 | 运算复杂度 | 1024点耗时（ms） | 是否实用 |
|---|---|---|---|
| DFT（手动循环） | $ O(N^2) $ | ~80 | ❌ 仅用于学习 |
| FFT（内置） | $ O(N \log N) $ | ~0.1 | ✅ 工程首选 |
差距超过800倍！这也说明了为什么工程实践中必须使用FFT。
不过，手动实现也有好处。比如我们可以用矩阵形式验证结果正确性：
$$
\mathbf{X} = \mathbf{W} \cdot \mathbf{x}
$$
其中 $ \mathbf{W} $ 是DFT矩阵，元素为 $ W_{kn} = e^{-j2\pi kn/N} $。
N = 8;
n = 0:N-1;
k = n';
W = exp(-1i*2*pi*k*n/N);
x = randn(1, N);
X_mat = W * x';
X_builtin = fft(x);
max_diff = max(abs(X_mat - X_builtin'));
fprintf('矩阵法与fft最大偏差: %.2e\n', max_diff); % 通常 < 1e-14
误差极小，说明两种方法高度一致 👏。
🌀 实信号 vs 复信号：频谱对称之美
你知道吗？实信号和复信号的频谱长得完全不同！
实验演示
% 实信号
x_real = cos(2*pi*100*(0:1/1000:1));
X_fft_real = fft(x_real);
% 复信号（解析信号）
x_complex = exp(1i*2*pi*100*(0:1/1000:1));
X_fft_comp = fft(x_complex);
N = length(x_real);
f = (-N/2:N/2-1)*(1000/N);
figure;
subplot(2,1,1);
plot(f, fftshift(abs(X_fft_real)));
title('实信号频谱（共轭对称）');
subplot(2,1,2);
plot(f, fftshift(abs(X_fft_comp)));
title('复信号频谱（非对称）');
📊 看图说话：
- 实信号：左右对称，负频率是正频率的镜像；
- 复信号：能量集中在正频带，适合通信系统中的IQ调制解调。
💡 这个性质在雷达、声纳、生物医学等领域至关重要。例如，在脉冲多普勒雷达中，正是利用复信号频谱的非对称性来判断目标是靠近还是远离。
📏 频率分辨率建模：参数之间的博弈
再来梳理一遍关键参数的关系：
$$
\Delta f = \frac{f_s}{N} = \frac{1}{T}, \quad T = \frac{N}{f_s}
$$
| 参数 | 调整影响 | 工程建议 |
|---|---|---|
| $ f_s $ | 提高可测最高频率 | 至少为信号带宽的2倍 |
| $ T $ | 决定最小可分辨间隔 | 尽可能延长 |
| $ N $ | 影响FFT点数选择 | 应匹配 $ f_s \cdot T $ |
📌 举个实际例子：你想区分50Hz和51Hz两个频率，那至少需要：
$$
T \geq \frac{1}{1\,\text{Hz}} = 1\,\text{秒}
$$
也就是说，采集时间不能少于1秒，否则无论你怎么补零都无济于事。
🧼 信号预处理全流程：别让脏数据毁了整个分析
现实中的信号哪有那么干净？直流偏置、趋势漂移、噪声干扰、异常跳变……这些都会让你的频谱变得面目全非。所以我们必须建立一套完整的预处理流水线。
数据导入：从文件到变量
常见格式处理方式如下：
CSV/TXT 文件读取
data = readmatrix('sensor_data.csv'); % 忽略标题
time = data(:,1);
signal = data(:,2);
% 或保留元信息
tbl = readtable('sensor_data.csv');
time_tbl = tbl.Time;
signal_tbl = tbl.Voltage;
MAT 文件加载
load('exp_data.mat'); % 加载所有变量
S = load('exp_data.mat', 'vib_sig', 'fs'); % 按需加载
signal = S.vib_sig;
fs = S.fs;
批量处理框架
files = dir('*.csv');
for k = 1:length(files)
data = readmatrix(files(k).name);
% 自动化分析...
end
时间对齐：多通道信号同步之道
多个传感器可能存在延迟，需要用互相关估计并补偿：
[corr, lags] = xcorr(signal1, signal2);
[~, idx] = max(abs(corr));
delay_samples = lags(idx);
if delay_samples > 0
aligned_sig2 = [signal2(delay_samples+1:end); zeros(delay_samples,1)];
else
aligned_sig2 = [zeros(-delay_samples,1); signal2(1:end+delay_samples)];
end
去趋势与去均值：清除低频干扰
% 去DC
x_clean = x - mean(x);
% 去线性趋势
x_clean = detrend(x);
% 去多项式趋势（如温漂）
p = polyfit(t, x, 2);
trend = polyval(p, t);
x_clean = x - trend;
| 处理方式 | 对频谱的影响 |
|---|---|
| 不处理 | 0Hz处巨大峰值，掩盖低频特征 |
| 去均值 | 抑制DC峰，暴露真实波动 |
| 去趋势 | 消除缓慢漂移，改善基线 |
异常值检测与平滑
IQR法识别离群点
Q1 = prctile(x, 25);
Q3 = prctile(x, 75);
IQR = Q3 - Q1;
lower_bound = Q1 - 1.5*IQR;
upper_bound = Q3 + 1.5*IQR;
outliers = (x < lower_bound) | (x > upper_bound);
x_clean = fillmissing(x, 'linear');
x_clean(outliers) = NaN;
平滑滤波选择
% 移动平均（去白噪声）
smoothed_ma = movmean(x, 5);
% Savitzky-Golay（保峰性能好）
smoothed_sg = sgolayfilt(x, 2, 11);
👉 推荐：振动/光谱类信号优先用sgolayfilt，避免破坏尖锐特征。
🪟 窗函数实战：抑制频谱泄漏的终极武器
当你截断一段非整周期信号时，等于乘上了矩形窗，导致能量扩散——这就是“频谱泄漏”。
泄漏原理简析
设信号 $ x(t) = \cos(2\pi f_0 t) $，若观测时间 $ T $ 不是周期整数倍，则其频谱不再是δ函数，而是sinc函数展宽：
$$
X(f) \propto \text{sinc}\left((f - f_0)T\right)
$$
后果就是主瓣展宽、旁瓣起伏，严重时会淹没弱信号。
解决方案？加窗！
graph TD
A[原始信号] --> B{是否整周期?}
B -- 是 --> C[窄谱线, 无泄漏]
B -- 否 --> D[能量扩散, 主瓣展宽]
D --> E[引入窗函数抑制旁瓣]
常见窗函数对比
| 窗函数 | 主瓣宽度（相对） | 旁瓣衰减（dB） | 适用场景 |
|---|---|---|---|
| 矩形窗 | 1.0 | -13 | 整周期、高分辨率 |
| 汉宁窗 | 1.5 | -31 | 通用平衡型 |
| 汉明窗 | 1.36 | -41 | 抑制近旁瓣 |
| 凯撒窗 | 可调（β参数） | 可达-58 | 自适应设计 |
win_hamming = hamming(N);
x_windowed = x .* win_hamming;
⚠️ 注意：加窗会降低信噪比，并使主瓣变宽。这是一种 分辨率与泄漏的权衡 。
能量补偿：别忘了恢复真实幅值
加窗后信号总能量减少，需进行幅度补偿：
coherent_gain = sum(win)/N;
compensation_factor = 1 / coherent_gain;
X_comp = fft(x_win) * compensation_factor / N;
封装成函数更方便复用：
function [f, P1] = compute_spectrum_with_window(x, fs, window_type)
N = length(x);
switch window_type
case 'rectangular'
win = ones(N,1);
case 'hann'
win = hann(N);
case 'hamming'
win = hamming(N);
case 'kaiser'
win = kaiser(N, 8);
end
x_win = x .* win;
cg = sum(win)/N;
X = fft(x_win) / (N * cg);
P2 = abs(X);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(N/2))/N;
end
调用起来超轻松：
[f, P1] = compute_spectrum_with_window(signal, fs, 'hamming');
plot(f, P1); grid on;
📊 归一化与能量守恒：让频谱更有物理意义
单边谱构造
X = fft(x);
P2 = abs(X)/N;
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1); % 中间点加倍
f = fs*(0:(N/2))/N;
📝 直流和Nyquist频率不加倍！
功率谱密度（PSD）
单位：$ V^2/Hz $
PSD = (abs(X).^2) / (fs * N);
PSD(2:end-1) = 2*PSD(2:end-1);
或者直接用MATLAB推荐的 pwelch ：
[PSD_welch, f_welch] = pwelch(x, hamming(256), 128, N, fs);
更适合噪声信号分析，方差更小 ✅。
Parseval定理验证：能量守恒检查
时域能量应等于频域能量：
$$
\sum |x[n]|^2 = \frac{1}{N} \sum |X[k]|^2
$$
验证代码：
energy_time = sum(abs(x).^2);
energy_freq = sum(abs(X).^2) / N;
fprintf('误差: %.2e\n', abs(energy_time - energy_freq)); % 应接近0
测试结果表明，几乎所有信号都能满足该定理（误差在浮点精度范围内），证明我们的流程是可靠的。
🖼 频谱可视化：画出专业级图表
构建频率轴（等效Python的fftfreq）
Fs = 1000; NFFT = 1024;
f = (0:NFFT/2)*Fs/NFFT; % 单边频率轴
支持补零插值，但不会提升真实分辨率。
绘制清晰频谱图
X = fft(x, NFFT);
X_mag = abs(X(1:NFFT/2+1))/N * 2;
X_psd = (abs(X(1:NFFT/2+1)).^2)/(N*Fs);
figure;
subplot(2,1,1);
plot(f, X_mag, 'b-', 'LineWidth', 1.2);
xlabel('频率 (Hz)'); ylabel('幅度 (V)');
title('单边幅度谱'); grid on;
subplot(2,1,2);
plot(f, 10*log10(X_psd), 'r-', 'LineWidth', 1.2);
xlabel('频率 (Hz)'); ylabel('功率密度 (dB/Hz)');
title('功率谱密度 (PSD)'); grid on;
时频联合分析：spectrogram上场！
对于非平稳信号（如启动过程、语音、地震波），要用短时傅里叶变换（STFT）：
window = hamming(256);
noverlap = 128;
nfft = 512;
[~, F, T, P] = spectrogram(x, window, noverlap, nfft, Fs);
figure;
imagesc(T, F, 10*log10(P));
axis xy; colorbar;
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('时频谱图（STFT）');
完整流程如下：
graph TD
A[原始时域信号] --> B{是否平稳?}
B -- 是 --> C[执行FFT]
B -- 否 --> D[分段加窗做STFT]
C --> E[构造频率轴]
D --> F[生成时频矩阵]
E --> G[绘制幅度/功率谱]
F --> H[绘制spectrogram]
G --> I[标注峰值频率]
H --> J[识别频率演化趋势]
I --> K[输出频域特征报告]
J --> K
这套流程已在工业监测、故障诊断、音频分析等多个领域得到广泛应用。
✅ 总结：构建你的频谱分析思维框架
经过这一番深入剖析，你应该已经明白：
- FFT ≠ 分辨率提升器 ：真正的分辨率取决于观测时间；
- 补零 ≠ 增加信息 ：只是频域插值，用于改善显示效果；
- 加窗是一把双刃剑 ：抑制泄漏的同时也会展宽主瓣；
- 预处理决定成败 ：去趋势、去噪、对齐缺一不可；
- 归一化才有物理意义 ：不然看到的只是“数字游戏”。
最终的目标是在时间、内存、精度之间找到最佳平衡点。而这，正是工程的艺术所在 🎨。
下次当你面对一段复杂的信号时，不妨问自己这几个问题：
1. 我的采集时间够长吗？
2. 是否存在趋势或异常值？
3. 需不需要加窗？选哪种窗？
4. 频谱是否做了正确归一化？
只要把这些细节都把控到位，你也能做出让人眼前一亮的专业级频谱分析报告！📈✨
简介：快速傅里叶变换（FFT）是信号处理中将时域信号转换为频域的关键技术，广泛应用于频谱分析和功率谱分析。本文以MATLAB为工具，深入讲解如何利用 fft 函数进行信号的频域变换，涵盖从数据读取、加窗处理到频谱与功率谱计算的完整流程。通过 fft.m 脚本实例，介绍采样率、频率分辨率、窗口函数应用及频率轴校正等核心概念，并结合 plot 、 specgram 等函数实现可视化，帮助用户掌握噪声分析、滤波设计和系统识别中的关键技术。本项目适用于信号处理初学者及工程实践者，具备较强的可操作性与实用价值。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
