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

關(guān)鍵詞:平穩(wěn)性檢驗;替代數(shù)據(jù);時頻分布;隨機(jī)信號

中圖分類號:TP274 文獻(xiàn)標(biāo)識碼: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 引 言

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

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

1 平穩(wěn)性檢驗原理

1.1 平穩(wěn)性定義及其檢驗的重要性

假設(shè)有一個高斯過程{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)

被稱為自相關(guān)函數(shù)(ACF)。

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

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

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

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

Sxx(ω)=∫∞-∞r(nóng)xx(τ)e-jωτdτ(5)

足以確定它的性質(zhì)[7]。

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

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

1.2 替代數(shù)據(jù)

替代數(shù)據(jù)的概念最初是由Theiler和其合作作者提出的[4]，這種技術(shù)是用來產(chǎn)生一種所謂的“替代數(shù)據(jù)”，這種替代數(shù)據(jù)是平穩(wěn)的，同時保持了原數(shù)據(jù)的一些相關(guān)的統(tǒng)計特性。

Theiler在文獻(xiàn)[4]中提出了一種具體的產(chǎn)生替代數(shù)據(jù)的方法。由這種方法產(chǎn)生的替代數(shù)據(jù)是平穩(wěn)的，同時保持了原數(shù)據(jù)的二階統(tǒng)計特性。具體地說，替代數(shù)據(jù)保持了原數(shù)據(jù)功率譜的幅度值不變。

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



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

(6)

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

1.3 時頻分布

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

具體來說，根據(jù)文獻(xiàn)[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函數(shù)，定義式為:

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

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

1.4 平穩(wěn)性檢驗

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

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



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

(9)

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

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

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



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

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

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

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

2 用戶界面生成

2.1 Matlab中用戶界面的生成

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

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

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

(3) 利用Layout Editor，完成用戶面板以及界面的制作，并對相應(yīng)的按鈕及控件屬性進(jìn)行適當(dāng)?shù)脑O(shè)置;

(4) 在Programme Editor中編輯各個對象的回調(diào)函數(shù)，實現(xiàn)各個對象的具體功能;

(5) 利用Mfile 編譯器生成客戶端，完成隨機(jī)數(shù)據(jù)仿真系統(tǒng)的設(shè)計。

2.2 用戶界面介紹

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

2.2.1 主窗口介紹

主窗口如圖1所示。主窗口主要用于繪制原數(shù)據(jù)和替代數(shù)據(jù)的各種波形以及顯示平穩(wěn)性檢驗結(jié)果。

圖1 主窗口示意圖

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

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

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

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

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

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

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

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

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

2.2.2 數(shù)據(jù)生成窗口

數(shù)據(jù)生成窗口如圖2所示。數(shù)據(jù)生成窗口主要用來產(chǎn)生實驗用數(shù)據(jù)或者打開已經(jīng)存在的實驗數(shù)據(jù)。

圖2 數(shù)據(jù)生成窗口示意圖

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

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

產(chǎn)生的調(diào)頻信號，因此需要確定參數(shù)P1和P2，另外還要確定t的區(qū)間和采樣周期。數(shù)據(jù)成功產(chǎn)生后會在新建窗口中顯示該數(shù)據(jù)時域波形。如果當(dāng)前存在數(shù)據(jù)波形，將會覆蓋它。

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

“Confirm”按鈕用于確認(rèn)新建窗口顯示的數(shù)據(jù)就是用戶想要的數(shù)據(jù)，并返回主窗口。它在創(chuàng)建或打開的數(shù)據(jù)存在的情況下才有效。

3 數(shù)據(jù)仿真和分析

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

圖3 原始數(shù)據(jù)時域和頻域示意圖

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

圖4 替代數(shù)據(jù)時域和頻域示意圖

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

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

圖5 原數(shù)據(jù)與替代數(shù)據(jù)時頻分布示意圖

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

由圖6可見，替代數(shù)據(jù)的平穩(wěn)度主要分布在0.02附近。數(shù)據(jù)的平穩(wěn)度落在0~0.04之間可以認(rèn)為是平穩(wěn)的，而在此之外可以認(rèn)為是非平穩(wěn)的。

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

圖6 替代數(shù)據(jù)平穩(wěn)度分布示意圖

4 結(jié) 語

利用替代數(shù)據(jù)法和時頻透視法，并采用Matlab的GUI開發(fā)環(huán)境，設(shè)計了一個信號平穩(wěn)性檢驗系統(tǒng)。該系統(tǒng)的用戶界面友好。利用該系統(tǒng)可以觀察信號及其替代數(shù)據(jù)的頻域和時頻域波形，檢驗信號的平穩(wěn)性。

參考文獻(xiàn)

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

[2]Xiao J，Borgnat P，F(xiàn)landrin P.Testing Stationarity with Time-frequency Surrogates[A].Proc.EUSIPCO-07[C].2007.

[3]Xiao J，Borgnat P，F(xiàn)landrin 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.