企业官网建设流程全解析

做信号处理的，谁没跟频谱较过劲呢。现场采集回来的数据，传感器噪声、设备振动、环境干扰叠在一起，一条扭来扭去的波形里藏着好几个频率成分，想把它们挨个儿揪出来，功率谱估计是第一道绕不开的门槛。这篇文章要聊的是最经典、也最容易被低估的两种方法——BT法（Blackman-Tukey法）和周期图法。它们一个从自相关函数出发，一个直接对数据做变换，表面上走的是两条路，背后却共享同一个数学根基：维纳-辛钦定理。我会把自相关估计里那个“k除不除”的细节、窗函数到底怎么选、MATLAB代码怎么写，以及为什么你的谱图总是抖个不停，逐条拆开讲。适合刚开始接触随机信号分析的同学，也适合做故障诊断、振动分析、音频特征提取的工程师，哪怕你只想弄明白手里一段数据的频率构成，这篇文章也能给你一个可以直接照着跑的起点。

1. 谱估计到底在估什么：从相关函数到功率谱

1.1 为什么不能只靠一条FFT

很多刚上手的人第一反应是：对信号做一次FFT，再把模值平方，不就是功率谱吗？听起来很顺，但这里有一个隐蔽的坑。FFT处理的是确定性信号或者有限能量信号，而随机信号是随机过程的一次实现，你拿到的那一段数据只是所有可能样本里的一条路径。直接对它做FFT，得到的是这条样本路径的频谱，它的相位是随机抖动的，幅度也带有明显的随机起伏，严格说并不具备统计意义上的“频谱”含义。

平稳随机信号属于功率信号而不是能量信号，它的总能量理论上发散，但平均功率是有限且稳定的。所以工程上要估计的是“功率随频率的分布”，也就是功率谱密度（PSD）。维纳-辛钦定理告诉我们，平稳随机信号的自相关函数和功率谱密度是一对傅里叶变换：

S(f) = Σ_{k=-∞}^{∞} r(k) · e^{-j2πfk}

这里 r(k) 是自相关函数。BT法的思路就非常清晰了：先用观测数据估计自相关 r(k)，再对它做离散时间傅里叶变换，得到功率谱估计。换句话说，BT法是“先相关、后谱”的路线，它把随机信号的处理转换成了一个确定性的数学运算。

周期图法走的则是另一条路。它直接对观测数据做FFT，然后取模平方，再除以数据长度：

I(f) = (1/N) · |Σ_{n=0}^{N-1} x(n) · e^{-j2πfn}|²

从表达式上看没有出现自相关函数，可实际上它隐含了数据中所有自相关的信息。FFT模平方展开以后正好是自相关的傅里叶变换形式，只是这里的自相关是由这一段数据直接堆出来的，后面没有再给你加窗和截断的自由度。

1.2 两种方法的血缘关系

BT法和周期图法并不是两个毫无关系的算法。如果令BT法里的最大滞后点数 M = N-1，并且自相关窗选矩形窗，BT法估计出来的结果和周期图法几乎一样。两者的差别主要体现在可调自由度上。

实际工程里，周期图法胜在简单直接，一条FFT链路就能跑起来，适合现场快速看全貌；BT法则多了一个 M 作为旋钮，你可以主动控制谱的平滑程度，这也是它在现代谱估计方法出现之前长期占据工具书页面的原因。

2. 数学细节：自相关怎么估、窗函数怎么选

2.1 有偏估计与无偏估计：括号里那个除法

写自相关估计式的时候，有一个细节特别容易被忽略，就是分母到底是 N 还是 N-k。对付一段长度为 N 的数据，滞后 k 的自相关估计可以写成两种形式：

有偏估计：

r_hat(k) = (1/N) · Σ_{n=1}^{N-k} x(n) · x(n+k)

无偏估计：

r_hat(k) = (1/(N-k)) · Σ_{n=1}^{N-k} x(n) · x(n+k)

两个式子只差一个分母。无偏估计的数学期望严格等于真实自相关 r(k)，听起来更“正确”，但工程里几乎总是用有偏估计。原因在于当 k 接近 N 的时候，比如 k = N-1，参与计算的有效样本只有一对数据，系统偏差很大，估计值就会在那个点上剧烈跳动。这个高方差的自相关值，经过傅里叶变换以后，会在功率谱底部激起一堆细碎的假峰。

