1、第7章 隨機(jī)信號處理與譜分析信號從廣義上可分為兩大類：確定性信號和隨機(jī)信號，上一章內(nèi)容所針對的信號形式就為確定性信號，本章將學(xué)習(xí)隨機(jī)信號的處理方法，其中譜分析是一個非常重要的部分。這里介紹的譜分析包括功率譜估計及高階譜估計兩大部分，各自又有很多種估計方法。由于隨機(jī)信號的處理方法是基于統(tǒng)計特性的分析，因此本章將首先給出描述統(tǒng)計特性的一些特征量，之后介紹功率譜估計方法，包括經(jīng)典譜估計法、改進(jìn)的直接法、AR模型法等，接下來將介紹高階譜估計方法，包括非參數(shù)化方法及參數(shù)模型法，最后給出一個綜合實例。學(xué)習(xí)目標(biāo)r 掌握隨機(jī)信號的統(tǒng)計描述和統(tǒng)計特性r 掌握經(jīng)典功率譜估計方法及改進(jìn)方法r 了解AR模型功率譜估計

2、方法r 理解高階累積量和高階譜定義r 了解高階譜估計方法練習(xí)案例r 估計給定隨機(jī)序列的均值和方差。r 產(chǎn)生一個相關(guān)正態(tài)隨機(jī)序列，并計算它的自相關(guān)函數(shù)。r 產(chǎn)生兩個被白噪聲污染的隨機(jī)序列，并分別計算二者的自協(xié)方差及互協(xié)方差。r 產(chǎn)生一頻率為300Hz的諧波信號（分實數(shù)和復(fù)數(shù)兩種情況），疊加高斯白噪聲，利用周期圖法估計該混合信號的功率譜。r 產(chǎn)生一頻率為300Hz的諧波信號（復(fù)數(shù)形式），疊加高斯白噪聲，利用Welch法估計該混合信號的功率譜。r 模擬一個非高斯ARMA隨機(jī)過程，并估計它的3階累積量，采用函數(shù)cum3x。r 分別利用直接法和間接法估計一個已知的非高斯信號的雙譜。r 估計一已知的非高斯

3、信號ARMA模型參數(shù)，并利用此模型估計該信號的雙譜。r 生成一混合信號，包含頻率分別為300Hz和320Hz兩個諧波信號，并疊加高斯白噪聲，分別利用周期圖法、Welch法、AR模型法、Burg法、特征向量法以及MUSIC法估計該混合信號的功率譜，并對結(jié)果進(jìn)行分析比較。7.1 概述信號常用一組變量值表示，如，如果序列在每個時刻的取值不是隨機(jī)的，而是服從某種固定函數(shù)的關(guān)系，則稱之為確定性信號，比如上一章介紹的階躍信號、方波信號、正弦波信號等。與確定性信號不同，如果序列在每個時刻的取值是隨機(jī)變量，則稱之為隨機(jī)信號。隨機(jī)信號也稱為隨機(jī)過程，具有以下特點：r 隨機(jī)信號在任何時間的取值都不能先驗確定，是隨

4、機(jī)變化的；r 雖然隨機(jī)信號取值不能先驗確定，但這些取值服從某種統(tǒng)計規(guī)律，或者說隨機(jī)信號可以用概率分布特性統(tǒng)計描述。例如，一個正弦波信號如果式中相位是一個隨機(jī)變量，假設(shè)在-，上均勻分布那么該正弦波信號就是一個隨機(jī)信號。隨機(jī)信號又分為平穩(wěn)隨機(jī)信號和非平穩(wěn)隨機(jī)信號兩大類，如果隨機(jī)信號的概率分布不隨時間推移而變化，即信號的統(tǒng)計特性與起始時間無關(guān)，只與時間間隔有關(guān)，這類隨機(jī)信號就稱為平穩(wěn)隨機(jī)信號，否則為非平穩(wěn)隨機(jī)信號。因而，平穩(wěn)信號常被稱作時不變信號，而非平穩(wěn)信號常被稱作時變信號。提示：時變與時不變是指隨機(jī)信號的某個或某幾個統(tǒng)計量是否隨時間變化，而非針對信號的波形而言。隨機(jī)信號的一個重要性質(zhì)是遍歷性，它

5、所關(guān)心的是從隨機(jī)信號的一次觀測值能否估計信號的統(tǒng)計量。如果一個隨機(jī)信號的均值和方差都具有遍歷性，則稱該信號具有各態(tài)歷經(jīng)性，此時只通過一次觀測值就可估計隨機(jī)過程的多個統(tǒng)計特性，這在實際當(dāng)中就非常有意義了。為了簡化分析和運(yùn)算，我們常常假設(shè)所獲得的隨機(jī)信號具有各態(tài)歷經(jīng)性。7.2 隨機(jī)信號的統(tǒng)計特性隨機(jī)信號可用其統(tǒng)計特性進(jìn)行描述，這些統(tǒng)計特性又可進(jìn)一步劃分為一階、二階和高階（三階及三階以上）統(tǒng)計特性。一階特性主要是指均值函數(shù)，即信號的數(shù)學(xué)期望，二階特性主要包括方差、相關(guān)函數(shù)、協(xié)方差函數(shù)和功率譜密度，高階特性則指代高階矩。對于平穩(wěn)信號而言，二階特性往往是最重要的。接下來將對隨機(jī)信號的幾個主要統(tǒng)計特性逐一

