功率谱密度的估计

首先功率谱密度是用来处理随机信号的。不同于确定性信号，随机信号不能通过一个确切的数学公式来描述，也不能准确地进行预测。因此，对随机信号一般只能在统计意义上来研究。

通常，随机信号是时间的函数，是无限长信号。在分析随机信号时，往往取某一段有限长（样本）信号用以研究随机信号的特征。对于各态遍历平稳随机信号，不同样本得到的时间平均和自相关函数是一样的，可用来表征总体信号的特征。

功率谱密度的估计包括非参数和参数化方法，本文主要介绍非参数化方法的功率谱估计及对比。
在 MATLAB 中，提供了 fft、 pwelch、periodogram、pmtm 等可用于非参数功率谱估计的函数，下面将对这些方法进行比较。

周期图法(periodogram)

[pxx, f] = periodogram(x, window, nfft, fs)

$P_{xx}(f)=\frac{1}{L{F}_s}{\left|\sum _{n=0}^{L-1}{x}_L(n)e^{-j2\pi fn/F_s}\right|}^2$
对于实数信号，除了0频率和Nyquist frequency频率点(fs/2)，其他点都要乘以2才能得到频率谱密度的估计。

• 周期图法 (periodogram) 与 fft 的比较

% 信号
fs = 1000;
t = (0:fs-1)/fs;
A = [1, 2];
f = [150;140];
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

% periodogram
subplot(211);
% periodogram 默认使用 Hamming 窗口，所以这里必须指定窗口类型为矩形窗，否则两种方式计算功率谱密度不完全一样
periodogram(xn, rectwin(length(xn)), length(xn), fs);

% FFT
subplot(212);
N = length(xn);
Y = fft(xn);
Y = Y(1:N/2+1);
xpsd = abs(Y).^2/(N*fs);
xpsd(2:end-1) = 2*xpsd(2:end-1);
f = 0:fs/N:fs/2;
plot(f, 10*log10(xpsd));
grid on;
xlabel('Frequency(Hz)');
ylabel('Power/Frequency(db/Hz)');
title('Periodogram Using FFT');

• 窗函数的作用
  前面用的是矩形窗(rectwin)，矩形窗口容易造成频谱泄露。对于非矩形窗口，被截断的信号的端点被平滑衰减，因此引入的杂散频率不那么严重。但另一方面，非矩形窗口会加宽主瓣，导致分辨率降低。

fs = 1000;
t = (0:fs/10)/fs;
A = [1, 2];
f = [150; 140];
nfft = 1024;
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

% 矩形窗
subplot(211);
periodogram(xn, rectwin(length(xn)), nfft, fs);

% Hamming窗，用汉明窗可达到的PSD分辨率大约是用矩形窗可达到的PSD分辨率的一半。
subplot(212);
periodogram(xn, hamming(length(xn)), nfft, fs);

可用的窗函数

周期图谱估计的方差不随数据长度N的增大而减小，而是近似于某一个值，说明了周期图谱估计不是一致估计。

提高(pwelch、pmtm)

1. [Pxx, f] = pwelch(x, window, noverlap, nfft, fs);

welch 方法将时间序列数据划分为多段（有重叠），用上面的方法计算每段的功率谱密度，然后平均PSD估计值。相对于整个数据的单个周期图估计，修改后的平均周期图估计方差要小。

fs = 1000;
t = (0:0.3*fs)/fs;
A = [2, 8];
f = [150; 140];
% 与之前不同，这里增大了噪声信号幅度
xn = A*sin(2*pi*f*t) + 5*randn(size(t));

% periodogram 矩形窗
subplot(211);
periodogram(xn, rectwin(length(xn)), 1024, fs);

% pwelch 矩形窗 分段越多，方差越小，频率分辨率越低，二者相互制约
subplot(212);
pwelch(xn, rectwin(150), 50, 1024, fs);

periodogram 方法中噪声和泄漏使正弦曲线中的一个无法分辨，而 pwelch 尽管分辨率低还是区分开了正弦信号与噪声。
用 pwelch 方法接着减小方差，频率分辨率减小，会发现另一个正弦信号也会「消失」。

pwelch(xn, rectwin(100), 75, 1024, fs);

2. [Pxx, f] = pmtm(x, nw, nfft, fs);

fs = 1000;
t = (0:fs)/fs;
A = [1, 2];
f = [150; 140];

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

pmtm(xn, 4, [], fs);

通过降低时间带宽乘积，可以以更大的方差为代价来提高分辨率。

pmtm(xn, 1.5, [], fs);