频谱分析(完整版)

1 2
对实信号，PSD 是关于直流信号对称的，所以 0 的 Pxx 就足够完整的描述 PSD 了。然而要获得整个 Nyquist 间隔上的平均功率，有必要引入单边 PSD 的概念：
0 0 Ponesided 2 Pxx 0
xL n x n wR n
因为时域的乘积等效于频域的卷积，所以上式的傅立叶变换是
1 XL f fs
fs / 2
fs / 2

X WR f d
前文中导出的表达式
ˆ f P xx
XL f fs L
2
说明卷积对周期图有影响。 正弦数据的卷积影响最容易理解。假设 x n 是 M 个复正弦的和
Parametric methods （参量类方法） 这类方法是假设信号是一个由白噪声驱动的线性系统的输 出。这类方法的例子是 Yule-Walker autoregressive (AR) method 和 Burg method 。这些方法先估计假设的 产生信号的线性系统的参数。 这些方法想要对可用数据相对较少的情况产生优于传统非参数 方法的结果。 Subspace methods （子空间类） 又称为 high-resolution methods （高分辨率法）或者 super-resolution methods （超分辨率方法） 基于对自相关矩阵的特征分析或者特征值分解产生信号的频率分量。 代表 方 法 有 multiple signal classificat ion (MUS IC) method 或 eigenvector (EV) method。这类方法对线谱（正弦信号的谱）最合适，对检测噪声下的正弦信号很有效，特 别是低信噪比的情况。 方法 周期图 Welch 多椎体 Yule-Walker AR Burg Covariance （协方差） 修正协方差 MUSIC 特征向量法 PSD 估计 重叠，加窗的信号段的平均周期图 多个正交窗（称为锥）的组合做谱估计 时间序 列的 估计 的自 相关函 数计 算自 回归 （AR）谱估计 通过最小化线性预测误差计算自回归（ AR） 谱估计 描述 函数 spectrum.periodogram, periodogram spectrum.welch, pwelch, cpsd, tfestimate, mscohere spectrum.mtm, pmtm spectrum.yulear, pyulear spectrum.burg, pburg
序列 xn 在整个 Nyquist 间隔上的平均功率可以表示为
Rxx 0
上式中的



s S xx S f d xx df 2 fs fs / 2
f /2
Pxx
S f S xx 以及 Pxx f xx 2 fs
被定义为平稳随机信号 xn 的 power spectral density (PSD)（功率谱密度） 一个信号在频带 1 , 2 ,0 1 2 上的平均功率可以通过对 PSD 在频带上积分 求出
S xx
m
R m e
xx

j m
注： S xx X ，其中 X lim
2
N
1 N
n N /2

N /2
xn e jn 。其 matlab
近似为 X=fft(x,N)/sqrt(N)，在下文中 X L f 就是指 matlab fft 函数的计算结果了 使用关系 2 f / f s 可以写成物理频率 f 的函数，其中 f s 是采样频率
-30 -40 -50 -60 -70 -80 -500
-400
-300
-200
-100 0 100 frequency/Hz
200
300
400
500
该图显示了一个主瓣和若干旁瓣，最大旁瓣大约在主瓣下方 13.5dB 处。这些旁瓣说 明了频谱泄漏效应。无限长信号的功率严格的集中在离散频率点 f k 处，而有限长信号在离 散频率点 f k 附近有连续的功率。 因为矩形窗越短，它的频率响应对 Dirac 冲击的近似性越差，所以数据越短它的频谱 泄漏越明显。考虑下面的 100 个采样的序列 randn('state',0) fs = 1000; t = (0:fs/10)/fs; A = [1 2]; f = [150;140]; % Sampling frequency % One-tenth of a second worth of samples % Sinusoid amplitudes % Sinusoid frequencies
% Sampling frequency % One second worth of samples % Sinusoid amplitudes (row vector)
f = [150;140]; % S inusoid frequencies (column vector) xn = A*sin(2*pi*f*t) + 0.1*randn(size(t)); 注意：最后三行表明了一个方便的表示正弦之和的方法，它等价于： xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t)); 对这个 PSD 的周期图估计可以通过产生一个周期图对象（periodogram object）来计算 Hs = spectrum.periodogram('Hamming'); 估计的图形可以用 psd 函数显示。 psd(Hs,xn,'Fs',fs,'NFFT',1024,'SpectrumType','twosided')
通过最小化前向预测误差做时间序列的自回 spectrum.cov, pcov 归（AR）谱估计 通过最小化前向及后向预测误差做时间序列 spectrum.mcov, pmcov 的自回归（AR）谱估计 多重信号分类 虚谱估计 spectrum.music, pmusic spectrum.eigenvector, peig

