企业官网建设流程全解析

说实话，我刚入行做信号处理那会儿，拿到一段振动数据下意识就是plot(abs(fft(x)))，看着谱峰一阵兴奋。直到有一次做功率谱估计的项目，同一个测点连续采了十段数据，每段幅值谱长得都不一样，毛刺多到下不去眼。带我的老师傅只说了一句：你要看的是功率在频率上怎么分布，不是瞬时的幅值。这句话直接把我推向了BT法和周期图法这条经典路线。

这篇是这个系列的第一篇，先把谱估计里两条最古老的路线讲透——BT法（Blackman–Tukey，也叫自相关法/间接法）和周期图法（也叫直接法）。这两兄弟构成了经典谱估计的骨架，后边不管你是学Welch法、加窗平均，还是搞现代参数谱估计、子空间方法，基本思维都能从这里顺下来。我会把原理、Matlab代码、尺度和工程坑一次讲清楚，代码我试过能直接跑，建议你亲手复现一遍。

1. 为什么幅值谱在随机信号面前会“失灵”

1.1 确定性信号做FFT没毛病，问题出在随机信号

先看一段理想正弦波。给它做FFT，频谱上只有一根干净谱线，幅值、频率、相位都能读出来。换成两个正弦叠加，峰照样干干净净。这类确定性信号有确定的数学表达式，FFT相当于把一个有限长的时域快照变换到频域，结果稳定、可复现。

但实际工程里的信号几乎全是随机过程：轴承振动、声学噪声、海面回波、通信信道干扰。随机信号没有确定的表达式，每次观测只是一个“样本实现”。你今天采一段、明天采一段，波形不同，直接做FFT得到的“幅值谱”当然也不同。

就好比你站在一条河边测水流速度，如果水流是平稳的，你舀起一瓢水测出的流速可以作为整条河的近似；如果河面波浪很大，一瓢水测出的数就毫无代表性。随机信号每一次观测的FFT都相当于只舀了一瓢水，瓢瓢都不一样的。

1.2 猜一下：周期的“幅值谱”方差有多大？

你可以自己做个实验。生成一段高斯白噪声，长度N=1024，连续做50次FFT，把50条幅值谱叠加画在一起。会发现每条谱线都像刺猬一样到处乱跳，峰值差异能有好几个dB。

更反直觉的是：把N从1024提高到2048、4096，谱线的毛刺并不会变少。这跟确定性信号的FFT完全不同，确定性信号加长数据后谱线只会更细更稳定，随机信号做直接FFT，谱估计的方差基本不随N减少，教科书里称这种现象为“不一致估计”。这就是为什么我们不能只靠一次FFT就去描述随机信号的频域特性。

1.3 功率谱密度：换一个统计量再出发

既然随机信号的幅值不稳定，那工程界的标准做法就是放弃看幅值，改看“功率在频率上怎么分布”，也就是功率谱密度（PSD，Power Spectral Density）。

定义很简单，一个平稳随机过程的功率谱密度是自相关函数的傅里叶变换，这就是著名的维纳–辛钦定理：

[ S_{xx}(f) = \int_{-\infty}^{+\infty} R_{xx}(\tau) e^{-j 2\pi f \tau} d\tau ]

其中自相关函数定义为：

[ R_{xx}(\tau) = E[x(t) \cdot x(t+\tau)] ]

它衡量的是信号与自身延迟复制品之间的相似程度。如果信号里有个周期性分量，那它在某个滞后处就会显著地“和自己相似”，自相关函数会在那个位置附近出现振荡特征，再一傅里叶变换，就还原成频域里的谱峰。

频域里的功率谱，它的积分等于信号平均功率：

[ P = \int_{-\infty}^{+\infty} S_{xx}(f) df = E[|x(t)|^2] ]

所以功率谱密度是“每赫兹里有多少功率”，单位通常是V²/Hz或者g²/Hz，而日常我们习惯取10倍对数画成dB，就得到熟悉的谱图。

明白了这一点，就自然引出两条路线：既然功率谱是自相关函数的傅里叶变换，那我们可以先估计自相关再做变换，这就是BT法；也可以直接对数据做FFT再取模平方去逼近功率谱，这就是周期图法。

2. BT法与周期图法的数学本质：从两个方向逼近同一个东西

2.1 BT法：先估计自相关，再傅里叶变换

BT法也叫间接法，步骤非常朴素：

1. 对观测数据 (x(0),...,x(N-1))，估计出自相关函数序列 (\hat{R}_{xx}(\tau))；
2. 对自相关序列做傅里叶变换（实际是DTFT或FFT），得到功率谱密度。