有偏估计虽然对长滞后部分有系统偏差，但它的好处是把所有滞后都用同一个 N 归一，尾部被压得比较光滑，最后计算出来的功率谱不容易出现负值。严格说，有偏自相关矩阵在Toeplitz意义下保持半正定，这保证了谱估计的非负性，这在物理上更合理。你自己也可以做一个简单实验：生成一段白噪声，分别用两种估计画 r(k)，会看到无偏估计在 k 稍大的位置尾巴像刺猬一样炸开，而有偏估计则干净得多。

2.2 窗函数不是玄学：主瓣与旁瓣的交换

BT法里，估计出频率范围为 -M 到 M 的自相关序列以后，还需要乘一个窗函数再送进FFT。很多初学者不理解这一步在干嘛，觉得窗口只是为了让曲线好看。其实窗函数解决的是谱估计里的核心矛盾：泄漏与分辨率的权衡。

任意有限长窗的频谱都由主瓣和旁瓣构成。主瓣越窄，两个相邻频率越容易被分开，分辨率越高；旁瓣越低，则强频率成分对周围频段的泄漏越少，动态范围越好。但这两个目标不能同时满足。以常用的窗函数为例：

从表格能看到，矩形窗的主瓣最窄，分辨率最好，但旁瓣只衰减 -13 dB，一个强正弦的泄漏足以淹没旁边的小信号。汉宁窗、汉明窗牺牲了主瓣宽度，换来更干净的旁瓣，所以做谱估计时，除非你明确知道自己要极致的分辨率，否则默认选择汉宁窗是稳妥的。

BT法里的窗函数作用对象是自相关估计值 r(k)，它在频域里会对真实功率谱做卷积平滑。这个平滑会导致谱峰变宽，但也杀掉了 r(k) 尾部的高方差成分。这就是为什么BT法可以通过调整 M 和窗函数，在谱估计的“偏差”和“方差”之间做交易。

2.3 关于M：BT法独有的旋钮

M 是BT法最关键的自由度，它的含义是你决定把自相关估计保留到第几个滞后。M 取小，尾部噪声被截掉，谱线平滑，但主瓣变宽，相近频率容易糊在一起；M 取大，自相关保留得长，分辨率提升，但尾部估计的方差又把谱线弄得毛刺丛生。

我建议的经验范围是 M = N/8 到 N/2 之间。具体取多少，要看你想分辨的最小频率间隔。如果已知两个频率间隔为 Δf（Hz），采样率为 fs，那么粗略估算需要 M ≥ fs / Δf。比如采样率 1000 Hz，想分辨 5 Hz 间隔的频率，M 至少要 200。先按这个公式算一个初值，再往小里调，观察谱线是否变得过于平滑；或往大里调，看是否出现许多细碎波动。记住一点：BT法里 M 增大带来的“假细节”，远没有现代谱估计里那种统计特性，别迷信大 M。

3. MATLAB代码实现：从零开始写BT法和周期图法

3.1 测试信号：老规矩，正弦叠噪声

为了检验算法，我们生成一个谱结构已知的信号：两个频率分量加高斯白噪声。之所以用这个信号，是因为结果一目了然——两个峰该出现在哪里、峰形该多宽、噪声本底该多高，心里都有数。代码里的一次性数据可以先用固定随机种子，方便对照现象。

clear; close all; fs = 1000; % 采样率 1000 Hz N = 1024; % 数据点数 t = (0:N-1) / fs; % 时间向量 f1 = 50; % 第一个正弦 50 Hz f2 = 120; % 第二个正弦 120 Hz xn = 1.5*sin(2*pi*f1*t) + 0.8*sin(2*pi*f2*t); xn = xn + 0.8*randn(1,N); % 加性高斯白噪声 figure; plot(t, xn); grid on; xlabel('时间/s'); ylabel('幅值'); title('测试信号时域波形');

注意 randn 每次执行结果不同，所以下面的谱图每次也会有细节差异，但两个主峰的位置和高度应该相对稳定。如果想严格复现，可以在开头加一行 rng(1)。