6、介紹，并給出MATLAB工具箱提供的相應(yīng)函數(shù)。7.2.1 均值和方差均值和方差是常用的描述隨機(jī)信號統(tǒng)計特性的特征量。假設(shè)表示一個各態(tài)歷經(jīng)的隨機(jī)信號，那么該信號的均值和方差的定義式分別為：式中表示信號的概率密度函數(shù)。在MATLAB工具箱中，分別提供了兩個函數(shù)來實現(xiàn)隨機(jī)序列的均值和方差的計算，即均值函數(shù)mean和方差函數(shù)var。l m = mean(X)該函數(shù)計算均值的公式為：格式中若輸入?yún)⒘縓是向量，則返回值m為該向量的均值；若輸入?yún)⒘縓為矩陣，則返回值m為一個行向量，其元素分別為矩陣X每一列的均值。l y = var(X)若輸入?yún)⒘縓是向量，則返回值y為該向量的方差；若輸入?yún)⒘縓為矩陣，則返回

7、值y為一個行向量，其元素分別為矩陣X每一列的均值。例7-1估計給定隨機(jī)序列的均值和方差。% 隨機(jī)序列的均值和方差估計x = normrnd(0.1,0.8,1,6) %產(chǎn)生正態(tài)分布的隨機(jī)序列xm = mean(x)xv = var(x)這里x是一個長度為6的正態(tài)分布隨機(jī)序列，結(jié)果如下：>>x = -0.2461 -1.2325 0.2003 0.3301 -0.8172 1.0527xm = -0.1188xv = 0.68407.2.2 相關(guān)函數(shù)相關(guān)函數(shù)是最為常用的描述平穩(wěn)隨機(jī)信號統(tǒng)計特性的二階統(tǒng)計量之一，可用來描述隨機(jī)信號不同時刻間的相關(guān)程度。我們首先定義單個平穩(wěn)隨機(jī)信號的自相

8、關(guān)函數(shù)，其定義式為：同理，對于兩個平穩(wěn)隨機(jī)信號和，其互相關(guān)函數(shù)可定義為：式中是和 的二維混合概率密度函數(shù)。MATLAB工具箱提供了函數(shù)xcorr來實現(xiàn)相關(guān)函數(shù)的估計，其調(diào)用格式為：l c = xcorr(x,y);返回隨機(jī)序列x和y的互相關(guān)函數(shù)，長為2N-1，假設(shè)x和y長度都為N。l c = xcorr(x,y,maxlags);輸入?yún)⒘縨axlags為x和y之間的最大延遲，返回值c長度為2maxlags+1，如缺省，c長度默認(rèn)為2N-1。l c = xcorr(x,y,maxlags,option);這里option包括以下幾個選擇項：r biased有偏估計，自相關(guān)函數(shù)和互相關(guān)函數(shù)的估計表

9、達(dá)式分別取：r unbiased無偏估計，自相關(guān)函數(shù)和互相關(guān)函數(shù)的估計表達(dá)式分別取：r coeff m=0的相關(guān)函數(shù)值歸一化為1。r none 不做歸一化處理。作為特例，函數(shù)xcorr也可用來計算自相關(guān)函數(shù)，調(diào)用格式為：l c = xcorr(x);l c = xcorr(x,maxlags,option);例7-2產(chǎn)生一個相關(guān)正態(tài)隨機(jī)序列，并計算它的自相關(guān)函數(shù)。% 相關(guān)函數(shù)估計程序a = 0.85;sigma = 1.5;N=200;x = randn(N,1);y(1) = sigma*x(1)/sqrt(1-a2);for i=2:Ny(i) = a*y(i-1)+sigma*x(i);

10、endc = xcorr(y,'coeff');subplot(121);plot(y,'LineWidth',2);grid onsubplot(122);plot(c,'LineWidth',2);grid on程序中的變量y就是產(chǎn)生的相關(guān)正態(tài)隨機(jī)序列，其相關(guān)函數(shù)表達(dá)式為程序運(yùn)行結(jié)果如圖7-1所示。圖7-1 隨機(jī)序列波形及其自相關(guān)函數(shù)7.2.3 協(xié)方差函數(shù)除了相關(guān)函數(shù)，協(xié)方差也是一個重要的描述隨機(jī)過程統(tǒng)計特性的二階統(tǒng)計量。同樣，我們也首先定義單個平穩(wěn)隨機(jī)信號的自協(xié)方差函數(shù)，其定義式為：自協(xié)方差函數(shù)用于表示它的兩個隨機(jī)變量和之間的相關(guān)程度。同理

