摘 要：功率譜估計是隨機(jī)信號分析中的一個重要內(nèi)容。從介紹功率譜的估計原理入手分析經(jīng)典譜估計和現(xiàn)代譜估計兩類估計方法的原理、各自特點(diǎn)及在Matlab中的實(shí)現(xiàn)方法。經(jīng)典功率譜估計的方差大、譜分辨率差，分辨率反比于有效信號的長度，但現(xiàn)代譜估計的分辨率不受此限制。給出了功率譜估計的應(yīng)用。

關(guān)鍵詞：功率譜估計;周期圖法;AR參數(shù)法;隨機(jī)信號分析

中圖分類號：TP274 文獻(xiàn)標(biāo)識碼：B 文章編號：1004373X(2008)1813203

Matlab Realization of Random Signal Spectrum Estimation

ZHANG Weiwei1，2，XING Wuqiang2

(1.Electronic Engineering College，Xidian University，Xi′an，710071，China;2.Xi′an Institute of Post Telecomunications，Xi′an，710120，China)

Abstract：Power Spectrum Estimation is an important research topic in the area of the random signal.The paper mainly introduces the principles of classical spectrum estimation and modern spectrum estimation，discusses characteristics of the methods of realization in Matlab.The variance of classical spectrum estimation is large and the spectral resolution is not good.Resolution is inverse to the length of the valid signal.But modern spectrum estimation doesn′t limit it.At last it gives the application of power spectrum estimation.

Keywords：power spectrum estimation;periodogram method;AR parameter method;random signal analysis

1 引 言

在一般工程實(shí)際中，隨機(jī)信號通常是無限長的，例如，傳感器的溫漂，不可能得到無限長時間的無限個觀察結(jié)果來獲得完全準(zhǔn)確的溫漂情況，即隨機(jī)信號總體的情況，一般只能在有限的時間內(nèi)得到有限個結(jié)果，即有限個樣本，根據(jù)經(jīng)驗來近似地估計總體的分布。有時，甚至不需要知道隨機(jī)信號總體地分布，而只需要知道其數(shù)字特征，如均值、方差、均方值、相關(guān)函數(shù)、功率譜的比較精確的情況即估計值。功率譜估計（PSD）是用有限長的數(shù)據(jù)估計信號的功率譜，它對于認(rèn)識一個隨機(jī)信號或其他應(yīng)用方面都是重要的，是數(shù)字信號處理的重要研究內(nèi)容之一。功率譜估計可以分為經(jīng)典譜估計（非參數(shù)估計）和現(xiàn)代譜估計（參數(shù)估計）。前者的主要方法有2種：一種是間接法；另一種是直接法。后者的主要發(fā)法有最大熵譜分析法（AR模型法）、Pisarenko諧波分解法、Prony提取極點(diǎn)法、Prony譜線分解法以及Capon最大似然法。其中周期突法和AR模型法是用得較多且最具代表性的方法。

2 周期圖法

由于對信號做功率譜的估計，需要用計算機(jī)實(shí)現(xiàn)，如果是連續(xù)信號，則需要變換為離散信號。下面討論離散隨機(jī)信號序列的功率譜估計問題。

連續(xù)時間隨機(jī)信號的功率譜密度與自相關(guān)函數(shù)是一對傅里葉變換對，即:Sx(Ω)=∫∞－∞Rx(τ)e－jΩτdτ

若Rx(m)是Rx(τ)的抽樣序列，由序列的傅里葉變換的關(guān)系，可得:Sx(ejω)=∑∞m=－∞Rx(m)e－jωn

即Sx(ejω)與Rx(m)也是一對傅里葉變換對。顯然，由序列傅里葉的頻譜特性可知Sx(ejω)是以2π為周期的。而實(shí)際計算只能從離散隨機(jī)信號序列x(n)的有限長（長度為N）的數(shù)據(jù)來對Sx(ejω)與Rx(m)進(jìn)行估計。設(shè)有限長離散序列為xN(n)，則:xN(m) = 1N［xN (m)*xN (-m)］

xN(ejω) =DFT［RxN(m)］

由DFT的下列卷積特性：

若X(ejω)=DFT［xN(m)］，則:X(ejω)=DFT［xN(－m)］

從而:DTFT［xN(m)］ = 1NDFT［xN (m)］DFT［xN (-m)］

即:xN(ejω)=1NX(ejω)X(ejω)=1N|X(ejω)|2