其中自相关用的是“有偏估计”：

[ \hat{R}{xx}(\tau) = \frac{1}{N} \sum{n=0}^{N-1-|\tau|} x(n+|\tau|) x(n) ]

之所以分母取N而不是 (N-|\tau|)，是为了让估计的均方误差更小，这件事在后面“工程坑”那一章还会细说。

得到自相关序列后，一般还要对它加一个窗 (w(\tau))，截取其中一段，于是BT法的谱估计写作：

[ S_{BT}(f) = \sum_{\tau=-(N-1)}^{N-1} \hat{R}_{xx}(\tau) w(\tau) e^{-j 2\pi f \tau / f_s} ]

窗函数在这里很关键。如果你不截断、不加窗，BT法就不是一个独立的方法，而是周期图法的另一种写法；一旦主动截短，比如只取滞后范围 (|\tau| \le M)，谱线就会被平滑，方差变小，但频率分辨率也被拉低。这就是BT法的核心思想：用可控的偏差去换更小的方差。

2.2 周期图法：直接对数据做FFT再模平方

周期图法直接得多。给定N点数据，功率谱估计是：

[ S_{per}(f) = \frac{1}{N f_s} \left| \sum_{n=0}^{N-1} x(n) e^{-j 2\pi f n / f_s} \right|^2 ]

在Matlab里实现就是一个fft取模平方，再归一化。注意分母里的 (N f_s) 是把“离散傅里叶变换的求和平方”换算成“单位赫兹的功率密度”，这一步很多初学者会漏，漏掉之后谱的高低就完全不对了。

周期图法胜在简单直观，一度是工程上用得最多的谱估计方法，尤其配合FFT之后计算效率极高。但它的缺点也很致命：方差大。前面说的随机信号谱毛刺的问题，在周期图法里体现得淋漓尽致。

2.3 两者其实是远房亲戚：周期图也是一种BT法

想深一层，把周期图公式展开：

[ S_{per}(f) = \frac{1}{N} \sum_{n=0}^{N-1} \sum_{m=0}^{N-1} x(n) x(m) e^{-j 2\pi f (n-m)/f_s} ]

令 (\tau = n-m)，交换求和顺序后，这个式子等价于“对所有滞后点的有偏自相关估计不做任何截断地做DTFT”。换句话说，周期图法就是BT法在窗函数取全通、自相关取全部滞后时的特殊情况。

那为什么实际算出来两条谱线长得不一样？因为工程里BT法一定会截断自相关，而截断本身等效于在频域用窗函数做卷积，谱被平滑了。所以你会看到：周期图法毛刺多、峰尖；BT法曲线更顺滑、峰更胖。

既然明白了这两种方法的本质，下面给出可运行的Matlab代码，你自己跑一遍就能直观感受它们的区别。

3. Matlab代码实战：用一组仿真信号看清两种方法的脾气

3.1 仿真信号为什么要这样设计

我做仿真的习惯是：先构造一个“我知道答案”的信号，再把算法丢上去，看它能还原几分。这里的信号设计是这样的：

• 采样率 (f_s = 1024) Hz，数据长度 (N = 2048)；
• 两个正弦分量：75 Hz（幅值1.0）和 78 Hz（幅值0.6），间隔只有3 Hz；
• 叠加高斯白噪声，标准差0.8，信噪比不算高。

为什么要让两个频率这么近？因为频率分辨率极限大约等于 (1/T)，这里 (T=N/f_s=2) 秒，所以极限分辨率约0.5 Hz。3 Hz的间隔理论上能分开，但加上噪声和窗效应以后，能不能稳定地识别出两个峰，就非常考验算法了。这个设计能把周期图法的“毛刺”和BT法的“平滑”同时暴露出来。

3.2 完整可运行的Matlab代码