11、，定義兩個平穩(wěn)隨機(jī)信號互協(xié)方差函數(shù)，其估計表達(dá)式為：互協(xié)方差函數(shù)表示隨機(jī)變量和之間的相關(guān)程度，實際上，也是表示隨機(jī)信號和之間的相關(guān)程度。MATLAB工具箱提供了函數(shù)cov計算隨機(jī)信號的協(xié)方差函數(shù)，其調(diào)用格式為：l C = cov(X);當(dāng)X為向量時，返回值C是一個標(biāo)量值，表示方差；當(dāng)X為矩陣形式，每一行表示一次觀測值，每列對應(yīng)一個變量，則返回值C是一個協(xié)方差矩陣。提示：方差函數(shù)等價于協(xié)方差函數(shù)的對角元素，即var(X) = diag(cov(X)。l C = cov(X,Y);計算輸入向量X和Y的互協(xié)方差矩陣，X和Y是具有相同長度的列向量。例7-3產(chǎn)生兩個被白噪聲污染的隨機(jī)序列，并分別計算二者

12、的自協(xié)方差及互協(xié)方差。% 計算自協(xié)方差及互協(xié)方差程序fs = 100;t =0:1/fs:.1;f1 = 15;f2 = 30;x0 = cos(2*pi*f1*t).'x = x0+randn(size(x0);y0 = cos(2*pi*f2*t).'y = y0+randn(size(y0);Cx = cov(x)Cy = cov(y)Cxy = cov(x,y)程序運(yùn)行結(jié)果如下：>> Cx = 1.5176Cy = 0.7867Cxy = 1.5176 0.5067 0.5067 0.7867同時，我們計算序列x和y的方差，得到：>> var(x

13、)ans = 1.5176>> var(y)ans = 0.7867可以看出，方差函數(shù)確實等價于協(xié)方差函數(shù)的對角元素。7.3 功率譜估計 功率譜估計是利用已觀測到的一定數(shù)量樣本數(shù)據(jù)估計一個平穩(wěn)隨機(jī)信號的功率譜密度，因其能夠分析信號的能量隨頻率變化的分布特性，在許多實際應(yīng)用中功率譜的分析與估計已變得越來越重要。例如，強(qiáng)噪聲背景下弱信號的檢測、雷達(dá)和聲納回波信號分析、故障診斷等多個領(lǐng)域都涉及到了功率譜的分析與估計。功率譜密度，簡稱功率譜，可定義為從表達(dá)式中可以看出，隨機(jī)信號的自相關(guān)函數(shù)與它的功率譜構(gòu)成一對離散Fourier變換對，二者分別從時域和頻域反映了隨機(jī)信號的二階統(tǒng)計特性。功率譜

14、估計方法有很多種，大體上可以分為三大類：r 非參數(shù)方法r 參數(shù)方法r 子空間方法非參數(shù)方法是直接基于觀測信號本身進(jìn)行功率譜估計的一類方法，而不涉及一些模型參數(shù)等。其中最簡單的一種方法是周期圖法，在其基礎(chǔ)上的改進(jìn)方法包括Welch方法等。參數(shù)方法是一類基于模型參數(shù)的方法，該類方法假設(shè)觀測信號為一個帶有加性白噪聲的線性系統(tǒng)的輸出，在此架設(shè)基礎(chǔ)上進(jìn)行功率譜的估計。代表算法包括Yule-Walker自回歸（AR）算法和Burg算法。這類算法首先要估計產(chǎn)生觀測信號的線性系統(tǒng)的一組參數(shù)（系數(shù)）。當(dāng)樣本數(shù)較少時，參數(shù)方法往往優(yōu)于經(jīng)典的非參數(shù)方法。子空間方法常常被稱作高分辨率方法或超分辨率方法，這類方法基于觀

15、測信號相關(guān)矩陣的特征分析或特征分解來對功率譜進(jìn)行估計，例如多信號分類（MUSIC）方法和特征向量（EV）法。這類方法非常適用于分析信號中含有線譜分量的情況，尤其在信噪比很低的情況下，子空間方法就顯示出了優(yōu)越性。MATLAB信號處理工具箱針對以上三類功率譜估計方法分別提供了相應(yīng)的函數(shù)，線將這些方法匯總在一個表中，見表7-1。表7-1 功率譜估計方法及MATLAB工具箱函數(shù)類別方法函數(shù)非參數(shù)方法周期圖法periodogram,Welch方法pwelch參數(shù)方法AR模型法spectrum.yulear, pyulearBurg方法spectrum.burg, pburg子空間法MUSIC方法spec

16、trum.music, pmusic特型向量法spectrum.eigenvector, peig7.3.1 非參數(shù)方法7.3.1.1 周期圖法非參數(shù)方法中最簡單的是周期圖法，它是把隨機(jī)序列x(n)的N個觀測數(shù)據(jù)視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，然后再取其幅值的平方，并除以N，作為序列x(n)真實功率譜的估計，即這里周期圖法的優(yōu)點是計算效率高，但缺點是頻率分辨率較低。MATLAB信號處理工具箱提供了函數(shù)periodogram來實現(xiàn)周期圖法的功率譜估計，其調(diào)用格式包括：l Pxx,w = periodogram(x,window,nfft,range);l Pxx,f =

17、periodogram(x,window,nfft,fs,range);l periodogram();返回值Pxx為利用周期圖法估計的功率譜；w和f分別是單位為弧度和Hz的頻率向量，與Pxx相同長度；輸入?yún)⒘縳為有限長信號向量；window為所選的窗函數(shù)，默認(rèn)矩形窗；nfft是FFT點數(shù)，一般取2的整數(shù)次冪，默認(rèn)值256；range包括兩種情況：r twosided在頻率區(qū)間0, fs)上計算雙邊PSD，當(dāng)樣本數(shù)據(jù)x為復(fù)數(shù)時，這是默認(rèn)情況；r onesided在指定的頻率區(qū)間上計算單邊PSD，當(dāng)樣本數(shù)據(jù)x為實數(shù)時，這是默認(rèn)情況。提示：range字段可以放在輸入?yún)⒘縲indow之后的任意位置。