3.2 BT法完整实现

下面给出一个自包含的MATLAB函数。我没有直接调用信号处理工具箱里的 xcorr，而是用循环实现自相关估计，这样没有工具箱的读者也能跑。顺带提一句，xcorr 默认会把结果补成 2N-1 长度，切出中间那一段其实也是等效的，只是新手容易在这里把长度搞错。

function [f, Pbt] = bt_psd(x, fs, M, winType) % BT法功率谱估计 % x : 输入信号（行向量） % fs : 采样率（Hz） % M : 自相关保留的最大滞后点数 % winType: 'rect' 、'hann' 或 'hamming' x = x(:)'; % 统一为行向量 x = x - mean(x); % 去均值，避免直流泄漏 N = length(x); % 1. 有偏自相关估计 r = zeros(1, M+1); for k = 0:M r(k+1) = sum( x(1:N-k) .* x(k+1:N) ) / N; end % 2. 自相关窗 switch lower(winType) case 'hann' w = 0.5*(1 - cos(2*pi*(0:M)/M)); case 'hamming' w = 0.54 - 0.46*cos(2*pi*(0:M)/M); otherwise w = ones(1, M+1); end % 3. 加窗并做谱估计 rw = r .* w; rw = rw / mean(w); % 窗能量归一化，保持相对幅度可比 nfft = max(1024, 2^nextpow2(2*M)); Pbt = abs(fft(rw, nfft)).^2; % 4. 整理为单边谱 Pbt = Pbt(1:nfft/2+1); f = (0:nfft/2) * fs / nfft; Pbt(2:end-1) = Pbt(2:end-1) * 2; % 单边谱补偿功率 end

主脚本里调用：

[f_bt, Pbt] = bt_psd(xn, fs, 256, 'hann'); figure; plot(f_bt, 10*log10(Pbt + eps), 'LineWidth', 1.2); grid on; xlim([0 300]); ylim([-40 20]); xlabel('频率/Hz'); ylabel('功率谱密度/dB'); title('BT法功率谱估计（M=256, Hanning窗）');

代码里除以 mean(w) 这步，是为了保证加窗后自相关的总能量不因为窗口系数而被整体压低。如果你只是看峰值相对高度，这步可以省；但如果后面还要做信号功率积分对比，就必须保留。

3.3 周期图法实现与频率轴校正

周期图法写起来比BT法短，但有两个细节容易出错：一是单边谱补偿，二是频率轴换算。下面的函数把两个问题都处理好了。

function [f, Pper] = periodogram_my(x, fs, winType) % 周期图法功率谱估计 % x: 输入信号（行向量） x = x(:)'; x = x - mean(x); N = length(x); switch lower(winType) case 'hann' w = 0.5*(1 - cos(2*pi*(0:N-1)/(N-1))); otherwise w = ones(1, N); end w = w / sqrt(sum(w.^2)) * sqrt(N); % 保持加窗后信号平均功率一致 xw = x .* w; X = fft(xw); Pper = abs(X).^2 / N; % 双边功率谱估计 % 取单边谱 Pper = Pper(1:N/2+1); f = (0:N/2) * fs / N; Pper(2:end-1) = Pper(2:end-1) * 2; end

主脚本调用：

[f_per, Pper] = periodogram_my(xn, fs, 'hann'); figure; plot(f_per, 10*log10(Pper + eps), 'LineWidth', 1.2); grid on; xlim([0 300]); ylim([-40 20]); xlabel('频率/Hz'); ylabel('功率谱密度/dB'); title('周期图法功率谱估计（Hanning窗）');

这里单独说明一下单边谱的功率折算。FFT结果是对称的，正频率和负频率各分了一半能量，所以要画单边谱就得把除直流以外的正频点乘以2。严格说，0频率和奈奎斯特频率这两个点不乘2，代码里用 Pper(2:end-1) 正好避开了这两个位置。很多错误的谱图，峰位置对但幅度比理论值低了约一倍，多半是在这一步漏乘了。

3.4 两张图怎么摆在一起看

把两条曲线画在同一张图上，可以直观比较两种法的差异。我建议用同一组数据、同样的窗类型，只让算法路径不同，这样差别才纯粹来自方法本身。