%% 功率谱估计（一）：BT法与周期图法对比 clear; close all; clc; rng(1); % 固定随机种子，方便复现 %% 1. 构造仿真信号 fs = 1024; % 采样率 Hz N = 2048; % 数据长度 t = (0:N-1) / fs; f1 = 75; A1 = 1.0; f2 = 78; A2 = 0.6; x = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t) + 0.8*randn(1, N); %% 2. 周期图法（直接法） NFFT = 4096; % 补零到4倍，只做谱线插值 X = fft(x, NFFT); % 对x补零后做FFT P_per = abs(X).^2 / (N * fs); % 双边功率谱密度，单位 V^2/Hz f = (0:NFFT-1) * fs / NFFT; % 频率轴 P_per = P_per(1:NFFT/2 + 1); % 只取单边 P_per(2:end-1) = 2 * P_per(2:end-1); % 单边化：0和Nyquist点不乘2 %% 3. BT法（自相关法/间接法） Mmax = N / 4; % 最大滞后点数，经验值是N/4左右 Rxx = xcorr(x, Mmax, 'biased'); % 有偏自相关估计，长度 2*Mmax+1 P_bt = abs(fft(Rxx, NFFT)) / fs; % 对自相关做FFT，得到双边PSD P_bt = P_bt(1:NFFT/2 + 1); P_bt(2:end-1) = 2 * P_bt(2:end-1); % 单边化 %% 4. 功率积分校验：PSD积分应接近信号平均功率 df = fs / NFFT; fprintf('信号平均功率 from x : %.3f\n', mean(x.^2)); fprintf('周期图法总功率积分 : %.3f\n', sum(P_per) * df); fprintf('BT法总功率积分 : %.3f\n', sum(P_bt) * df); %% 5. 画图对比 figure('Color', 'w', 'Position', [100 100 900 500]); plot(f, 10*log10(P_per), 'LineWidth', 1); hold on; plot(f, 10*log10(P_bt), 'LineWidth', 1); xline(f1, 'k--'); xline(f2, 'k--'); xlim([0 150]); ylim([-80 20]); xlabel('频率 / Hz'); ylabel('功率谱密度 / dB'); legend('周期图法', 'BT法', '真实频率', 'Location', 'best'); grid on; title('BT法与周期图法对比（N=2048, M=N/4, 补零到4096点）');

3.3 运行结果怎么解读

这段代码跑完，你会看到两条典型的谱线：

• 周期图法的曲线毛刺非常明显，整体像锯齿状，但75 Hz和78 Hz两个峰尖依然清晰，两个峰之间有一个明显的“谷”。这说明周期图法对频率的定位更锐利，代价是噪声带来的毛刺。
• BT法的曲线明显更顺滑，毛刺被“熨平”了很多，但75 Hz和78 Hz两个峰的轮廓也变胖了，小的那个峰（78 Hz）可能被大峰的裙边盖住一部分，需要仔细看才能认出是两个峰。

这就是分辨率与方差的跷跷板：BT法用M=N/4截断自相关，等效窗变短，主瓣变宽，所以两个靠近的峰变得难分，但代价是估计方差变小，谱图更干净。

控制台还会输出三行功率校验数值。正常情况下总功率应该比较接近，比如信号平均功率在1.8左右，周期图法和BT法的积分也会落在附近。BT法如果M取得太小，总功率会明显偏低，因为长滞后的自相关被截掉了，细碎的高频功率丢失了一部分。看到这种偏差不必惊慌，这是BT法的固有偏差。

4. 工程里真正暗藏杀机的是尺度和窗函数

4.1 频率轴映射：f、NFFT、df到底是什么关系

初学者最容易栽在频率轴上。记住这套规则：

• (\Delta f = f_s / N_{\text{FFT}})，这是频率轴的最小间隔；
• 频率轴写成f = (0:NFFT-1) * fs / NFFT；
• 如果只画单边谱，取前 (N_{\text{FFT}}/2 + 1) 个点，并且中间点要乘2，因为负频率的功率被折到正频率来了；
• 0 Hz（直流）和 Nyquist频率（(f_s/2)）这两个点不能再乘2，因为它们没有对应的负频率分量。

为什么乘2？因为实信号的双边谱关于0 Hz对称，正负频率各占一半功率，用单边谱表示时自然要把另一半加回来。漏乘2是新手最常见的尺度错误，谱图形状看着对，但功率校验时对不上。

4.2 xcorr的 'biased' 和 'unbiased'：无偏的那个其实很坑

在BT法里，自相关估计我强调用xcorr(x, Mmax, 'biased')。很多人在这一步会困惑：有偏估计听起来不严谨，要不要用'unbiased'？

答案是：不要。xcorr的'unbiased'模式在每个滞后点除以 (N-|\tau|)，也就是用“真实可用的样本数”做平均，这样每个滞后点的估计都是无偏的，但滞后越大、参与平均的样本数越少，方差被放得越大。你会在自相关序列的尾巴上看到剧烈的抖动，这些抖动在FFT之后变成谱里的高频毛刺。

而'biased'模式始终除以N，长滞后的估计虽然引入了一定的偏差，但方差更稳定。数学上，有偏自相关估计的一致性更好。在数据长度有限时，用有偏估计做谱估计是工程上的默认选择。