18、例7-4產(chǎn)生一頻率為300Hz的諧波信號（分實數(shù)和復(fù)數(shù)兩種情況），疊加高斯白噪聲，利用周期圖法估計該混合信號的功率譜。% 周期圖法估計信號功率譜Fs = 1000; t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); % 產(chǎn)生實數(shù)混合信號periodogram(x,'twosided',512,Fs) % 使用默認(rèn)窗函數(shù)y = exp(j*2*pi*t*300)+0.1*randn(size(t); % 產(chǎn)生復(fù)數(shù)混合信號periodogram(y,'onesided',512,Fs) 程序運(yùn)行結(jié)果如圖7-2所

19、示。可以看出，當(dāng)信號形式為實數(shù)時，功率譜圖形會產(chǎn)生關(guān)于奈奎斯特頻率對稱的兩個波峰，而當(dāng)信號形式為復(fù)數(shù)時，僅保留300Hz位置的波峰。圖7-2 周期圖法估計功率譜密度圖（左圖為實數(shù)信號，右圖為復(fù)數(shù)信號）7.3.1.2 Welch方法周期圖法估計出的功率譜性能并不好，當(dāng)數(shù)據(jù)長度過大時，譜的曲線起伏劇烈，而當(dāng)數(shù)據(jù)長度過小時，譜的分辨率又不好，因此在其基礎(chǔ)上又產(chǎn)生了一系列的改進(jìn)方法，共同的原則是將周期圖法進(jìn)行平滑，使得估計方差減小，其中以Welch方法為代表。Welch方法是一種最為廣泛的經(jīng)典功率譜估計方法，該法在周期圖法的基礎(chǔ)上主要進(jìn)行了兩處改進(jìn)：其一是將采樣數(shù)據(jù)進(jìn)行分段，并允許各段數(shù)據(jù)有交疊的部分

20、；其二是每段數(shù)據(jù)加窗處理時，窗函數(shù)不局限于矩形窗，也可使用漢寧窗或海明窗等，從而可以改善譜失真的現(xiàn)象，具體步驟如下： 步驟r 先將N個數(shù)據(jù)分為L段，每段M個數(shù)據(jù)；r 選擇適當(dāng)?shù)拇昂瘮?shù)w(n)，對每段數(shù)據(jù)依次加權(quán)，然后確定每段的周期圖；r 對分段周期圖進(jìn)行平均，得到功率譜估計。根據(jù)以上步驟，最終得到的功率譜估計表達(dá)式為式中為第段的采樣數(shù)據(jù)，為選取的窗函數(shù)，是每段數(shù)據(jù)的長度，是歸一化因子，用于保證所得到的譜是漸進(jìn)無偏估計，表達(dá)式為MATLAB信號處理工具箱提供了函數(shù)pwelch來實現(xiàn)Welch改進(jìn)周期圖法的功率譜估計，其調(diào)用格式為：l Pxx,w = pwelch (x,window,noverl

21、ap,nfft);l Pxx,f = pwelch (x,window,noverlap,nfft,fs);l Pxx,f = pwelch (x,window,noverlap,f,fs);l pwelch (x, );格式中，noverlap為分段重疊點數(shù)，必須小于窗長度，默認(rèn)值為窗長的一半，其他各參量同periodogram。例7-5產(chǎn)生一頻率為300Hz的諧波信號（復(fù)數(shù)形式），疊加高斯白噪聲，利用Welch法估計該混合信號的功率譜。% Welch法估計信號功率譜Fs = 1000; t = 0:1/Fs:.3;x = exp(j*2*pi*t*300)+0.1*randn(size(t

22、); pwelch(x,128,256,Fs) 程序運(yùn)行結(jié)果如圖7-3所示，與圖7-2相比較可以看出，Welch方法估計得到的功率譜圖比周期圖法估計結(jié)果要平滑很多，曲線起伏較為平緩，優(yōu)于周期圖法。圖7-3 Welch法估計信號功率譜圖7.3.2 參數(shù)方法任何具有有理功率譜密度的隨機(jī)信號都可以看作一白噪聲激勵一個線性系統(tǒng)而形成，如果根據(jù)已觀測到得數(shù)據(jù)估計出這個線性系統(tǒng)地模型參數(shù)，就不必假設(shè)N個觀測值之外的數(shù)據(jù)都為零，從而克服經(jīng)典功率譜估計的缺點。因而，參數(shù)方法并不直接從觀測信號來估計功率譜，而是首先建立一個線性系統(tǒng)模型，將觀測信號視為該系統(tǒng)的輸出，而后對線性系統(tǒng)的參數(shù)進(jìn)行估計，具體步驟如下： 步

23、驟r 選擇合適模型；r 用已觀測到的數(shù)據(jù)估計模型參數(shù)；r 將模型參數(shù)代入相應(yīng)公式，得到功率譜估計值。這類方法的優(yōu)點之一是，當(dāng)信號采樣長度較短時，參數(shù)方法能夠獲得比非參數(shù)方法更高的功率譜分辨率。這里我們介紹兩種較為典型的參數(shù)方法，即AR模型方法和Burg法。7.3.2.1 AR模型法AR模型（自回歸模型）使用最為廣泛，這是因其參數(shù)計算為線性方程，比較簡單，同時適用于表示很窄的頻譜。AR模型是一個全極點的模型，該模型當(dāng)前輸出等于當(dāng)前輸入與過去輸出的加權(quán)和，利用差分方程可表示為：式中是高斯白噪聲序列，為AR模型階數(shù)，為AR模型系數(shù)，。由上面的差分方程，進(jìn)一步得到系統(tǒng)函數(shù)最終得到AR模型的功率譜估計式

24、：待估計參數(shù)包括，共個，這些參數(shù)可由Yule-Walker方程求得提示：任何ARMA和AR模型均可用一個無限階的AR模型表示，即使模型選擇不當(dāng)，只要模型階數(shù)足夠高，仍能比較好地逼近建模的隨機(jī)過程。MATLAB信號處理工具箱提供了函數(shù)aryule來估計AR模型參數(shù)，其調(diào)用格式為：l a,e = aryule(x,p)返回值a包含AR模型系數(shù)參量，返回值e為方差參數(shù)，x為輸入信號向量，p為AR模型階數(shù)。這里該函數(shù)采用的是Yule-Walker方法估計模型的參數(shù)。例7-6估計例7-5中信號的10階AR模型參量。% AR模型參數(shù)估計Fs = 1000;t = 0:1/Fs:.3;x = cos(2*p

25、i*t*300)+0.1*randn(size(t); a,e = aryule(x,10)程序運(yùn)行結(jié)果如下：>> a = 1.0000 0.1433 0.4348 -0.3201 -0.1534 0.1377 -0.1429 -0.0693 0.0136 0.0276 -0.0281e = 0.0218同時，MATLAB工具箱提供了函數(shù)pyulear利用AR模型估計功率譜，即首先由Yule-Walker方法估計AR模型參數(shù)，而后將參數(shù)代入AR模型的功率譜估計公式，便得到觀測信號的功率譜估計值。函數(shù)pyulear的具體調(diào)用格式包括：l Pxx = pyulear(x,p,nfft,

26、fs)l Pxx,f = pyulear(x,p,nfft,fs,'range')l Pxx,w = pyulear(x,p,nfft,'range')l pyulear(.)格式中，輸入?yún)⒘縫為AR模型階數(shù)，其他參量同函數(shù)periodogram例7-7采用AR模型法估計例7-5中信號的功率譜密度。% AR模型法估計功率譜Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); pyulear(x,10);程序運(yùn)行結(jié)果如圖7-4所示，注意到，圖形中的橫坐標(biāo)為歸一化之后的頻率，這與上面幾個函數(shù)如p