figure; plot(f_bt, 10*log10(Pbt + eps), 'b', 'LineWidth', 1.2); hold on; plot(f_per, 10*log10(Pper + eps), 'g', 'LineWidth', 1.2); legend('BT法', '周期图法', 'Location', 'northeast'); grid on; xlim([0 300]); xlabel('频率/Hz'); ylabel('功率谱密度/dB'); title('BT法与周期图法对比');

跑完以后，两个峰应该都在 50 Hz 和 120 Hz 附近出现。BT法因为多做了一个自相关加窗平滑，峰形通常更圆润一点；周期图法单条曲线的毛刺会多不少，尤其在两峰之间的噪声区域，起伏更明显。这不是代码写错了，而是周期图法方差大的固有表现，下一章专门讲这个。

4. 分辨率与方差：为什么结果老在抖

4.1 分辨率的下限：瑞利极限与补零误区

很多资料里会看到“频率分辨率约等于 fs/N”，这是有限长数据谱估计的理论下限。两个正弦频率相差不到 fs/N 时，无论用什么谱估计方法，谱峰都会糊成一个包，这一点与算法优劣无关，是数据长度决定的。

有一个常见误区我必须强调：给FFT补零并不能提高真实分辨率。补零只是把离散谱的横坐标插值得更密，让峰形看起来更顺滑，但两个靠得很近的频率依旧无法分开。比如采样率 1000 Hz、数据长度 N = 256，频率分辨率大约是 3.9 Hz。如果你想分辨 80 Hz 和 82 Hz，理论上就很难做到；即使把FFT补到 4096 点，也只是把原来的包络画得更细，不会凭空多出第二个峰。想真正提高分辨率，只有延长观测时间 N 或采用现代参数化谱估计方法，这条路留给下一篇。

如果用BT法，由于自相关加窗，实际分辨率还要进一步受 M 制约。M 越小，分辨率越差。所以在工程里我不会把BT法和周期图法做“谁更高级”的比较，它们只是给了你不同位置的控制旋钮——周期图法的分辨率锁死在 N 上，BT法的分辨率由 M 决定，代价是方差特性也随之改变。

4.2 周期图法的方差问题：无偏但不一致

随机信号分析里有一句经典表述：周期图法是渐近无偏的，但不是一致估计。意思是当数据长度 N 增加时，周期图在每一个频率点的期望会逼近真实谱，但它的方差并不会收敛到 0，甚至可能保持在和谱本身同量级的水平。

生活化的类比是：一个人每天量一次血压，单次测量值总有波动；你今天量出的值是 120，明天可能是 128，后天是 115。你并不能因为某天量出 120，就认定血压精准等于 120。周期图法就是把“某一次测量”直接当结论，所以谱线会明显地上下跳动。

你可以自己做一个蒙特卡洛实验，把上一章的信号生成过程重复几十次，每次都画一条周期图叠加在一起。结果就是图上一团“意大利面”，谱线包围的区域非常宽，这代表单次周期图的置信度很低。BT法通过选择合适的 M，相当于对多次测量做了隐式平滑，方差明显下降，代价是谱峰变宽。这就是谱估计里最核心的偏差-方差权衡：你不可能在固定数据长度下同时拥有高分辨率和低方差，只能在两者之间找平衡。

4.3 M怎么选：从数据“摸”出来的经验

M 的选取并没有一个放之四海而皆准的公式，但有一些可以落地的经验法则：

• 宽带噪声为主、关心整体包络：M 选小一些，比如 N/8 到 N/4。这时谱线很平滑，底噪不会乱跳。
• 窄带信号、需分辨多个峰：M 先取 N/2 附近起步，再逐步减小，观察峰数量和峰位是否稳定。
• 想分辨 Δf Hz 的频率间隔：估算 M ≈ fs / Δf，这是窗主瓣宽度决定的底线，低于它会直接糊峰。
• 避免直接取 M = N-1：此时BT法退化成接近周期图法的形态，长滞后自相关的尾部方差又被引进来，失去了加窗平滑的意义。

我自己的习惯是写一个小循环，把 M 从 N/8 到 N/2 按 1.5 倍递增，画三条谱线放在同一张图里。哪个 M 下峰的数量不变化、底噪不过分抬升，就选那一个。这个方法粗暴但很有效，比我背任何经验公式都可靠。

