【正文】 方便，因?yàn)閷?duì)一個(gè)樣本過程進(jìn)行長(zhǎng)時(shí)間統(tǒng)計(jì)比同時(shí)對(duì)許多樣本進(jìn)行統(tǒng)計(jì)要容易實(shí)現(xiàn) 。本文計(jì)算所涉及的隨機(jī)數(shù)據(jù)均看作廣義平穩(wěn)和各態(tài)歷經(jīng)的隨機(jī)過程。 隨機(jī)信號(hào)的相關(guān)函數(shù) 平穩(wěn)隨機(jī)信號(hào)的自相關(guān)函數(shù) ()xRm定義為 *( ) [ ( ) ( ) ]xR m E x n x n m?? (46) 式中 ()xn 是平穩(wěn)隨機(jī)信號(hào)， * 代表取共軛，如果是 ()xn 實(shí)隨機(jī)信號(hào)，則上式成為 ( ) [ ( ) ( ) ]xR m E x n x n m?? (47) 平穩(wěn)隨機(jī)信號(hào)的相關(guān)函數(shù)具有如下性質(zhì)： ① 若 ()xn 為實(shí)信號(hào)，則 ()xRm為實(shí)偶函數(shù)，即 *( ) ( )xxR m R m? (48) ( ) ( )xxR m R m?? (49) 若 ()xn 為復(fù)信號(hào)，則 ()xRm共軛偶對(duì)稱， 即 *( ) ( )xxR m R m?? (410) 本文所用到的隨機(jī)信號(hào)均為實(shí)隨機(jī)信號(hào)。 ② 在 0m? 時(shí)， ()xRm取得最大值，即 (0) ( )xxR R m? 且 (0)xR 就是序列的平均功率，即 2(0) [ ( )]xR E x n? (411) ③ 維納 辛欽定理 隨機(jī)信號(hào)自相關(guān)函數(shù)的傅立葉變換是信號(hào)的功率譜密度。如果用 ()jwxSe 表示隨機(jī)信號(hào)序列 ()xn 的功率譜密度，則有 ( ) [ ( ) ] ( )jw jw mx x xmS e F R m R m e? ?? ? ??? ? (412) 1 1( ) [ ( ) ] ( )2jw jw jw mx x xR m F S e S e e d?? ??? ??? ? (413) 在式中令 0m? ，可得 1( 0 ) ( )2 jxxR S e d? ?? ?? ?? ? (414) 因?yàn)?2(0) [ ( )]xR E x n? 為信號(hào)的平均功率，所以將 ()jxSe? 稱為隨機(jī)信號(hào)的功率譜密度 (PSD)或功率譜密度譜，簡(jiǎn)稱功率譜。 隨機(jī)信號(hào)的功率譜 對(duì)確定性能量信號(hào)，可以用 FFT（快速傅立葉變換） 作頻譜分析，得到頻域特性。對(duì)于平穩(wěn)隨機(jī)信號(hào)，因?yàn)槭菬o限能量的信號(hào)，故其傅立葉變換不存在（在 z平面的單位圓上不滿足絕對(duì)可和的條件）。由于自相關(guān)函數(shù)序列 ()xRm是一個(gè)能量有限的確定性序列，故能滿足序列傅 立葉變換絕對(duì)可和的條件。由維納 辛欽定理可知，對(duì)序列 ()xRm求傅立葉變換得到的就是序列的功率譜 ()jxSe? 。 平穩(wěn)隨機(jī)信號(hào)的功率譜具有如下性質(zhì)： ① 不論 ()xn 是實(shí)序列還是復(fù)序列， ()jxSe? 都是 ? 的實(shí)函數(shù)。 ② 如果 ()xn 是實(shí)序列， ()jxSe? 具有偶對(duì)稱性。 ③ ()jxSe? 對(duì)所有的 ? 都是非負(fù)的，且是 ? 的周期函數(shù)，周期為 2? 。 功率譜估計(jì)的各種方法 根據(jù)維納 辛欽定理，廣義平穩(wěn)隨機(jī)過程的功率譜是自相關(guān)函數(shù)的傅立葉變換，它取決于無限多個(gè)自相關(guān)函數(shù)值。但對(duì)于許多實(shí)際應(yīng)用問題，可資利用的觀測(cè)數(shù)據(jù)往往是有限的，所以要準(zhǔn)確計(jì)算功率譜通常是不可能的。比較合理的目標(biāo)是 設(shè)法得到功率譜的一個(gè)好的估計(jì)值，這就是功率譜估計(jì)。也就是說，功率譜估計(jì)是根據(jù)平穩(wěn)隨機(jī)過程一個(gè)實(shí)現(xiàn)的有限個(gè)觀測(cè)值，來估計(jì)隨機(jī)過程的功率譜密度（ PSD）。這里涉及到兩個(gè)問題，怎樣評(píng)價(jià)一個(gè)估計(jì)是好的估計(jì)？怎樣得到好的估計(jì)？ 功率譜估計(jì)評(píng)價(jià)指標(biāo)包括客觀度量和統(tǒng)計(jì)度量。在客觀度量中，譜分辨率特性是一個(gè)主要指標(biāo)。譜分辨率是指估計(jì)譜對(duì)真實(shí) 譜 中兩個(gè)靠的很近的譜峰的分辨能力。統(tǒng)計(jì)度量是指估計(jì)的偏差、方差、均方誤差、一致性等評(píng)價(jià)指標(biāo)。但需要注意的是，對(duì)統(tǒng)計(jì)特性的分析方法只適用于長(zhǎng)數(shù)據(jù)記錄。所以，利用統(tǒng)計(jì)度量對(duì)不同的譜估計(jì)方法進(jìn) 行比較是不妥當(dāng)?shù)?，只能用來?duì)某種譜估計(jì)方法進(jìn)行描述，并且一般只用來描述古典譜估計(jì)方法，因?yàn)楝F(xiàn)代譜估計(jì)方法往往用于短數(shù)據(jù)記錄情況。 至于怎樣得到好的估計(jì)，這就是后面將要介紹的各種譜估計(jì)方法。這些方法主要分為兩大類。通常，將以傅立葉分析為理論基礎(chǔ)的譜估計(jì)方法叫做古典譜估計(jì)方法或經(jīng)典譜估計(jì)方法；而把不同于傅立葉分析的新的譜估計(jì)方法叫做現(xiàn)代譜估計(jì)或近代譜估計(jì)。 本文主要介紹了古典譜估計(jì)和最大熵譜估計(jì)，分述如下： 古典譜估計(jì) 古典譜估計(jì)主要有相關(guān)法（間接法）和周期圖法（直接法）兩種，以及由此派生出來的各種 改進(jìn)算法。 相關(guān)法譜估計(jì) 相關(guān)法譜估計(jì)是以相關(guān)函數(shù)為媒介來計(jì)算功率譜，所以又叫間接法。它的理論基礎(chǔ)是維納 辛欽定理，因?yàn)槭怯?Blackman 和 Turkey 提出的，所以又叫 BT法。其具體步驟如下： 第一步：由獲得的 N 點(diǎn)數(shù)據(jù)構(gòu)成的有限長(zhǎng)序列 ()Nxn來估計(jì)自相關(guān)函數(shù) ()xRm?序列，即 101( ) ( ) ( )Nx NNnR m x n x n mN?????? (415) 這一步要注意 三個(gè)問題： ① ()xRm? 是雙邊序列，自變量的取值范圍應(yīng)為( 1 ) , 1 , 0 , 1 , , ( 1 )m N N? ? ? ???? ??? ?。對(duì)實(shí)數(shù)序列 ()Nxn，由自相關(guān)函數(shù)的偶對(duì)稱性，只需求出 0,1, , ( 1)mN? ??? ?的 ()xRm? ，另一半也就知道了。 ② 因?yàn)榈玫降闹皇?()xn 的 N個(gè)觀測(cè)值 ( 0) , (1) , ( 1)N N Nx x x N??? ?，所以對(duì) nN? 時(shí)的 ()xn 的值只能假設(shè)為零。 ③ 為了較少地使用已知數(shù)據(jù)段之外的數(shù)據(jù)，通常取 mN?? 。所以，計(jì)算自相關(guān)函數(shù)時(shí)只需求出 0,1, , ( 1)mm? ??? ?的 ()xRm? （這里 MN?? ），然后利用 ( ) ( )xxR m R m????求出另一半。 第二步：由自相關(guān)函數(shù)的傅立葉變換求功率譜， 即 1( 1 )( ) ( )Mj j mxxmMS e R m e????? ?? ? ?? ? (416) 以上兩步中經(jīng)歷了兩次截?cái)?，一次是估?jì) ()xRm? 時(shí)僅利用了 ()xn 的 N 個(gè)觀測(cè)值()Nxn，這相當(dāng)于對(duì) ()xn 加矩形窗截?cái)唷Ｔ摯笆羌釉跀?shù)據(jù)上的，一般稱為加數(shù)據(jù)窗。另一次是估計(jì) ()jxSe?? 時(shí)僅利用了從 ( 1)M??到 ( 1)M? 的 ()xRm? ，這相當(dāng)于對(duì) ()xRm加矩形窗截?cái)?，?()xRm截成 (2 1)M? 長(zhǎng)，這稱為加延遲窗。式中的 ()xRm? 和 ()jxSe??分別表示對(duì) ()xRm和 ()jxSe? 的估值。由于 MN?? ，使得式 (416)的運(yùn)算量不是很大 ， 因此，在 FFT 問世之前，相關(guān)法是最常用的譜估計(jì)方法。 當(dāng) FFT 問世后情況有所變化。因?yàn)榻財(cái)嗪蟮?()Nxn可視作能量信號(hào)，由相關(guān)卷積定理，可將式 (415)寫為 1( ) [ ( ) ( ) ]x NNR m x m x mN? ? ? ? (417) 對(duì)上式兩邊取 (2 1)N? 點(diǎn)的 DFT，并將時(shí)域卷積變?yōu)轭l域乘積，有 2*2 1 2 1 2 111( ) [ ( ) ( ) ] ( )jx N N NS e X k X k X kNN?? ? ? ?? ? ? (418) 于是可得用 FFT 求 ()xRm? 的完整方案如下： ① 對(duì) N 長(zhǎng)的 ()Nxn充 ( 1)N? 個(gè)零，成為 (2 1)N? 長(zhǎng)的 21()nxn? 。 ② 求 (2 1)N? 點(diǎn)的 FFT，得 222 1 2 1 2 10( ) ( )N nkN n NnX k x n W?? ? ??? ? (419) ③ 求 2211 ()NXkN ?。 ④ 求 (2 1)N? 點(diǎn)的 IFFT，得 1 22 1 2 1( 1 )11( ) ( )21N nkxNNkNR m X k WNN?? ???? ? ?? ? ? (420) ⑤ 由式 (416)便可求出 ()jxSe?? 。 上面的相關(guān)運(yùn)算中，充零的目的是為了能用圓周卷積代替線性卷積，以便采用快速卷積算法。 周期圖法譜估計(jì) 周期圖法又稱直接法，其具體步驟如下： 第一步：由獲得的 N 點(diǎn)數(shù)據(jù)構(gòu)成的有限長(zhǎng)序列 ()Nxn直接求傅立葉變換，得頻譜 ()jNXe? ，即 10( ) ( )Nj j mNNnX e x n e??? ??? ? (421) 第二步：取頻譜幅度的平方，并除以 N ，以此作為對(duì) ()xn 真實(shí)功率譜 ()jxSe? 的估計(jì)，即 21( ) ( )jjx NS e X eN??? ? (422) 古典譜估計(jì)的改進(jìn) 相關(guān)法和周期圖法都是用獲得的 N 個(gè)數(shù)據(jù)對(duì)隨機(jī)過程進(jìn)行譜估計(jì)，它隱含著對(duì)無限長(zhǎng)數(shù)據(jù)加了一個(gè)矩形窗截?cái)?。時(shí)域中與矩形窗函數(shù)相乘對(duì)應(yīng)于頻域中與矩形窗頻譜相卷積，就這一點(diǎn)來說，估計(jì)譜相當(dāng)于真實(shí)譜與矩形窗頻譜相卷積的結(jié)果。矩形窗頻譜等于矩形序列 ()NRn的傅立葉變換 ，即 1120sin ( )1 2( ) ( ) 1sin ( )2NjNN jj j n j nRN jnnNeW e R n e e ee? ?? ? ???????? ????? ? ? ??? ? ? ???? (423) 對(duì)功率譜有影響的是矩形窗頻譜的幅度譜 sin( ) / sin( )22N??，該函數(shù)與功率譜相卷積必然使所得的估計(jì)譜不同于真實(shí)的功率譜，因?yàn)楹瘮?shù)不同于 ? 函數(shù)（任何函數(shù)與 ? 函數(shù)相卷積形狀不變）。 為了使估計(jì)譜逼近真實(shí)譜，應(yīng)設(shè)法使窗譜幅度函數(shù)逼近 ? 函數(shù)。矩形窗頻譜的幅度函數(shù)為 sin( ) / sin( )22N??，它與 ? 函數(shù)相比存在兩方面的差異，一是主瓣不是無限窄，二是有旁瓣。由于主瓣不是無限窄，真實(shí)的功率譜與主瓣卷積后將使功率向附近頻域擴(kuò)散，使信號(hào)模糊，降低了譜分辨率，主瓣越寬分辨率越低。由于存在旁瓣，又會(huì)產(chǎn)生兩個(gè)后果，一是功率譜主瓣內(nèi)的能量“泄漏”到旁瓣使譜估計(jì)的方差增大，二是旁瓣卷積后得到的功率譜完全屬于干擾。嚴(yán)重情況下，強(qiáng)信號(hào)與旁瓣的卷積可能大于弱信號(hào)與主瓣的卷積，使弱信號(hào)淹沒在強(qiáng)信號(hào)的干擾中，無法檢測(cè)出來，這是古典譜估計(jì)的主要缺點(diǎn)。 為了減少泄漏和提高譜估計(jì)的分辨率，改善窗的形狀是必要的 。對(duì)窗函數(shù)幅度譜的總要求是希望它的主瓣盡可能窄，同時(shí)旁瓣盡可能小，以減少能量外泄，改善方差，提高譜估計(jì)的分辨率。但實(shí)現(xiàn)時(shí)兩者往往是矛盾的，任何減小旁瓣的努力都要以犧牲主瓣寬度為代價(jià)，反之亦然。兩者只能互換，不可兼得。因此古典譜估計(jì)無法克服譜分辨率低的缺點(diǎn)。如果用統(tǒng)計(jì)量度量進(jìn)行評(píng)價(jià)，可以得出相關(guān)法譜估計(jì)和周期圖法譜估計(jì)都不是對(duì) ()jxSe? 的一致估計(jì)，主要問題是方差大。于是人們研究各種改進(jìn)方法以減少方差，尤其是周期圖法。改進(jìn)周期圖法的主要途徑是平滑和平均。 平滑是用一個(gè)適當(dāng)?shù)拇昂瘮?shù) ()jWe? 與算出的功率譜 ()jxSe?? 進(jìn)行卷積，使譜線平滑。這種方法得出的譜估計(jì)是無偏的，方差也小，但分辨率下降。平均就是將截取的數(shù)據(jù)段 ()Nxn再分成 L 個(gè)小段，分別計(jì)算功率譜后取功率譜的平均，這種方法叫 Bartlett 方法。因?yàn)?L 個(gè)平均的方差比隨機(jī)變量的單獨(dú)方差小 L 倍， 所以當(dāng)L?? 時(shí)， L 個(gè)平均的方差趨于零，可以達(dá)到一致譜估計(jì)的目的。 現(xiàn)在比較常用的方法是 Welch 法，又叫加權(quán)交疊平均法，簡(jiǎn)記為 WOSA 法。這種方法以加窗（加權(quán)）求取平均，以分段重疊求得平均，因此集平均與平滑的優(yōu)點(diǎn)于一體，但同時(shí)也不可避免地帶有兩者的缺點(diǎn)。其主要步驟如下： ① 將 N 長(zhǎng)的數(shù)據(jù)段分成 L 個(gè)小段，每小段 M 點(diǎn)，相鄰小段間交疊 /2M 點(diǎn)，于是段數(shù)為 /2/2NML M?? (424) ② 對(duì)各小段加同樣的平滑窗 ()wn 后求傅立葉變換，得 10( ) ( ) ( ) 1 , ,Mj j niinX e x n w n e i L??? ??? ? ???? ， (425) ③ 求各小段功率譜的平均，得 2111( ) ( )Ljjx iiS e X eL M U????? ? (426) 這里， 1 201 ()MnU w nM??? ?代表窗函數(shù)平均功率，所以 1 20 ()MnMU w n????是 M 長(zhǎng)窗函數(shù)()wn 的能量。 Welch 法對(duì)數(shù)據(jù)分段加非矩形窗，因?yàn)榇昂瘮?shù)在其邊沿為零，從而減小了分段時(shí)的截?cái)嘈?yīng)。分段平均減小了由數(shù)據(jù)樣本本身的隨機(jī)性帶來的方差，段數(shù)越多方差越小，但分辨率下降。分段多了每段的點(diǎn)數(shù)必然減少，分段時(shí)允許數(shù)據(jù)有部分重疊，可以在段數(shù)較多