27、eriodogram、pwelch是不同的。圖7-4 AR模型法估計功率譜圖7.3.2.2 Burg法Burg法估計功率譜又被稱為最大熵譜估計，這里首先定義了譜熵的概念，其表達(dá)式為：式中為功率譜密度函數(shù)。Burg最大熵譜估計可以表述為，求功率譜，使其在約束條件下之下，能夠使譜熵達(dá)到最大。對這一具有約束條件的優(yōu)化問題采用Lagrange乘子法求解，得到式中為Lagrange乘子。MATLAB信號處理工具箱提供了函數(shù)pburg實現(xiàn)功率譜估計，其調(diào)用格式為：l Pxx = pburg (x,p,nfft,fs)l Pxx,f = pburg (x,p,nfft,fs,'range')

28、l Pxx,w = pburg (x,p,nfft,'range')l pburg (.)格式中，Pxx為功率譜估計值，x是輸入信號向量，p為AR模型階數(shù)，nfft是FFT點數(shù)，fs為采樣頻率，range同函數(shù)periodogram。例7-8采用Burg法估計例7-5中信號的功率譜密度。% Burg法估計功率譜Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); pburg(x,10);程序運(yùn)行結(jié)果如圖7-5所示，注意到，圖形中的橫坐標(biāo)也為歸一化之后的頻率，這與函數(shù)pyulear相同。圖7-5 Burg法估

29、計功率譜圖7.3.3 子空間法子空間法主要包括MUSIC法和特征向量法，這類方法都是以自相關(guān)矩陣的特征分析為基礎(chǔ)，將系統(tǒng)自相關(guān)矩陣分成兩個子空間：信號子空間和噪聲子空間。其中特征向量法主要適用于混有高斯白噪聲的正弦信號的功率譜估計，而MUSIC法適合更為普遍情況下正弦信號的功率譜估計。子空間法功率譜估計步驟一般如下： 步驟r 根據(jù)N個觀測樣本值，估計自相關(guān)函數(shù)，假設(shè)和分別為第個復(fù)正弦信號的功率和頻率，是白噪聲方差；r 對進(jìn)行特征分解，得到個特征值，對應(yīng)特征向量；r 估計功率譜7.3.3.1 特征向量法在上面的功率譜估計表達(dá)式中，如果令，則所得的功率譜估計就稱為特征向量估計，此時表達(dá)式寫作在MA