綜上所述，先用FFT求出宿疾隨機(jī)離散信號N點(diǎn)的DFT，再計算幅頻特性的平方，然后除以N，即得出該隨機(jī)信號得功率譜估計。由于這種估計方法在把Rx(τ)離散化的同時，使其功率譜周期化，故稱之為“周期圖法”，也稱為經(jīng)典譜估計方法。周期圖法進(jìn)行譜估計，是有偏估計，由于卷積的運(yùn)算過程會導(dǎo)致功率譜真實(shí)值的尖峰附近產(chǎn)生泄漏，相對地平滑了尖峰值，因此造成譜估計的失真。另外，當(dāng)N→∞時，功率譜估計的方差不為零，所以不是一致性估計。并且功率譜估計在ω等于2π/N整數(shù)倍的各數(shù)字頻率點(diǎn)互不相關(guān)。其譜估計的波動比較顯著，特別是當(dāng)N越大、2π/N越小時，波動越明顯。但如果N取得太小，又會造成分辨率的下降。

Matlab示例：

Fs=2000;

N=1024;

Nfft=1024;

n1=0:N-1;

t1=n1/Fs;

f1=50;

f2=100;

xn1=sin(2*pi*f1*t1)+2*sin(2*pi*f2*t1)+randn(1，N);

Pxx1=10*log10(abs(fft(xn1，Nfft).^2)/N);

f1=(0:length(Pxx1)-1)*Fs/length(Pxx1);

subplot(2，1，1)

plot(f1，Pxx1)

圖1 Matlab示例13 AR模型法

經(jīng)典譜的主要缺點(diǎn)是頻率分辨率低。這是由于周期圖法在計算中把觀測到的有限長的N個數(shù)據(jù)以外的數(shù)據(jù)認(rèn)為是零，這顯然與事實(shí)不符。如果把以觀察到的數(shù)據(jù)估計出一白噪聲激勵一物理網(wǎng)絡(luò)所形成，就不必認(rèn)為N個以外的數(shù)據(jù)全為零，就有可能克服經(jīng)典譜估計的缺點(diǎn)。

一個實(shí)際中的隨機(jī)過程總是可以用以下模型很好的表示:H(z)=∑pi=1biz－i1+∑pk=1anz－k

當(dāng)除b0外的所有bi均為零時的形式稱為p階自回歸模型即AR模型，又稱為全極點(diǎn)模型。

當(dāng)方差為σ2的白噪聲通過AR模型時，輸出的功率譜密度為:Pxx(ω)=σ2ω1+∑pk=1ake－jωk2

若已知參數(shù)a1，a2，…，ap及σ2ω，就可以得到信號的功率譜估計。它們之間是YuleWalker方程。解這個方程是一個復(fù)雜的數(shù)學(xué)問題，這里不做討論。

Matlab示例：

Fs=2000;

n=0:1/Fs:1;

xn=sin(2*pi*50*n)+2*sin(2*pi*100*n)+randn(size(n));

nfft=1024;

window=boxcar(100);

noverlap=20;

range=′half′;

［Pxx，f］=pwelch(xn，window，noverlap，nfft，F(xiàn)s，range);

plot_Pxx=10*log10(Pxx);

figure(1)

plot(f，plot_Pxx);

圖2 Matlab示例24 功率譜估計的應(yīng)用

功率譜估計可應(yīng)用于估計系統(tǒng)傳遞函數(shù)，在理想情況下，系統(tǒng)輸入x(t)，輸出y(t)及單位沖激響應(yīng)h(t)有如下關(guān)系:y(t)=h(t)x(t)

上述各量的傅里葉變換分別為X(Ω)，Y(Ω)和H(Ω)，則系統(tǒng)的頻響為:H(Ω)=Y(Ω)X(Ω)

利用Ω與s的關(guān)系，即可得到系統(tǒng)傳遞函數(shù)。

但在實(shí)際中，輸入與輸出都可能存在噪聲的干擾，即其實(shí)際輸入xsn(t)、實(shí)際輸出ysn(t)應(yīng)表示為:xsn(t)=x(t)+nx(t)

ysn(t)=yx(n)+yxn(t)+ny(t)

