基于Matlab的信號平穩性檢驗系統

摘 要:信號的平穩性檢驗在信號處理中起著十分重要的作用。介紹Matlab環境下設計和實現信號平穩性檢驗系統。該系統主要利用替代數據的平穩性特點，通過在時頻域中分別計算原數據和對應替代數據的平穩度并相互比較，以實現對信號平穩性的檢驗。它可以求出輸入數據的替代數據，并分析原始數據和替代數據的頻域和時頻域性質，同時還可以通過計算原始數據和替代數據各自時頻域的變化程度來判斷原始數據的平穩性。該系統提供了友好的用戶界面。實驗表明，該系統可方便地完成信號的平穩性檢驗，對于測試數據給出了較好的檢驗結果。

關鍵詞:平穩性檢驗;替代數據;時頻分布;隨機信號

中圖分類號:TP274 文獻標識碼:A

文章編號:1004-373X(2010)03-083-04

System of Stationarity Test Based on Matlab

SUI Ye，LI Ming

(School of Information Science Technology，East China Normal University，Shanghai，200241，China)

Abstract:Stationarity test remains a challenge problem in the field of signal processing.Due to the importance of stationa-rity test，to find a solution to that problem is greatly desired.A stationarity test system designed and implemented based on Matlab is introduced.The system realizes stationarity test via comparing characteristics of original data and surrogate in time-frequency distribution which represent their stationarity.And the comparison on the basis of stationarization property of surrogate which is first explored in time-frequency perspective.The system can make surrogates from original data and display their characteristics in both frequency domain and time-frequency domain.The system can also compare original data with its surrogates in time-frequency domain to determine whether it is stationary or non-stationary.The system gives a friendly interface and is convenient to use.The experiment indicates our system can test stationarity of signal well.

0 引 言

信號的平穩性檢驗在隨機信號處理中起著十分基礎的作用。由于平穩信號和非平穩信號的性質差別顯著，因此在處理信號之前先行判斷它的平穩性就顯得尤為重要。雖然信號平穩性的定義十分明確，但是實際判斷過程卻是復雜的，例如觀察尺度對信號平穩性判斷就有很大的影響。

這一領域的研究已經取得了一定的成果。一些人提出了受限和帶參數的非平穩性判定方法，而另一些人則將他們的平穩性判定建立在對原始數據的一些假設上[1]。而對于更一般信號的平穩性檢驗的研究還沒有取得太多成果。文獻[2，3]中又提到了這一問題，并且提出了一種新的檢驗平穩性框架。這一框架混合了時頻透視法和有名的替代數據法[4，5]。它的基本思想是引入“可控噪聲”，即替代數據。并且由于替代數據的一些特性，它可以作為平穩性的評判標準。本文參考了文獻[6]中的平穩性檢驗方法，設計了一個信號平穩性檢驗系統，并在Matlab的GUI 開發環境下實現了圖形用戶界面的設計。實踐表明，本系統不但提供了友好的用戶界面，并且可以方便地完成信號的平穩性檢驗。

1 平穩性檢驗原理

1.1 平穩性定義及其檢驗的重要性

