功率谱估计方法综述

  • 格式:doc
  • 大小:199.00 KB
  • 文档页数:7

下载文档原格式

  / 7
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

功率谱估计方法综述:

简介:随机信号的持续时间是无限长的,因此随机信号的总能量是无限的,因而随机过程的任意一个样本寒暑都不满足绝对可积条件,所以其傅里叶变换不存在。尽管随机信号的总能量是无限的,但其平均功率却是有限的,因此,要对随机信号的频域进行分析,应从功率谱出发进行研究才有意义。信号的功率谱密度描述随机信号的功率在频域随频率的分布。功率谱估计(PSD)是用有限长的数据来估计信号的功率谱,即利用给定的N个样本数据估计一个平稳随机信号的功率谱密度。

背景:功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。功率谱估计方法主要分为2大类:非参数化方法(又称经典功率谱估计)和参数化方法(又称现代功率谱估计)。非参数化方法有相关函数法(BT法)、周期图法、平均周期图法、平滑平均周期图法等;而参数化谱估计有R模型法、移动平均模型法(简称MA模型法)、自回归移动平均模型法(简称ARMA模型法)、最大熵谱分析法(AR模型法)、Pisarenko谐波分解法、Prony 提取极点法、Prony谱线分解法以及capon最大似然法等,由于涉及许多复杂数学计算,在此未作详细数学推导,以下介绍几种常用的功率谱估计方法

一、非参数化方法(经典法)

经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。

1、自相关法

又称相关函数法(BT法),根据维纳—辛钦定理:平稳随机过程的自相关函数和功率谱函数是一傅里叶变换对,对于平稳随机信号来说,其相关函数是确定性函数,故其功率谱也是确定的.这样可由平稳随机离散信号的有限个离散值,求出自相关函数,然后作Fourier变换,得到功率谱。由于随机序列{X(n)}的自相关函数R(n)=E[X(n)X(n+m)]定义在离散点m上,设取样间隔为错误!未找到引用源。,可将随机序列的自相关函数用连续时间函数表示为错误!未找到引用源。

等式两边取傅里叶变换,则随机序列的功率谱密度

错误!未找到引用源。

错误!未找到引用源。

错误!未找到引用源。

BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值错误!未找到引用源。。即

错误!未找到引用源。

其中错误!未找到引用源。可有式得到。

Fs=500; %采样频率

n=0:1/Fs:1; %产生含有噪声的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+ randn(size(n));

nfft=512;

cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数,matlab函数 xcorr(求自相关函数)unbiased无偏

CXk=fft(cxn,nfft); %对cxn(即自相关函数)进行快速傅里叶变换,nfft为周期Pxx=abs(CXk); %对CXk(频谱)取绝对值(为什么取绝对值)

index=0:round(nfft/2-1); %计算出各点对应的功率谱

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

figure(1)

plot(k,plot_Pxx);

2、周期图法

周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换,得x(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。Matlab代码示例2:

Fs=600;%采样频率

n=0:1/Fs:1;%产生含有噪声的序列

xn=cos(2*pi*40*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);%直接法计算功率谱密度,xn为功率谱密度信号,window为窗口,nfft为采样点数, fs采样频率

plot(f,10*log10(Pxx));

window=boxcar(length(xn));%矩形窗

nfft=1024;

[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法

figure(1)

plot(f,10*log10(Pxx));

3、平均法

即Bartlett平均周期图法,是将N点的有限长序列x(n)分段求周期图再平均.将长度为N的数据分为L段,先对每段数据用周期图法进行谱估计,然后对L段求平均得到长度为N的数据的功率谱.

平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少,达到一致估计的目的,有一个定理:如果错误!未找到引用源。是不相关的随机变量,且都有个均值错误!未找到引用源。及其方差错误!未找到引用源。,则可以证明它们的算术平均的均值为错误!未找到引用源。。

即:平均法将错误!未找到引用源。的N个数据分成L段(N=ML),若各数据段相互独立,则平方后估计量的方差是原来不分段估计量方差的错误!未找到引用源。。所以当错误!未找到引用源。时,估计量的方差趋于0,达到一致估计的目的。但是,随着分段数L的增加,M点数减少,分辨率减少,使估计变成有偏估计。相反,若L减少,M增加,虽偏差减少,但方差增大。所以,在应用中,必须兼顾分辨率和方差的要求来适当选择M和L的值。Matlab代码

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]);