在上式中，nx(t)，ny(t)分別為輸入/輸出端引入的噪聲；yxn(t)是nx(t)引起的響應(yīng)。這時的系統(tǒng)頻響就相應(yīng)變?yōu)?Hsn(Ω)=Yx(Ω)+Yxn(Ω)+Ny(Ω)X(Ω)+Nx(Ω)

上式中的各量是前面2式中各量的傅里葉變換。這樣，由于輸入/輸出端的噪聲影響的不確定性將帶來計算傳遞函數(shù)的誤差。為此，可以應(yīng)用功率譜計算的方法估計系統(tǒng)傳遞函數(shù)，先在時域上進(jìn)行相關(guān)運(yùn)算，由相關(guān)特性可知，由于噪聲的隨機(jī)性，在做時域相關(guān)計算時只要τ足夠長，輸入/輸出的噪聲自相關(guān)函數(shù)就為零。而隨機(jī)噪聲與有用信號之間一般不存在相關(guān)關(guān)系，它們的互相關(guān)函數(shù)也為零。根據(jù)功率譜和頻響函數(shù)間的數(shù)學(xué)關(guān)系，有以下關(guān)系成立：H(Ω)=Sxy(Ω)Sx(Ω)；|H(Ω)|2=Sy(Ω)Sx(Ω)Matlab示例：

N=1024;

Nfft=256;

window=hanning(256);

noverlap=128;

Fs=1000;

n=0:N-1;

t=n/Fs;

randn(′state′，0)

xn=sin(2*pi*50*t)+randn(1，N);

h=ones(1，20)/20;

yn=filter(h，1，xn);

［Txy，f］=tfe(xn，yn，Nfft，F(xiàn)s，window，noverlap);

H=freqz(h，1，f，F(xiàn)s);

subplot(2，1，1)

plot(f，abs(H));

axis(［0，500，0，1］);

grid

subplot(2，1，2)

plot(f，abs(Txy));

axis(［0，500，0，1］);

grid

5 結(jié) 語

通過實(shí)驗仿真可以直接看出以下特性：經(jīng)典功率譜估計的方差大，譜分辨率差，分辨率反比與有效信號的長度，但現(xiàn)代譜估計的分辨率不受此限制。這是因為對于給定的N點(diǎn)有限長序列x(n)，雖然其估計出的自相關(guān)函數(shù)也是有限長的，但是現(xiàn)代譜估計的一些隱含著數(shù)據(jù)和自相關(guān)函數(shù)的外推，使其可能的長度超過給定的長度，不像經(jīng)典譜估計那樣受窗函數(shù)的影響。因而現(xiàn)代譜的分辨率較高，而且現(xiàn)代譜線要平滑得多，這些從圖3可以清楚地看出。

圖3 Matlab示例3

參 考 文 獻(xiàn)

［1］陸華光，彭學(xué)愚.隨機(jī)信號處理\\.西安：西安電子科技大學(xué)出版社，2003.

［2］周浩敏.信號處理技術(shù)基礎(chǔ)\\.北京：北京航空航天大學(xué)出版社，2002.

［3］胡廣書.數(shù)字信號處理\\.北京：清華大學(xué)出版社，2003.

［4］劉松強(qiáng).數(shù)字信號處理系統(tǒng)及其應(yīng)用\\.北京：清華大學(xué)出版社，2000.

［5］樓順天.基于Matlab的系統(tǒng)分析與設(shè)計信號處理\\.西安：西安電子科技大學(xué)出版社，1998.

［6］飛思科級.Matlab 6.5輔助圖像處理\\.北京：電子工業(yè)出版社，2003.

［7］朱仁峰.精通Matlab 7\\.北京：清華大學(xué)出版社，2006.

［8］黃志宇，劉保化，陳高平，等.隨機(jī)信號的功率譜估計及Matlab的實(shí)現(xiàn)\\.現(xiàn)代電子技術(shù)，2002，25(3):2123.

［9］佚名.經(jīng)典譜估計方發(fā)的Matlab分析\\.武漢：華中理工大學(xué)，2000.

［10］魏鑫，張平.周期圖法功率譜估計中的窗函數(shù)分析\\.現(xiàn)代電子技術(shù)，2005，28(3):2021.

作者簡介 張薇薇 女，1978年出生，陜西戶縣人，西安電子科技大學(xué)電子工程學(xué)院研究生，西安郵電學(xué)院繼續(xù)教育學(xué)院助教。

注：本文中所涉及到的圖表、注解、公式等內(nèi)容請以PDF格式閱讀原文