假設有一個高斯過程{xl(t)}(-∞

μx(t)=E[xl(t)](1)

為任意確定t時刻的全體平均。同時:

r(t1，t2)=E[xl(t1)xl(t2)](2)

被稱為自相關函數(ACF)。

對于一個弱平穩過程，它的μx(t)和r(t1，t2)都是時不變的或者說與時間無關的。因此有:

μx(t)=const(3)

r(t1，t2)=E[xl(t+τ)xl(t)]=r(τ)(4)

式中:τ=t1-t2被稱為時延。因此，對于平穩高斯過程{xl(t)}，它的自相關函數或者它的功率譜密度函數(PSD)為:

Sxx(ω)=∫∞-∞rxx(τ)e-jωτdτ(5)

足以確定它的性質[7]。

另一方面，如果{xl(t)}是非平穩的，它的μx(t)和r(t1，t2)就是時變的或者說和時間相關的。這樣它的PSD就應該放在時頻域分析[8]。

由此可見，平穩性檢驗是任何信號處理前必不可少的一步，它決定了后續處理可以使用何種方法。

1.2 替代數據

替代數據的概念最初是由Theiler和其合作作者提出的[4]，這種技術是用來產生一種所謂的“替代數據”，這種替代數據是平穩的，同時保持了原數據的一些相關的統計特性。

Theiler在文獻[4]中提出了一種具體的產生替代數據的方法。由這種方法產生的替代數據是平穩的，同時保持了原數據的二階統計特性。具體地說，替代數據保持了原數據功率譜的幅度值不變。

根據Wiener-Khintchin理論，信號的功率譜等于其傅里葉變換的幅值平方。因此保持信號的功率譜幅度值不變，就是保持其傅里葉變換的幅度值不變。因此，假設原數據為x(t)，它的傅里葉變換為X(f)=∫e-i2πtfx(t)dt。則替代數據s(t)由:



s(t)=∫ei2πtf|X(f)|eiφfdf

(6)

產生。其中，φf是在[-π，π]上均勻分布的隨機相位。這樣就保證了s(t)和x(t)有相同的傅里葉變換幅值。在下面的例子中也可以看到，這樣產生的s(t)也是平穩的。

1.3 時頻分布

時頻分布主要用于分析非平穩隨機信號的功率譜。由于非平穩隨機信號的功率譜是時變的，因此在原來功率譜的基礎上再引入時間軸，成為時頻分布(TFD)。TFD可以顯示出信號的功率譜隨時間的變化情況。

具體來說，根據文獻[9]中的定義，信號x(t)的時頻分布Sx，K(t，f)可以表示為:



Sx，K(t，f)=1K∑Kk=1∫+∞-∞x(s)hk(s-t)e-i2πfsds2(7)



式中:hk(t)是k階Hermite函數，定義式為:

hk(t)=1k!2kπe-t2/2Hk(t)(8)

式中:Hk(t)是k階Hermite多項式。

1.4 平穩性檢驗

平穩性可以體現在頻譜隨時間的波動上。具體來說，對于平穩信號，其頻譜不隨時間變化;而對于非平穩信號，其頻譜會隨時間改變。因此，可以通過比較不同時間點上頻譜的相似程度來判斷信號的平穩性。

按照文獻[4]中的檢驗方法，定義不同時間點上的頻譜與頻譜平均值的距離c(x)n為:



c(x)n=κ(Sx，K(tn，#8226;)，n=1，2，…，N)

(9)

式中:符號“<>”表示求取平均。這里采用文獻[10]中的距離定義:

κ(G，H)=1+∫logG(f)H(f)df#8226;

∫(G~(f)-H~(f))logG~(f)H~(f)df(10)



式中:符號“~”表示對應函數的歸一化函數，例如f~(x)=f(x)/max(f(x))。

距離c(x)n隨時間的波動情況Θ1可以被定義為c(x)n的方差，即:

Θ1=1N∑Nn=1(c(x)n-)2(11)

根據替代數據的定義，它是平穩的。將替代數據的c(x)n隨時間的波動情況記作Θ0。通過比較Θ1和Θ0，可以確定原數據的平穩性。具體來說，將Θ0的概率密度函數記作f(Θ0)，選定適當的門限γ，若f(Θ1)<γ，則判為非平穩隨機信號，反之則判為平穩隨機信號。

2 用戶界面生成

2.1 Matlab中用戶界面的生成

Matlab為用戶設計圖形界面提供了一個高效、方便的集成環境。在Matlab中，基本的圖形對象主要包括坐標軸、控件、下拉菜單和內容菜單。用戶可以通過這些對象設計出界面友好，功能強大，操作簡單的圖形用戶界面。圖形用戶界面的生成主要分為以下幾個步驟:

(1) 規劃所設計的圖形用戶界面，主要包括:確定需要哪些窗口，每個窗口怎樣布局，窗口中的各個對象各有什么功能，對象之間如何配合工作，以及相應的異常處理;

(2) 在Matlab提示行下輸入GUIDE，載入用戶界面開發環境;

(3) 利用Layout Editor，完成用戶面板以及界面的制作，并對相應的按鈕及控件屬性進行適當的設置;

(4) 在Programme Editor中編輯各個對象的回調函數，實現各個對象的具體功能;

(5) 利用Mfile 編譯器生成客戶端，完成隨機數據仿真系統的設計。

2.2 用戶界面介紹

本文所實現的用戶界面主要包括兩個窗口，分別是主窗口和數據生成窗口。由于Matlab對保存繪圖區域有限制，因此設計時沒有在窗口中設置固定的繪圖區域。窗口只相當于一個命令菜單，所有的繪圖將會以獨立窗口的形式根據用戶需求動態產生。這樣便于用戶對比和保存圖片。下面對主要窗口分別加以介紹。

2.2.1 主窗口介紹

主窗口如圖1所示。主窗口主要用于繪制原數據和替代數據的各種波形以及顯示平穩性檢驗結果。

圖1 主窗口示意圖

其中，“Create / Open Original Data”按鈕用來打開數據生成窗口。

“View / Change Parameters”按鈕用來查看或改變當前仿真參數，它在原始數據存在的情況下才有效。主要的仿真參數有:

“Time Scale of TFD”和”Frequency Scale of TFD”用于確定繪制TFD圖片時的時間/頻率軸采樣周期，由于計算和顯示時頻分布圖比較耗費時間，將采樣周期設大，可以提高速度，但是相應的時頻分布圖的分辨率會下降。

“Max Level of Hermite Function”用于確定求TFD時所使用Hermite函數的最高階數。最高階數越高，則分辨率越高，但是相應的計算時間會加長。

“The Number of Surrogates”用于確定平穩性檢驗時所用的參考替代數據個數。個數越多，則檢驗結果越精確，但是會極大地延長計算時間。

“Create Surrogate”按鈕用于產生替代數據，其在原始數據存在的情況下才有效。由于替代數據具有隨機性，因此用戶可以多次產生不同的替代數據，觀察它們的性質。

右上方的下拉菜單用于選擇需要繪圖或者保存數據的對象，主要包括原始數據的時域、頻域和時頻域圖，替代數據的時域、頻域和時頻域圖，以及替代數據的平穩度分布。它在原始數據存在的情況下才有效。

“Show Selected Plot”按鈕用于在新窗口中繪制下拉菜單所選圖線，它在原始數據存在的情況下才有效。

“Save Selected Data”按鈕用于保存下拉菜單所選圖線對應的數據，它在原始數據存在的情況下才有效。

2.2.2 數據生成窗口

數據生成窗口如圖2所示。數據生成窗口主要用來產生實驗用數據或者打開已經存在的實驗數據。

圖2 數據生成窗口示意圖

“Creat Original Data”按鈕用來產生測試用數據。按下此按鈕后會提示輸入產生數據用的參數。由于數據是通過公式:

x(t)=sin[P1sin(P2t)t](12)

產生的調頻信號，因此需要確定參數P1和P2，另外還要確定t的區間和采樣周期。數據成功產生后會在新建窗口中顯示該數據時域波形。如果當前存在數據波形，將會覆蓋它。

“Open Original Data”按鈕用來打開已經存在的數據文件。選擇好文件后會提示輸入參數。主要包括讀入數據的時間起點、時間采樣周期和數據長度。數據成功讀入后會在新建窗口中顯示該數據時域波形。如果當前存在數據波形，將會覆蓋它。

“Confirm”按鈕用于確認新建窗口顯示的數據就是用戶想要的數據，并返回主窗口。它在創建或打開的數據存在的情況下才有效。

3 數據仿真和分析

將實驗數據取為調頻信號x(t)=sin(sin(t/8)πt)。t起始為0，采樣周期為0.1 s，數據長度為400個點，其時域波形和頻域波形如圖3所示。

圖3 原始數據時域和頻域示意圖

由式(6)產生的替代數據s(t)的時域波形和頻域波形如圖4所示。

圖4 替代數據時域和頻域示意圖

由圖3和圖4不難看出，替代數據與原數據的傅里葉變換幅值相同，但替代數據傅里葉變換的相位是隨機的。

圖5顯示了由式(7)計算得到的原數據和替代數據的時頻分布圖。由圖5中可見，原數據的時頻分布圖有明顯的結構性。它表明是非平穩的，而替代數據的時頻分布圖的結構性較原數據有明顯減弱，表明替代數據的平穩性增加。

圖5 原數據與替代數據時頻分布示意圖

由式(11)計算得到的Θ0的概率密度函數f(Θ0)如圖6所示，其中一共計算了1 000次替代數據。

由圖6可見，替代數據的平穩度主要分布在0.02附近。數據的平穩度落在0~0.04之間可以認為是平穩的，而在此之外可以認為是非平穩的。

由式(11)計算得到的Θ1=0.046。位于上述區間之外，因此判為非平穩。這一結果也與圖5所示的結果相吻合。

圖6 替代數據平穩度分布示意圖

4 結 語

利用替代數據法和時頻透視法，并采用Matlab的GUI開發環境，設計了一個信號平穩性檢驗系統。該系統的用戶界面友好。利用該系統可以觀察信號及其替代數據的頻域和時頻域波形，檢驗信號的平穩性。

參考文獻

[1]Hobijin B，Franses P H，Ooms M.Generalization of the KPSS-test for Stationarity[J].Statitica Neerlandica，2004，58(4):483-502.

[2]Xiao J，Borgnat P，Flandrin P.Testing Stationarity with Time-frequency Surrogates[A].Proc.EUSIPCO-07[C].2007.

[3]Xiao J，Borgnat P，Flandrin P，et al.Testing Stationarity with Surrogates a One-class SVM Approach[A].Proc.IEEE Stat.Sig.Proc.Workshop SSP-07[C].Madison，2007:720-724.

[4]Theiler J，Eubank S，Longtin A，et al.Testing for Nonlinea-rity in Time Series:The Method of Surrogate Data[J].Physica D，1992:58-77.

[5]Schreiber T，Schmitz A.Surrogate Time Series[J].Physica D，2000，142:346-382.

[6]Pierre Borgnat，Patrick Flandrin.Stationarization via Surrogates[J].Journal of Statistical Mechanics，2009(1):1-14.

[7]Papoulis A.Probability Random Variables and Stochastic Process[M].New York:McGraw-Hill，1984.

[8]Priestley M B.Evolutionary Spectra and Non-Stationary Process[J].Roy.Stat.Soc.Series B，1965，27(2):204-237.

[9]Bayram M，Baraniuk R G.Multiple Window Time-varying Spectrum Estimation[A].Nonlinear and Nonstationary Signal Processing ed W J Fitzgerald et al[C].Cambridge:Cambridge University Press，2000:292-316.

[10]Proestley M B，Rao T S.A Test of Non-stationarity of Time-series[J].Royal Statictic Society，1969:140-149.