30、TLAB信號處理工具箱中提供了函數(shù)peig來實現(xiàn)特征向量法估計功率譜，其調(diào)用格式為：l S,w = peig(x,p);l S,f = peig(x,p,f,fs)l S,f = peig(x,p,nfft,fs,nwin,noverlap);l peig()格式中，返回值S為偽譜估計值；w和f分別是單位為弧度和Hz的返回頻率值；x為輸入信號向量；p用來指定信號子空間中特征向量的數(shù)目；nfft為FFT算法的長度，默認(rèn)值為256；nwin用于指定矩形窗的寬度，信號向量x被分割成長度為nwin的各段，默認(rèn)值為2*p；noverlap是每段之間重疊的長度，默認(rèn)值為nwin-1。例7-9試用特征向量法

31、估計例7-5中信號的功率譜密度。% 特征向量法估計功率譜Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+cos(2*pi*t*250)+0.1*randn(size(t); p=10;peig(x,p)程序運(yùn)行結(jié)果如圖7-6所示。7.3.3.2 MUSIC法在上面的功率譜估計表達(dá)式中，如果令，則所得的功率譜估計就稱為MUSIC法估計，此時表達(dá)式寫作在MATLAB信號處理工具箱中提供了函數(shù)pmusic來實現(xiàn)MUSIC法功率譜估計，其調(diào)用格式為：l S,w = pmusic(x,p);l S,f = pmusic(x,p,f,fs)l S,f = pmusi

32、c(x,p,nfft,fs,nwin,noverlap);l pmusic()格式中各參數(shù)說明同函數(shù)peig。例7-10試用MUSIC法估計例7-5中信號的功率譜密度。% MUSIC法估計功率譜Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+cos(2*pi*t*250)+0.1*randn(size(t); p=10;pmusic(x,p)程序運(yùn)行結(jié)果如圖7-7所示。 圖7-6 特征向量法估計功率譜圖 圖7-7 MUSIC法估計功率譜圖7.4 高階譜估計前面介紹的功率譜估計過程中，通常需要假設(shè)觀測到的數(shù)據(jù)是由高斯白噪聲激勵線性系統(tǒng)產(chǎn)生的，而實際當(dāng)中信號

33、或噪聲往往并不服從高斯分布，并且信號的自相關(guān)函數(shù)只包含幅度信息，而不包含相位信息，這樣就使得二階統(tǒng)計量（時域的相關(guān)函數(shù)，頻域的功率譜）存在一定的局限性，難以滿足實際需要。因而研究三階及更高階的統(tǒng)計量是非常必要的，將它們統(tǒng)稱為高階統(tǒng)計量。目前，高階統(tǒng)計量已成為信號處理的一種有力工具，并且在信號分析與處理領(lǐng)域中扮演著越來越重要的角色。本節(jié)將首先給出一些相關(guān)概念，如高階累積量和高階譜，之后具體介紹高階譜估計的方法，包括非參數(shù)化的估計方法以及參數(shù)模型法。7.4.1 高階累積量和高階譜定義高階累積量和高階譜是兩個最為常用的高階統(tǒng)計量，這里首先給出二者的定義。為了方便表示，假設(shè)隨機(jī)過程的均值為0，那么零均

34、值平穩(wěn)隨機(jī)過程的二階、三階、四階累積量定義式分別為：這里假設(shè)為實信號。根據(jù)以上定義式，我們知道二階累積量其實就是信號的自協(xié)方差函數(shù)，當(dāng)均值為零時，也等價于信號的自相關(guān)函數(shù)。這里我們僅討論到第四階累積量，因為實際應(yīng)用中三階和四階累積量是最為常用的，而再高階的累積量應(yīng)用場合較少，這里就不再討論。對上面二階、三階、四階累積量做Fourier變換，有上面三個表達(dá)式分別是功率譜、雙譜和三譜的定義式。所謂雙譜是具有兩個頻率變量的函數(shù)，而三譜是具有三個頻率變量的函數(shù)。從而，我們可以了解高階譜的含義。高階譜也稱為多譜，就是有多個頻率變量的譜。習(xí)慣上，功率譜用表示，雙譜用表示，而三譜用表示。高階累積量和高階譜具

35、有以下特點：r 假設(shè)，其中和是兩個相互獨(dú)立的隨機(jī)過程，則 ；r 若是高斯隨機(jī)信號，則，這里，此時如果是非高斯信號，且與獨(dú)立，則利用高階累積量可從混合信號中恢復(fù)出來；r 若是線性過程，即，其中是獨(dú)立同分布隨機(jī)過程，那么相應(yīng)的功率譜、雙譜和三譜的表達(dá)式分別如下，其中。從表達(dá)式可以看出，功率譜并不包含任何的相位信息，而高階譜則包含了相位信息。當(dāng)是非高斯隨機(jī)過程時，其相位信息就可從高階譜中恢復(fù)。根據(jù)以上特點，我們知道相對于二階統(tǒng)計量，高階累積量和高階譜具有其自身的優(yōu)點，能夠彌補(bǔ)二階統(tǒng)計量在某些方面的不足。具體來說，高階統(tǒng)計量適用于以下兩方面的情況：1) 加性噪聲是高斯的，而信號是非高斯的，此時就可利用

36、高階統(tǒng)計量將信號從混合信號中恢復(fù)出來；2) 線性系統(tǒng)不是最小相位情況，或者說隨機(jī)過程是非線性的，此時高階譜可以恢復(fù)出信號的相位信息。7.4.2 高階累積量估計對于高階累積量和高階譜的估計，我們這里所應(yīng)用的是高階譜分析工具箱（High-order Spectral Analysis Toolbox）。該工具箱功能比較強(qiáng)大，包括了很多先進(jìn)的信號處理算法，如譜估計、高階譜估計、時頻分布計算等，具體應(yīng)用領(lǐng)域又包括參數(shù)和非參數(shù)的盲系統(tǒng)辨識、時延估計、諧波恢復(fù)、波達(dá)方位估計、Volterra模型參數(shù)估計、自適應(yīng)線性預(yù)測。一般觀測到的樣本數(shù)都為有限長度，如，那么利用觀測到的樣本值來估計二階、三階和四階累積量