你可以自己对比一下：把Mmax设成接近N，分别用'biased'和'unbiased'算谱，'unbiased'那版尾部噪声会大到一个可怕的程度。

4.3 BT法里M到底取多大

BT法的截断点数M是一个绕不开的参数。M太大，接近N，等于没截断，方差问题回到周期图法；M太小，自相关只留了一小段，频域窗变得很胖，峰分辨率急剧恶化。

我个人的经验值参考：

• 如果谱图以看“有没有峰、峰在哪”为目标：M取 (N/5) 到 (N/4) 比较适合，谱线平滑，还能保住主要峰位；
• 如果目标是把两个很近的峰彻底分辨出来：M取 (N/2) 以上，甚至直接用周期图法，不要为了平滑牺牲分辨率；
• M取 (N/10) 就会明显看到谱峰变宽，只适合看低频的大致包络。

记住一点：M本质上决定了频域窗的主瓣宽度，窗口越长，主瓣越窄、分辨率越高，但方差也越大。没有免费的午餐。

4.4 加窗的隐蔽代价：主瓣展宽和幅度损失

无论周期图还是BT法，只要数据不是整周期截断，就逃不过窗函数的影响。很多人以为加窗只是为了减少频谱泄漏，却忽略了两个副作用。

第一个副作用是主瓣展宽：加窗等于时域乘窗、频域卷积，任何非矩形窗都会把理想谱线“抹”成一片有一定宽度的包络，于是原来尖锐的谱峰变胖。两个间隔很近的频率就可能糊成一坨。

第二个副作用是幅度损失：窗函数会把数据两端的样本压低，整体能量被衰减。比如Hamming窗的相干增益约为0.54，意味着如果你不加修正，直接算出来的谱峰值会比真实值低一大截。要修正的话，需要把谱除以窗函数的平均增益：

[ G_c = \frac{1}{N} \sum_{n=0}^{N-1} w(n) ]

这个修正值很容易被忘记。你用统一幅度的仿真正弦去标定系统时，算出来的峰总是偏低，大概率就是没除相干增益。

5. 经典谱估计的边界，以及我学到的两件事

5.1 分辨率墙是1/T，不是1/NFFT

很多刚接触谱估计的人会误以为，把FFT点数增大（补零）就能看到更细的峰。这是错觉。补零只是在已知的DTFT采样点上做插值，让曲线更平滑、读数更方便，但真正的“分辨能力”是由数据长度决定的：

[ \Delta f_{\text{res}} \approx \frac{1}{T} = \frac{f_s}{N} ]

两个频率分量的间隔小于这个值，再补多少零也救不回来。我见过不少新人在仿真里拼命补零，以为看到两个峰分开，换一段数据又是一锅粥，就是因为没理解分辨率墙在这里。

5.2 周期图方差大，不是FFT的锅

还有一个常见误解：周期图毛刺多是FFT算法精度不够。不是的，FFT算得又快又准，问题出在随机信号本身。周期图估计的方差大是它的固有属性，而且不随N增大而降低。

想压方差，思路就是“平均”。把一段长数据切成了K段，每一段做周期图，再把K个周期图平均起来，这就是Welch法的雏形。平均的段数越多，方差越小，但单段数据变短，分辨率又降下来了。说白了，还是在方差和分辨率之间找平衡。

5.3 工程上什么时候用经典法，什么时候该上现代法

我的习惯是：拿到一段信号，先跑一版周期图法和BT法，用它们快速把握频谱结构——有几个峰、大致在什么频段、噪声底多高。经典法的优势是简单、稳健、不需要设置模型阶数，几乎不会因为参数偏差输出离谱的结果。

只有碰到经典法明显无能为力的时候，比如两个正弦间距小于1/T、或者要求极高频域分辨率，我才会考虑现代谱估计路线，比如AR模型法、Capon法、MUSIC类子空间方法。这些方法用模型先验换分辨率，代价是模型阶数选错时结果会很夸张，需要更多经验和校验手段。

我自己在项目里的体会是：经典谱估计就像一把普通的尺子，量不准时候会让你抓狂，但它从不会骗你；现代谱估计是一把高精度卡尺，用得好能测出微米级差异，用不好量出来的数可能没有任何物理意义。所以建议你先踏踏实实把这把“普通尺子”用明白，再去碰更高级的手段。

下一篇我打算把Welch法展开讲，包括分段、重叠率、窗函数选择三个参数到底怎么搭配，以及在做故障诊断时怎么用Welch法稳定地提取边频带。到时候把这一篇里“方差与分辨率要平衡”的结论落到实处。