首先功率谱密度是用来处理随机信号的。不同于确定性信号,随机信号不能通过一个确切的数学公式来描述,也不能准确地进行预测。因此,对随机信号一般只能在统计意义上来研究。
通常,随机信号是时间的函数,是无限长信号。在分析随机信号时,往往取某一段有限长(样本)信号用以研究随机信号的特征。对于各态遍历平稳随机信号,不同样本得到的时间平均和自相关函数是一样的,可用来表征总体信号的特征。
功率谱密度的估计包括非参数和参数化方法,本文主要介绍非参数化方法的功率谱估计及对比。
在 MATLAB 中,提供了 fft、 pwelch、periodogram、pmtm 等可用于非参数功率谱估计的函数,下面将对这些方法进行比较。
周期图法(periodogram)
[pxx, f] = periodogram(x, window, nfft, fs)
P
x
x
(
f
)
=
1
L
F
s
∣
∑
n
=
0
L
−
1
x
L
(
n
)
e
−
j
2
π
f
n
/
F
s
∣
2
P_{xx}(f)=\color{red}\frac{1}{LF_s}\color{black}\left|\displaystyle\sum_{n=0}^{L-1} x_L(n)e^{-j2\pi fn/F_s}\right|^2
Pxx(f)=LFs1∣∣∣∣∣n=0∑L−1xL(n)e−j2πfn/Fs∣∣∣∣∣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);