S xx f
m
R m e
xx
2 jfm / f s
相关序列可以从功率谱用 IDFT 变换求得：
Rxx m



s S xx e jm S f e2 jfm / f s d xx df 2 f s fs / 2
f /2
k 1 k 1
M
M
所以在有限长信号的频谱中，Dirac 函数被替换成了形式为 WR f f k 的项，该项对 应于矩形窗的中心在 f k 的频率响应。 一个矩形窗的频率响应形状是一个 sinc 信号，如下所示
矩形窗在物理频率上的功率谱密度 0 -10 -20
PSD dB watt/Hz
信号在频带 1 , 2 ,0 1 2 上的平均功率可以用单边 PSD 求出
2 1
P 1 ,2
P
onesided
d
频谱估计方法
Mat lab 信号处理工具箱提供了三种方法 Nonparametric methods （非参量类方法） PSD 直接从信号本身估计出来。最简单的就是 periodogram（周期图法） ，一种改进 的周期图法是 Welch's method 。更现代的一种方法是 multitaper method （多椎体法） 。
Matlab 信号处理工具箱 帮助文档 谱估计专题
翻译：无名网友 & Lyra
频谱分析
Spectral estimation（谱估计）的目标是基于一个有限的数据集合描述一个信号的功 率（在频率上的）分布。功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信 号的检测。 从数学上看，一个平稳随机过程 xn 的 power spectrum（功率谱）和 correlation sequence（相关序列）通过 discrete-time Fourier transform（离散时间傅立叶变换） 构成联系。从 normalized frequency（归一化角频率）角度看，有下式
Nonparametric Methods 非参数法
下面讨论 periodogram, modified periodogram, Welch, 和 multitaper 法。同时 也讨论 CPSD 函数，传输函数估计和相关函数。
Periodogram 周期图法
一个估计功率谱的简单方法是直接求随机过程抽样的 DFT ，然后取结果的幅度的平方。 这样的方法叫做周期图法。 一个长 L 的信号 xL n 的 PSD 的周期图估计是
M
x n Ak e jk n
k 1
其频谱是
X f f s Ak f f k
k 1
M
对一个有限长序列，就变成了
XL f
1 fs
fs / 2
fs / 2

f s Ak f k WR f d AkWR f f k
P 1 ,2
2 1
Pxx d
1 2

Pxx d
从上式中可以看出 Pxx 是一个信号在一个无穷小频带上的功率浓度，这也是为什么 它叫做功率谱密度。 PSD 的单位是功率（e.g 瓦特）每单位频率。在 Pxx 的情况下，这是瓦特/弧度/抽 或只是瓦特/弧度。在 Pxx f 的情况下单位是瓦特/赫兹。PSD 对频率的积分得到的单位是 瓦特，正如平均功率 P , 所期望的那样。
ˆ f P xx k
其中
X L fk fs L
2
， fk
kf s , k 0,1,, N 1 N
X L f k xL n e2 jkn / N
n 0
L 1
选择 N 是大于 L 的下一个 2 的幂次是明智的， 要计算 X L f k 我们直接对 xL n 补零到 长度为 N。假如 L>N，在计算 X L f k 前，我们必须绕回 xL n 模 N。 作为一个例子，考虑下面 1001 元素信号 xn ，它包含了 2 个正弦信号和噪声 randn('state',0); fs = 1000; t = (0:fs)/fs; A = [1 2];
Power Spectral Density Estimate via Periodogram 0 -10 -20
Power/frequency (dB/Hz)
-30 -40 -50 -60 -70 -80
0
0.1
0.2
0.3
0.4 0.5 0.6 Freqwenku.baidu.comency (kHz)
0.7
0.8
0.9
ˆ f P xx
XL f fs L
2
注：这里 X L f 运用的是 matlab 里面的 fft 的定义不带归一化系数，所以要除以 L 其中
X L f xL n e2 jfn / f s
n 0
L 1
实际对 X L f 的计算可以只在有限的频率点上执行并且使用 FFT。实践上大多数周期 图法的应用都计算 N 点 PSD 估计
平均功率通过用下述求和去近似积分 求得 [Pxx,F] = psd(Hs,xn,fs,'twosided'); Pow = (fs/length(Pxx)) * sum(Pxx) Pow = 2.5059 你还可以用单边 PSD 去计算平均功率 [Pxxo,F] = psd(Hs,xn,fs,'onesided'); Pow = (fs/(2*length(Pxxo))) * sum(Pxxo) Pow = 2.5011 周期图性能 下面从四个角度讨论周期图法估计的性能：泄漏，分辨率，偏差和方差。 频谱泄漏 考虑有限长信号 xL n ，把它表示成无限长序列 x n 乘以一个有限长矩形窗 wR n 的 乘积的形式经常很有用：