简介：这份MATLAB代码包聚焦快速傅里叶变换（FFT）在信号处理中的应用，并通过两个独立脚本对比傅里叶变换与小波变换在信号消噪中的实际效果，适合MATLAB初学者以及需要处理非平稳信号的工程人员。代码包内共2个M文件，分别用于对比FFT与WT消噪流程和供用户调整参数实践，压缩包仅2KB，结构精简易读。已有1887人学习，是信号处理领域中较为常见的参考实现。通过代码可学习MATLAB中fft函数的使用、fftshift中心对齐、abs与平方运算求幅度谱，理解频域图中横坐标为角频率、纵坐标为幅值的含义。同时掌握wavedec小波分解、wthresh阈值处理与waverec重构的完整消噪流程，便于直观比较两种方法在瞬态噪声下的性能差异，为实际工程中选择去噪方案提供参考。
1. 从一段含噪信号说起：FFT到底解决了什么，解决不了什么
拿到一段振动、电流或声学数据，绝大多数人的第一步是plot看时域波形，第二步就是fft看频谱。这个流程本身没问题，但很多人在第二步就停了：画完幅值谱，确认了主频，然后就没有然后了。真正让人头疼的问题往往在后面——信号带噪声怎么消、非平稳信号怎么分析、FFT 给出的幅值到底怎么换算才对。标题里的matlab代码_fft_这一组练习代码，正好覆盖了这两条线索：Untitle_practice.m是 FFT 频谱分析与幅值谱绘制的基本功练习，Untitled_compareWTwithFT.m则把傅里叶变换和小波变换放在同一个消噪任务下做对比。这篇博文就按这两条线展开，从fft函数的真实参数行为讲起，到 FT 与 WT 消噪各自适用的信号类型，最后给出一套可以自己改参数的对比实验模板。适合已经能跑通基础 MATLAB、但对频率轴换算和消噪方法选型还不太确定的读者。
2. FFT 在 MATLAB 中的实现：频率轴、幅值谱与功率谱的换算
2.1fft函数的行为与点数对齐：先搞清楚长度和采样率
MATLAB 里fft(X)和fft(X, N)行为不同，这一点在实际处理时经常被忽略。fft(X)返回长度与X相同的结果；fft(X, N)会在N大于信号长度时补零，在N小于长度时截断信号。补零不会提高频率分辨率——分辨率由真实信号时长决定，补零只是让频谱曲线更「光滑」。这一点在对比不同长度信号时尤其重要。
fs = 1000; % 采样率 1000 Hz t = (0:999) / fs; % 1 秒，1000 个点 x = 0.8*sin(2*pi*50*t) + 0.4*sin(2*pi*120*t); N = length(x); % FFT 点数 X = fft(x, N); % N 点 FFT f = (0:N-1) * fs / N; % 频率轴，单位 Hz plot(f, abs(X)); xlabel('Frequency (Hz)'); ylabel('Magnitude');这段代码的关键在于频率轴的构造：(0:N-1) * fs / N把 FFT 输出的第k个点映射到物理频率k * fs / NHz。如果只画abs(X)不管横轴，你会看到两个在 50 Hz 和 120 Hz 处的对称尖峰，但横轴的数字是完全错的。离散傅里叶变换输出的频谱关于N/2对称，N/2对应奈奎斯特频率fs/2，高于这个频率的部分在物理上是负频率的镜像。fft函数的输出中包含直流分量和正负频率的完整信息，直接取绝对值看到的是双边谱。
2.2 单边幅值谱的换算：为什么要乘 2
工程分析时我们通常只关心正频率，所以要把双边谱折叠成单边谱。常见做法是对除直流和奈奎斯特点外的所有正频率幅值乘 2，再统一除以N做归一化。
X = fft(x, N); A = abs(X) / N; % 归一化幅值 A_single = A(1:N/2+1); % 取正频率部分 A_single(2:end-1) = 2 * A_single(2:end-1); % 负频率能量并回正频率 f_single = (0:N/2) * fs / N; plot(f_single, A_single); xlabel('Frequency (Hz)'); ylabel('Amplitude');2 * A_single(2:end-1)这一行是很多初学者最容易漏掉的。FFT 是线性变换，原始信号的能量被平均分到了正负频率两个镜像上，单边谱要把负半轴的能量加回来，所以除直流点A_single(1)和奈奎斯特点A_single(end)之外都要乘 2。乘完之后，50 Hz 分量的幅值应该约为 0.8，120 Hz 分量约为 0.4，直接对应原始信号的幅值。如果看到幅值只有预期的一半，不用怀疑代码逻辑，先查这一行有没有写对。
| 参数 | 含义 | 计算方式 |
|---|---|---|
| 频率分辨率 Δf | 相邻两个频点的间隔 | fs / N，只由真实信号时长决定 |
| 最大分析频率 | 奈奎斯特频率fs/2 | FFT 能表示的最高物理频率 |
| 单边谱乘 2 | 把负频率能量并回正频率 | A_single(2:end-1) = 2 * A_single(2:end-1) |
| 角频率换算 | 从 Hz 转 rad/s | w = 2 * pi * f |
表格里有一项值得单独说明：如果你看到的参考书或论文里写「横坐标为角频率，纵坐标为幅值」，那是在用ω = 2πf画图，纵轴幅值不变，横轴数值全乘2π。比如 50 Hz 对应约 314 rad/s。实际做信号分析时用 Hz 更直观，做理论推导或阶次分析时用角频率更多，二者只是坐标缩放，不影响频谱峰值的位置关系。
2.3fftshift中心对齐与功率谱
fftshift用于把零频分量移到频谱中央，便于观察直流成分和频谱的对称结构。Hilbert 变换、调制解调分析里经常要用到这种视图。
X_shift = fftshift(fft(x)); f_shift = (-N/2:N/2-1) * fs / N; plot(f_shift, abs(X_shift));fftshift之后横轴要以(-N/2:N/2-1)重新构造，前面(0:N-1)的坐标轴已经不再适用。另一个容易被忽略的点是：fftshift只改变排列顺序，不改变数值，所以它对幅值谱和相位谱的处理方式完全相同。
功率谱的获取也常被一笔带过。直接对X取模平方得到的是信号的周期图，若需要功率谱密度估计，要除以fs * N；若只需要各频率成分的相对能量强弱，用abs(X).^2 / N^2即可，单位是信号幅值的平方。功率谱的峰值位置和幅值谱一致，但高低频分量的相对差距会被平方放大，这在消噪场景中判断「哪些频段值得保留」时更直观。
提示：
fft输出的第一个点是直流分量，即信号均值乘N。做频谱分析前先看时域信号是否去均值，趋势项不去掉会把低频段整体抬高，小波消噪时也会被误判成有效成分。
3. FT 消噪的局限与小波分解消噪的实现
3.1 为什么纯 FFT 路线不适合消噪
很多人在了解了 FFT 之后，自然会想「把噪声频段的系数置零，再ifft回来不就能消噪了」。这个思路在教科书上成立，实际用起来却处处受限。白噪声在频域里是平坦的，和信号的频谱在整个频带上重叠，简单地把高频系数置零，相当于让一个截止频率以上的所有成分全部丢失，信号里的陡峭边沿和瞬态脉冲也会被同时抹平。更麻烦的是，直接在频域做硬截断，重构时会在断点处产生 Gibbs 现象，表现在时域就是信号两端出现不衰减的振铃，这个振铃并不是真实信号，而是截断本身带来的伪迹。
标准的频域消噪路径应该是设计一个带通滤波器，而不是手动把 FFT 系数清零。比如用butter设计巴特沃斯滤波器，再用filtfilt做零相位滤波：
[b, a] = butter(4, [45 125] / (fs/2), 'bandpass'); x_filt = filtfilt(b, a, x_noisy);这里[45 125] / (fs/2)是把通带频率归一化到奈奎斯特频率，butter的参数 4 表示滤波器阶数，阶数越高过渡带越窄，但相位失真也更严重。filtfilt做了双向滤波，零相位、没有群延迟，代价是计算量翻倍。这套路线的核心问题不在实现，而在参数选择：通带边界必须由你先从频谱上判断出来——一旦信号里混着多根谱线，或者噪声根本不是白噪声，你根本不知道该保留哪个频段。
3.2wavedec分解与噪声标准差估计
小波消噪走的是另一条路：把信号分解成不同尺度（对应不同频段）的细节系数和近似系数，再对细节系数做阈值处理。wavedec返回的[c, l]结构里，c是拼接在一起的系数向量，l记录每一段长度。分解结构的顺序是[近似系数, 最高层细节, ..., 最底层细节]，所以第一层（最高频）细节系数位于c的末尾l(1)个点。
% 生成含噪信号 rng(0); x0 = 0.8*sin(2*pi*50*t) + 0.4*sin(2*pi*120*t); x_noisy = x0 + 0.15*randn(size(x0)); wname = 'db4'; level = 5; [c, l] = wavedec(x_noisy, level, wname); % 用最高频细节系数估计噪声标准差 d1 = c(end-l(1)+1:end); % 第一层细节系数 sigma = median(abs(d1)) / 0.6745; % 鲁棒标准差估计 thr = sigma * sqrt(2 * log(length(x_noisy))); % 通用阈值 c_soft = wthresh(c, 's', thr); % 软阈值处理 x_wt = waverec(c_soft, l, wname); % 重构这段代码里有两个值得展开的参数。第一，median(abs(d1)) / 0.6745是利用标准正态分布的性质做鲁棒标准差估计，0.6745 是标准正态分布 75% 分位数。中位数比均值抗离群点，所以即使第一层细节里混有少量真实信号的瞬态成分，这个估计也不会被带偏。第二，thr = sigma * sqrt(2 * log(N))是 Donoho 提出的通用阈值，来源于极值理论，意思是「白噪声在 N 个采样点里产生的最大幅值大概率不会超过这个界限，超过的部分才被认为是有效信号」。
3.3 软阈值与硬阈值的选择依据
wthresh的参数's'和'h'分别对应软阈值和硬阈值，行为差异对重构结果影响很大。
| 阈值方式 | 系数的处理规则 | 优点 | 缺点 |
|---|---|---|---|
'h'硬阈值 | 绝对值小于阈值的置零，其余保留原值 | 峰值幅度保留完整，适合信号本身有尖峰的场景 | 阈值处不连续，重构信号易出现局部抖动 |
's'软阈值 | 绝对值小于阈值的置零，其余向零收缩thr | 重构波形连续光滑，噪声抑制更彻底 | 所有保留的系数都被压缩，幅值偏小 |
| 多级阈值 | thr改为向量，逐层传入 | 每层噪声水平不同，可分别处理 | 参数数量多，需要观察每层系数分布 |
硬阈值重构的信号在突变位置更接近原始信号，但会在阈值边界附近出现不连续的小锯齿；软阈值整体更平滑，代价是真实信号的高频分量也被均匀压缩了一点。比较稳妥的做法是先看各层细节系数的分布：如果某一层系数的直方图有明显双峰（一个峰在零附近，另一个在远处），用硬阈值；如果只在零附近有一团密集的小值，用软阈值。
提示：
wthresh支持thr为向量，长度等于分解层数。更精细的做法是用wdencmp逐层传阈值，但前提是你对每层噪声水平有直观认识——先用wavedec分解一次，分别画出每层细节系数再决定。
4. 对比 FT 与 WT 消噪效果的脚本框架与实验设计
4.1 一套可复现的对比实验
Untitled_compareWTwithFT.m这类脚本的核心价值不在消噪算法本身，而在它把两种方法放到同一个评价体系下比较。对比实验最怕的是「FT 用了这个参数、WT 用了那个参数，最后结果不可比」。我一般把对比脚本设计成三个固定：固定同一段含噪信号、固定同一条评价链路、固定可复现的随机种子。
rng(1); fs = 1000; t = (0:999) / fs; x0 = 0.8*sin(2*pi*50*t) + 0.4*sin(2*pi*120*t); x0(500) = x0(500) + 2; % 在 0.5s 处加一个瞬态脉冲 x_noisy = x0 + 0.2*randn(size(x0)); % FT 路线：先看频谱再设计带通滤波 [b, a] = butter(4, [45 125] / (fs/2), 'bandpass'); x_ft = filtfilt(b, a, x_noisy); % WT 路线：小波分解 + 软阈值 wname = 'db4'; level = 5; [c, l] = wavedec(x_noisy, level, wname); d1 = c(end-l(1)+1:end); sigma = median(abs(d1)) / 0.6745; thr = sigma * sqrt(2 * log(length(x_noisy))); x_wt = waverec(wthresh(c, 's', thr), l, wname); % 评价：SNR（信噪比） snr_ft = 10*log10(sum(x0.^2) / sum((x0 - x_ft).^2)); snr_wt = 10*log10(sum(x0.^2) / sum((x0 - x_wt).^2)); fprintf('FT: %.2f dB, WT: %.2f dB\n', snr_ft, snr_wt);这个脚本里最关键的设计是在信号里加了一个x0(500) = x0(500) + 2的瞬态脉冲。这个脉冲在频域里展布在整个频谱上，幅值又小，FT 路线很难针对它单独处理；而小波变换的细节系数在脉冲位置会出现明显的局部极大值，阈值处理后这个点依然被保留。忽略瞬态成分只会得到一个「两种方法差不多」的平淡结论，加上这个脉冲，FT 和 WT 的优势区间立刻分化出来。
4.2 参数扫描：不同噪声强度下的对比结果
单看一组 SNR 数字还不够，噪声强度变化时两种方法的表现会交叉。一般会扫一组噪声标准差sigma_noise，看 SNR 提升幅度随噪声强度的变化趋势：
| 噪声标准差 σ | 输入 SNR（dB） | FT 输出（dB） | WT 输出（dB） | 结论 |
|---|---|---|---|---|
| 0.05 | 22.5 | 31.0 | 33.6 | 噪声小，两种方法差距不大 |
| 0.2 | 10.5 | 19.4 | 24.1 | WT 明显占优，FT 带通后残留大量带外噪声 |
| 0.5 | 2.1 | 8.6 | 13.8 | WT 优势进一步拉大，FT 开始丢信号细节 |
表格里的数值是说明趋势的示意值，具体结果会随随机种子浮动，但趋势是稳定的。噪声小时 FT 和 WT 都够用，因为信噪比高，带通滤波残留的噪声影响有限；噪声增大后，白噪声频谱变高，固定通带的 FT 路线无能为力，而 WT 的阈值是依据当前信号噪声水平自适应算出来的，所以优势越来越大。这也解释了为什么实践里没人用纯 FFT 做消噪——固定频带对非平稳噪声完全不设防。
4.3 SNR 之外还要看什么
只比较 SNR 也会被误导。SNR 把整个时间段的所有误差平均成一个标量，脉冲位置有没有保住、端点有没有振铃、相频特性有没有畸变，这些局部现象在 SNR 里体现不出来。实际操作中至少再加两个观察维度：时域残差图和分段时间的误差分布。
第二个维度是运算形态，这直接关系到你能不能把这个脚本用到更大的数据上。wavedec对 1000 点数据是一瞬间的事，但当信号长度到百万量级、分解层数到 8 层以上时，循环处理每层系数才会有可感知的时间开销。我在对比脚本里一般会把wavedec拆成逐层detcoef取出来看每一层的能量占比，判断哪些层该保留、哪些层该全部置零，这比整段套用一个固定阈值更贴合实际信号的频带分布。
第三个维度是方法的时间复杂度差别。FFT 的复杂度是O(N log N)，小波分解同样是O(N log N)，两者在纯粹计算量上没有代差。真正拉开差距的是参数获取成本：FFT 路线需要你先做频谱分析、手工定通带，而小波阈值可以通过噪声统计自动算出。对自动化的批处理任务来说，WT 路线明显更合适。
提示：对比实验一定要用同一份含噪信号跑完所有方法，不能每次重新生成噪声。否则测出来的 SNR 差 0.5 dB，你无法判断是方法差异还是随机种子差异。
5. 边界条件、分解层数设定与实测排错技巧
5.1 端点效应与延拓方式的影响
小波消噪重构出的信号，两端往往比中间差。这是因为小波变换在信号边界处没有足够的数据计算卷积，默认延拓方式和真实信号的趋势不匹配，重构时就会在边界区留下振铃。wavedec默认的延拓模式是周期延拓'per'，对两端不连续的信号效果较差；改成对称延拓'sym'通常能压低边界伪迹。
% 对比延拓方式对端点重构误差的影响 [c_sym, l_sym] = wavedec(x_noisy, 5, 'db4', 'sym'); x_sym = waverec(c_sym, l_sym, 'db4', 'sym'); err_sym = norm(x_sym - x0) / norm(x0); [c_per, l_per] = wavedec(x_noisy, 5, 'db4', 'per'); x_per = waverec(c_per, l_per, 'db4', 'per'); err_per = norm(x_per - x0) / norm(x0);如果err_sym明显小于err_per，说明你的信号两端不连续，对称延拓更适合。但对称延拓也有代价：它假设边界上信号是镜像对称的，如果信号本身是正弦类周期信号，周期延拓误差更小。延拓方式没有绝对好坏，要用重构误差说话。
5.2 分解层数的上限与阈值偏移的排查
分解层数不是越大越好。每一层将频率范围减半，分解过深时最低层的近似系数已经几乎不包含有用信号，只剩整个信号的趋势项。经验法则是：目标信号的最低频率f_min与层数level满足fs / 2^(level+1)接近f_min。对 1000 Hz 采样、希望保留 7 Hz 以上成分的信号，1000 / 2^6 ≈ 15.6 Hz、1000 / 2^7 ≈ 7.8 Hz，取 6 层比较合理，7 层会把 7 Hz 以下的趋势也收进近似系数里。
阈值偏移是最常见的排错对象。如果阈值thr算得偏大，重构信号的幅值整体偏低，原因是软阈值把所有细节系数都压缩了一个thr量，用sum(x_wt.^2) / sum(x0.^2)算能量比就能看出来；如果阈值偏小，重构信号里还有明显的毛刺，看残差信号x_noisy - x_wt的频谱就能确认它和白噪声频谱形状是否一致。残差频谱如果还有明显谱峰，说明阈值把真实信号成分一起抑制了——这时应该降低层数或改用硬阈值，而不是调整阈值本身。
5.3 从 MATLAB 到嵌入式 FFT 的换算注意点
把这套流程移植到 STM32F4 这类 MCU 上做实时频谱分析时，最常踩的坑是定点数与浮点数的换算。MATLAB 里fft直接返回浮点复数，而嵌入式端常用 CMSIS-DSP 的arm_cfft_f32或 FFT IP 核，输入输出做了一次定点缩放，幅值谱的数值比 MATLAB 结果多一个固定比例系数，直接用 MATLAB 的阈值去对应嵌入式结果必然失效。一般在嵌入式端不做归一化，只比较各频点幅值的相对大小，阈值靠实测标定而不是理论计算。MATLAB 端算出的fs/N频率分辨率和fs/2奈奎斯特频率在两个平台上完全一致，这两项可以放心沿用到固件参数里。
本文还有配套的精品资源，点击获取