5. 常见问题与排查技巧实录

5.1 谱峰为什么总在“抖动”

初次跑周期图法的朋友经常会来问我：为什么同样的代码，每次画出来谱线都不一样？因为信号里的 randn 噪声每次都在变，频谱估计本身就有随机性，这是正常的。你只需要先固定随机种子 rng(1)，就能复现同一张图，用来调试。如果你需要更稳定的谱估计结果，单靠增加 N 并不会改善方差，正确做法是把数据分段做平均，比如Welch方法，分段越多，平均后的曲线越平滑，代价是频率分辨率会下降。

另外，如果你画的是功率谱密度而不是累积功率谱，曲线本身就很“毛”，这不是bug。谱密度关注的是每个频率点附近的功率集中程度，单点值跳动大是统计特性决定的。

5.2 自相关尾部“毛刺”破坏谱图

BT法特有的一个现象：功率谱底部出现大量细碎假峰，尤其在高频段明显。问题多半出在自相关函数的尾部估计上。滞后 k 越大，参与乘加的样本对越少，估计方差越大，这些抖动经过FFT以后就变成了底部噪声。

排查方法是把自相关序列 r(0:M) 直接画出来看。如果尾部不是平滑衰减到 0，而是在零附近剧烈振荡，就说明 M 取得太大。把 M 减小一半再跑，假峰通常会明显减少。你也可以在 r(k) 上加一个巴特利特窗（三角窗），让自相关尾部自然衰减到 0，这个窗对抑制尾部方差很有效。

5.3 直流、趋势项与低频大峰

实际采集的数据很少是完全零均值的，传感器漂移、放大器温漂都会带来直流成分甚至缓慢变化的趋势项。这些低频成分在谱估计里会表现为一个很大的低频峰，足以把附近真正的信号峰淹没。处理办法很简单：先对信号去均值，再做一次去趋势处理。

% 去除直流 x = x - mean(x); % 去除线性趋势，假设x是行向量 n = 0:length(x)-1; p = polyfit(n, x, 1); x = x - polyval(p, n);

如果你的环境里有 detrend 函数，也可以直接用 x = detrend(x(:)')。这一步对周期图法尤其重要，因为FFT的周期性默认处理的是零均值信号，直流分量一旦泄漏，整个低频段的谱都会被抬高。

5.4 频率轴对不齐与坐标单位换算

代码报错里最常见的是画图时 f 和 P 长度不一致。写代码前先想清楚：单边谱的频率轴长度应该是 N/2+1 个点，而不是 N 个点。用 f = (0:N/2) * fs / N 得到的正好是 N/2+1 个点。

另外特别提醒一句：功率谱转对数坐标时用 10log10，因为已经是功率量纲；而幅度谱幅值转 dB 时才用 20log10。这个低级错误我在项目里见得太多了。

5.5 代码运行报错的快速定位

给几个排查思路，遇到问题可以按顺序查：

• 自相关函数 r 的长度是 M+1，别拿 M 当长度用。
• fft 的长度 nfft 至少要大于等于自相关序列长度，否则 MATLAB 会截断数据，输出结果完全不对。
• 用 abs(fft(x)).^2 时，先确认平方符号没有丢；丢平方以后画出来的是幅度谱开方量纲，峰都在但数字全不对。
• 多写几行打印语句，直接在命令行里看 length(f) 和 length(P)，长度对齐是画图不出错的第一步。
• 如果两个函数写在同一个脚本里，注意MATLAB版本要求脚本尾部函数支持，建议还是把 bt_psd 和 periodogram_my 分别保存成 .m 文件，结构更清晰，调试也方便。

最后说点我自己的习惯。我在实际项目里很少只用一种估计方法下结论，通常是先用周期图法快速看一遍全貌，再用BT法或者后续的Welch法针对感兴趣频段精调。两个结果如果差别很大，往往不是算法选错了，而是窗长度或者M没调好。你可以在手边准备一个简单的已知信号测试函数，每次改代码之后先跑一遍，确认谱峰位置和相对幅度都合理，再拿到真实数据上，能省掉很多排查时间。