隨機信號譜估計方法的Ｍａｔｌａｂ實現

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

關鍵詞：功率譜估計;周期圖法;AR參數法;隨機信號分析

中圖分類號：TP274 文獻標識碼：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 引 言

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

2 周期圖法

由于對信號做功率譜的估計，需要用計算機實現，如果是連續信號，則需要變換為離散信號。下面討論離散隨機信號序列的功率譜估計問題。

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

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

即Sx(ejω)與Rx(m)也是一對傅里葉變換對。顯然，由序列傅里葉的頻譜特性可知Sx(ejω)是以2π為周期的。而實際計算只能從離散隨機信號序列x(n)的有限長（長度為N）的數據來對Sx(ejω)與Rx(m)進行估計。設有限長離散序列為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求出宿疾隨機離散信號N點的DFT，再計算幅頻特性的平方，然后除以N，即得出該隨機信號得功率譜估計。由于這種估計方法在把Rx(τ)離散化的同時，使其功率譜周期化，故稱之為“周期圖法”，也稱為經典譜估計方法。周期圖法進行譜估計，是有偏估計，由于卷積的運算過程會導致功率譜真實值的尖峰附近產生泄漏，相對地平滑了尖峰值，因此造成譜估計的失真。另外，當N→∞時，功率譜估計的方差不為零，所以不是一致性估計。并且功率譜估計在ω等于2π/N整數倍的各數字頻率點互不相關。其譜估計的波動比較顯著，特別是當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模型法

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

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

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

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

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

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，Fs，range);

plot_Pxx=10*log10(Pxx);

figure(1)

plot(f，plot_Pxx);

圖2 Matlab示例24 功率譜估計的應用

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

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

利用Ω與s的關系，即可得到系統傳遞函數。

但在實際中，輸入與輸出都可能存在噪聲的干擾，即其實際輸入xsn(t)、實際輸出ysn(t)應表示為: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)引起的響應。這時的系統頻響就相應變為:Hsn(Ω)=Yx(Ω)+Yxn(Ω)+Ny(Ω)X(Ω)+Nx(Ω)

上式中的各量是前面2式中各量的傅里葉變換。這樣，由于輸入/輸出端的噪聲影響的不確定性將帶來計算傳遞函數的誤差。為此，可以應用功率譜計算的方法估計系統傳遞函數，先在時域上進行相關運算，由相關特性可知，由于噪聲的隨機性，在做時域相關計算時只要τ足夠長，輸入/輸出的噪聲自相關函數就為零。而隨機噪聲與有用信號之間一般不存在相關關系，它們的互相關函數也為零。根據功率譜和頻響函數間的數學關系，有以下關系成立：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，Fs，window，noverlap);

H=freqz(h，1，f，Fs);

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 結 語

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

圖3 Matlab示例3

參 考 文 獻

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

［2］周浩敏.信號處理技術基礎\\.北京：北京航空航天大學出版社，2002.

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

［4］劉松強.數字信號處理系統及其應用\\.北京：清華大學出版社，2000.

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

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

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

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

［9］佚名.經典譜估計方發的Matlab分析\\.武漢：華中理工大學，2000.

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

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

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