功率谱密度的估计

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

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

功率谱密度的估计包括非参数参数化方法,本文主要介绍非参数化方法的功率谱估计及对比。
在 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)=LFs1n=0L1xL(n)ej2πfn/Fs2
对于实数信号,除了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);

在这里插入图片描述

  • 7
    点赞
  • 67
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值