37、的公式分別如下：這里 式中和是選取的樣本數(shù)目，要求當(dāng)時，只對求和，選取合適的是為了獲取無偏估計。令公式中，就會得到自累積量的估計式。工具箱提供了函數(shù)cum2x、cum3x、cum4x分別估計二階、三階和四階累積量，函數(shù)cumest用來估計自累積量。調(diào)用格式分別如下：l cvec = cum2x(x,y,maxlag,samp_seg,overlap,flag);l cvec = cum3x(x,y,z,maxlag,samp_seg,overlap,flag,k1)l cvec = cum4x(w,x,y,z,maxlag,samp_seg,overlap,flag,k1,k2)l cvec

38、= cumest(y,norder,maxlag,samp_seg,overlap,flag,k1,k2);各參數(shù)的說明如下：r y是樣本矩陣或向量；r norder指定了累積量階數(shù)，可以取2、3、4三個值，默認(rèn)值為2；r maxlag指定了累積量允許的最大延遲，默認(rèn)值為0；r samp_seg規(guī)定每個數(shù)據(jù)段的樣本數(shù)，默認(rèn)值是樣本序列長度；r overlap指定數(shù)據(jù)段間重疊的百分比，允許的最大值為99，默認(rèn)值是0；r flag包括biased和unbiased兩種情況，如果為biased，則計算有偏樣本估計，否則計算無偏樣本估計，默認(rèn)為biased；r 如果y為矩陣，則每一列單獨(dú)進(jìn)行累積量估計

39、運(yùn)算，此時overlap需設(shè)為零，samp_seg為行維數(shù)，最終累積量為每列估計值的平均結(jié)果；r 參數(shù)k1和k2決定計算累積量函數(shù)的哪個一維切片，默認(rèn)值都為零。r 產(chǎn)生一頻率為300Hz的諧波信號（復(fù)數(shù)形式），疊加高斯白噪聲，試用特征向量法估計混合信號的功率譜密度。例7-11模擬一個非高斯ARMA隨機(jī)過程，并估計它的3階累積量，采用函數(shù)cum3x。% 函數(shù)cum3x估計3階累積量u = rpiid(1024,'exp');n = 25;y = filter(1,-2,1,-1.5,0.8,u);for k = -n:n cmat(:,k+n+1) = cum3x(y,y,y,n

40、,128,20,'biased',k);endfigure;contour(-n:n,-n:n,cmat,8)函數(shù)rpiid產(chǎn)生一個獨(dú)立同分布隨機(jī)序列。這里樣本序列被分割成長度128的數(shù)據(jù)段，重疊百分比為20%，對每個數(shù)據(jù)段計算有偏的3階累積量，再計算平均值。程序運(yùn)行結(jié)果如圖7-8(a)所示，等高線圖顯示了三階累積量的基本對稱性。例7-12采用函數(shù)cumest估計例11中隨機(jī)過程的3階累積量。% 函數(shù)cumest估計3階累積量u = rpiid(1024,'exp');n = 25;y = filter(1,-2,1,-1.5,0.8,u);for k = -n

41、:n cmat(:,k+n+1) = cumest(y,3,n,128,20,'biased',k);endfigure;contour(-n:n,-n:n,cmat,8)程序運(yùn)行結(jié)果如圖7-8(b)所示，與圖(a)相比較可以看出，函數(shù)cum3x和函數(shù)cumest的估計結(jié)果基本一致，二者都體現(xiàn)了三階累積量的基本對稱性。(a) 函數(shù)cum3x估計結(jié)果 (b) 函數(shù)cumest估計結(jié)果圖7-8 一ARMA(2,1)隨機(jī)過程的三階累積量估計7.4.3 高階譜估計最基本的互雙譜估計方法是對三階累積量做Fourier變換，即這里是隨機(jī)序列，的Fourier變換。上面的估計方法并非一致性估

42、計，通常需要采取平滑的方法獲取一致性。當(dāng)時，就是雙譜估計表達(dá)式。所謂平滑處理，一般通過窗函數(shù)乘以三階累積量來完成。假設(shè)是一個二維窗函數(shù)，其二維Fourier變換有界且非負(fù)，并且假設(shè)同時，窗函數(shù)也要滿足三階累積量的對稱性。平滑后的互雙譜估計式為：高階譜分析工具箱提供了不同方法實現(xiàn)高階譜估計的函數(shù)，包括：r 函數(shù)bispecd，直接法雙譜估計r 函數(shù)bispeci，間接法雙譜估計r 函數(shù)bispect，ARMA模型的理論雙譜估計對于這些函數(shù)及相應(yīng)的算法，接下來依次進(jìn)行介紹。7.4.3.1 函數(shù)bispecd 函數(shù)bispecd實現(xiàn)直接法雙譜估計，算法的實現(xiàn)步驟如下： 步驟r 將樣本數(shù)據(jù)分成段，每段

