随机信号谱估计方法的Matlab实现
摘要：
功率谱估计是随机信号分析中的一个重要内容。从介绍功率谱的估计原理入手分析经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。经典功率谱估计的方差大、谱分辨率差，分辨率反比于有效信号的长度，但现代谱估计的分辨率不受此限制。给出了功率谱估计的应用。
关键词：功率谱估计；周期图法；AR参数法；
1 引言
在一般工程实际中，随机信号通常是无限长的，例如，传感器的温漂，不可能得到无限长时间的无限个观察结果来获得完全准确的温漂情况，即随机信号总体的情况，一般只能在有限的时间内得到有限个结果，即有限个样本，根据经验来近似地估计总体的分布。有时，甚至不需要知道随机信号总体地分布，而只需要知道其数字特征，如均值、方差、均方值、相关函数、功率谱的比较精确的情况即估计值。功率谱估计(PSD)是用有限长的数据估计信号的功率谱，它对于认识一个随机信号或其他应用方面都是重要的，是数字信号处理的重要研究内容之一。功率谱估计可以分为经典谱估计(非参数估计)和现代谱估计(参数估计)。
2 .平均周期图法和平滑平均周期图法
对于周期图的功率谱估计, 当数据长度N 太大时, 谱曲线起伏
加剧, 若N 太小, 谱的分辨率又不好,因此需要改进。两种改进的估
计法是平均周期图法和平滑平均周期图法。
(1)Bartlett 法:
Bartlett 平均周期图的方法是将N 点的有限长序列x(n)分段求周期图再平均。
Matlab 代码示例1:
fs=600;
n=0:1/fs:1;
xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n));
nfft=512;
window=hamming(nfft); %矩形窗
noverlap=0;%数据无重叠
p=0.9;%置信概率
[Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p);
index=0:round(nfft/2- 1);
k=index*fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);
figure(2)
plot(k,[plot_Pxx plot_Pxx- plot_Pxxc
plot_Pxx+plot_Pxxc]);
matlab调试图下图
(2)Welch 法:
Welch 法对Bartlett 法进行了两方面的修正, 一是选择适当的窗函数w(n), 并在周期图计算前直接加进去, 加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时, 可使各段之间有重叠,这样会使方差减小。
Matlab 代码示例4:
Fs=600;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));
nfft=512;
window=boxcar(100);%矩形窗
window1=hamming(100);%海明窗
window2=blackman(100);%blackman 窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2], 计算一半的频率[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range); [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range); [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range); plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(1)
plot(f,plot_Pxx);
title('加矩形窗');
figure(2)
plot(f,plot_Pxx1);
title('加海明窗');
figure(3)
plot(f,plot_Pxx2);
title('加blackman 窗');
matlab调试如下
(3)AR模型法
经典谱的主要缺点是频率分辨率低。这是由于周期图法在计算中把观测到的有限长的N个数据以外的数据认为是零，这显然与事实不符。如果把以观察到的为N个以外的数据全为零，就有可能克服经典谱估计的缺点。
一个实际中的随机过程总是可以用以下模型很好的表示：
当除外的所有，均为零时的形式称为p阶自回归模型即AR模型，又称为全极点模型。当方差为的白自噪声通过AR模型时，输出的功率谱密度为：
在Matlab 仿真中可调用Pburg函数直接画出基于burg 算法的功率谱估计的曲线图。
Matlab 代码示例:
用周期图法求出的功率谱曲线和burg 算法求出的AR 功率谱曲线( p=50)
fs=200;
n=0:1/fs:1;
xn =cos (2*pi*40*n)+cos (2*pi*41*n)+3*cos
(2*pi*90*n)+0.1*randn(size(n));
window=boxcar(length(xn));
nfft=512;
[pxx,f]=periodogram(xn,window,nfft,fs);
subplot(121)
plot(f,10*log10(pxx))
xlabel(' frequency(hz)' );
ylabel(' power spectral density(Db/Hz)' );
title(' periodogram psd estimate' );
order1=50;
range=' half' ;
magunits=' db' ;
subplot(122)