43、含有個樣本值，記作，其中，這里允許兩段相鄰數(shù)據(jù)之間有重疊；r 計算離散Fourier變換（DFT）系數(shù)，其中，；r 計算DFT系數(shù)的三重相關(guān)，其中，和應(yīng)選擇滿足的值；r 段雙譜估計的平均值即為樣本數(shù)據(jù)的最終雙譜估計值；，函數(shù)bispecd的調(diào)用格式為：l bspec,waxis = bispecd(y)l bspec,waxis = bispecd(y,nfft,wind,samp_seg,overlap)各參數(shù)的說明如下：r y是樣本矩陣或向量；r nfft為計算FFT的長度，一般取2的整數(shù)次冪，常規(guī)默認(rèn)值為128；r wind指定頻域平滑窗參數(shù)；1) 如果wind是標(biāo)量，則采用長度為此標(biāo)量

44、值的Rao-Gabr窗；2) 如果wind等于1，則不采用平滑窗；3) 如果wind為向量，則計算二維窗；4) 如果wind為矩陣，則直接指定二維平滑窗；r samp_seg規(guī)定每個數(shù)據(jù)段的樣本數(shù)，默認(rèn)值為8；r overlap指定數(shù)據(jù)段間重疊的百分比，默認(rèn)值是50；r 返回值bspec為雙譜的估計結(jié)果，大小為nfft×nfft，其原點在中心位置，軸的方向是向下與向右；r 返回值waxis是與bspec矩陣行和列相關(guān)的頻率向量；7.4.3.2 函數(shù)bispeci 函數(shù)bispeci實現(xiàn)間接法雙譜估計，算法的實現(xiàn)步驟為： 步驟r 將原個樣本數(shù)據(jù)分成段，每段含有個樣本值；r 設(shè)為第段數(shù)據(jù)

45、，計算各段的三階累積量估計值，其中，；r 取所有段的三階累積量的平均值作為整個觀測數(shù)據(jù)組的三階累積量估計；r 計算雙譜估計，其中，為二維滯后窗函數(shù)；函數(shù)bispeci的調(diào)用格式為：l bspec,waxis = bispeci(y,nlag)l bspec,waxis = bispeci(y,nlag,samp_seg,overlap,flag,nfft,wind)參數(shù)說明如下：r y是樣本矩陣或向量；r nlag規(guī)定了累積量計算延遲數(shù)，假設(shè)計算三階累積量，則要求， ，該參數(shù)必須指定，一般可將該值取為樣本序列 長度的十分之一。r samp_seg規(guī)定每個數(shù)據(jù)段的樣本數(shù)，默認(rèn)值為觀測序列長度；r

46、 overlap指定數(shù)據(jù)段間重疊的百分比，取值范圍0,99，默認(rèn)值為0；r flag包括biased和unbiased兩種情況，如果為biased，則計算有偏樣本估計，否則計算無偏樣本估計，默認(rèn)為biased；r nfft為計算FFT的長度，一般取2的整數(shù)次冪，常規(guī)默認(rèn)值為128；r wind指定平滑窗參數(shù)，當(dāng)wind = 0時，采用Parzen窗，否則采用六角形窗，默認(rèn)Parzen窗。r 返回值bspec為雙譜的估計結(jié)果，大小為nfft×nfft，其原點在中心位置，軸的方向是向下與向右；r 返回值waxis是與bspec矩陣行和列相關(guān)的頻率向量；例7-13分別利用直接法和間接法估計

47、一個已知的非高斯信號的雙譜。% 直接法雙譜估計load qpc;figure(1)bspecd = bispecd(zmat,128,5,64,50);% 間接法雙譜估計figure(2)bspeci = bispeci(zmat,25,64,50,'biased',128,5);程序中采用了MATLAB自帶的非高斯檢測信號qpc.mat，直接法估計函數(shù)bispecd的參數(shù)設(shè)置中，將FFT長度確定為128，每段數(shù)據(jù)長64，并且重疊百分比為50%，這些與間接法估計函數(shù)bispeci的設(shè)置相同，同時間接法采取了有偏估計。程序運(yùn)行結(jié)果如圖7-9所示，(a)圖為直接法估計結(jié)果，(b)圖為間接法估計結(jié)果，可以看出，兩種方法得到的峰值位置相同，結(jié)果基本一致。(a)直接法 (b)間接法圖7-9 直接法和間接法雙譜估計圖7.4.3.3 函數(shù)bispect 函數(shù)bispect采用ARMA模型實現(xiàn)理論雙譜的估計。ARMA模型經(jīng)常用于描述一個隨機(jī)序列，表達(dá)式為式中是獨(dú)立同分布序列，方差。自回歸（AR）多項式和滑動平均（MA）多項式分別定義為，利用ARMA模型估計理論雙譜，關(guān)鍵是對模型參數(shù)和的估