無需參考信號(hào)的單晶爐溫度干擾抑制新方法

張新雨，劉丁，梁軍利，汪姣，楊文

(西安理工大學(xué) 自動(dòng)化與信息工程學(xué)院，陜西 西安 710048)

單晶爐是利用直拉法生產(chǎn)集成電路硅材料的重要設(shè)備。在其直拉單晶硅的工作過程中，熱是根本驅(qū)動(dòng)力，影響著單晶硅的整個(gè)生長(zhǎng)過程。隨著集成電路線距達(dá)到低納米級(jí)，對(duì)單晶爐固液界面處的溫度檢測(cè)與控制提出更為苛刻的要求，因?yàn)闇囟劝l(fā)生細(xì)微的波動(dòng)都會(huì)導(dǎo)致晶體發(fā)生位錯(cuò)而造成拉晶失敗[1-2]。然而在單晶爐的工作過程中，坩堝轉(zhuǎn)動(dòng)、加熱器紋波、晶體轉(zhuǎn)動(dòng)等變化都會(huì)使得熱場(chǎng)溫度測(cè)量值存在周期性低頻干擾，較高精度的溫度控制變得較為困難。因此，如何在存在多個(gè)周期性低頻干擾的情況下，對(duì)單晶爐熱場(chǎng)溫度測(cè)量值中的干擾實(shí)現(xiàn)抑制以達(dá)到熱場(chǎng)溫度控制系統(tǒng)的需求，就顯得極為關(guān)鍵了。

目前，國(guó)內(nèi)外大部分單晶爐采用紅外傳感器檢測(cè)熱場(chǎng)溫度，然后運(yùn)用中值濾波、均值濾波等工程方法進(jìn)行處理，這些方法雖對(duì)溫度測(cè)量值中的干擾具有一定抑制作用，但要滿足高精度的控制仍存在不足，檢測(cè)精度偏低，結(jié)果仍具有一定弱周期性。近年來，國(guó)內(nèi)學(xué)者對(duì)單晶爐熱場(chǎng)溫度測(cè)量數(shù)據(jù)進(jìn)行了研究分析，發(fā)現(xiàn)單晶爐熱場(chǎng)溫度測(cè)量值中存在多個(gè)周期性的低頻干擾[3]。自適應(yīng)濾波方法能夠抑制干擾，但如何確定自適應(yīng)濾波器的參考信號(hào)是一個(gè)重要問題[4]。目前，自適應(yīng)濾波器的參考信號(hào)主要有兩種[5-10]：一種采用正弦信號(hào)，即自適應(yīng)陷波器，能夠針對(duì)某一固定頻率的干擾進(jìn)行抑制，但需要已知干擾頻率[5-6]；另一種采用與干擾相關(guān)的函數(shù)，即自適應(yīng)對(duì)消法[7]，利用自適應(yīng)算法根據(jù)參考信號(hào)抵消干擾[3,8-10]。以上方法均需要給定參考信號(hào)，但在實(shí)際應(yīng)用過程中，會(huì)存在以下問題：

1)參考信號(hào)需要額外給定，一般很難準(zhǔn)確獲得；

2)在實(shí)際的工作環(huán)境中，由于條件限制，很難布置測(cè)量參考信號(hào)的傳感器；

3)在測(cè)量參考信號(hào)的實(shí)際過程中，傳感器難以僅采集到參考信號(hào)(例如干擾)，往往還包含了部分有用信號(hào)(例如溫度)，因此在對(duì)消時(shí)，有用信號(hào)很容易被當(dāng)作干擾對(duì)消掉。

近幾年，迭代自適應(yīng)法(iterative adaptive approach, IAA)作為一種非參數(shù)化的譜估計(jì)方法被提出[11]，為自適應(yīng)對(duì)消法的發(fā)展提供了新的思路。該方法以其估計(jì)單個(gè)分量時(shí)不會(huì)被其它分量影響、幅值估計(jì)精度高、可用于少量數(shù)據(jù)等優(yōu)點(diǎn)，被廣泛應(yīng)用于陣列信號(hào)處理[12]、水下通信[13]和圖像處理[14]等領(lǐng)域。

為此，本文針對(duì)單晶爐熱場(chǎng)溫度測(cè)量值中存在多個(gè)周期性低頻干擾的問題，提出了一種無需參考信號(hào)的新的自適應(yīng)對(duì)消法。該方法將迭代自適應(yīng)法引入自適應(yīng)對(duì)消法中，直接利用迭代自適應(yīng)法從輸入信號(hào)中逐個(gè)獲得低頻干擾分量，然后將其消除。該方法相比傳統(tǒng)的自適應(yīng)對(duì)消法具有以下優(yōu)點(diǎn)：1)當(dāng)干擾為多個(gè)窄帶干擾時(shí)，無需參考信號(hào)，降低了測(cè)量參考信號(hào)所需的硬件代價(jià)；2)估計(jì)每個(gè)干擾分量的過程中能夠減小其它干擾分量的影響；3)迭代次數(shù)少。最后通過仿真和工程實(shí)驗(yàn)對(duì)該方法進(jìn)行了有效性驗(yàn)證。實(shí)驗(yàn)結(jié)果表明，在實(shí)際的運(yùn)行過程中，該方法在無需參考信號(hào)的條件下，對(duì)溫度信號(hào)進(jìn)行檢測(cè)和處理時(shí)，能夠有效地抑制單晶爐熱場(chǎng)溫度測(cè)量值中的干擾。

1 問題描述

在單晶爐直拉單晶硅的過程中，爐膛內(nèi)部高溫、負(fù)壓的環(huán)境導(dǎo)致無法對(duì)整個(gè)熱場(chǎng)進(jìn)行測(cè)量。因此，一般采用紅外傳感器通過爐壁開孔獲得固液界面平面處加熱器的溫度，進(jìn)而間接獲得熱場(chǎng)溫度，用于單晶爐的熱場(chǎng)溫度控制。單晶爐熱場(chǎng)溫度測(cè)量原理如圖1所示。

圖1 單晶爐熱場(chǎng)溫度測(cè)量示意圖

直拉單晶硅的過程中，熱場(chǎng)溫度信號(hào)變化較為緩慢。由于坩堝轉(zhuǎn)動(dòng)、晶體轉(zhuǎn)動(dòng)、加熱器電源紋波等都是周期性的干擾，加之紅外傳感器本身的測(cè)量噪聲，因此，溫度傳感器獲得的測(cè)量值不能直接用于單晶爐熱場(chǎng)溫度控制系統(tǒng)，需要對(duì)上述干擾進(jìn)行濾波。根據(jù)上述分析，溫度傳感器獲得的測(cè)量值可以寫成：

(1)

式中，y(n)為傳感器測(cè)量得到的溫度信號(hào)，x(n)為緩慢變化的實(shí)際溫度信號(hào)，ωi為干擾頻率，ai和bi為干擾幅值，Δt為采樣周期，k為低頻干擾個(gè)數(shù),v(n)為測(cè)量噪聲，根據(jù)檢測(cè)理論[15]和實(shí)驗(yàn)證明，v(n)近似服從高斯分布。

2 新的自適應(yīng)對(duì)消法及單晶爐熱場(chǎng)溫度干擾抑制

2.1 建 模

設(shè)傳感器測(cè)量的N個(gè)數(shù)據(jù)組成一個(gè)觀察序列，則令:

由于一個(gè)頻率分量對(duì)應(yīng)一對(duì)正弦和余弦，因此，現(xiàn)將每一個(gè)頻率分量所對(duì)應(yīng)正弦和余弦的幅值定義為一個(gè)列向量，即:

αi=[a(ωi),b(ωi)]T

將ci和si寫成一個(gè)矩陣Di=[ci,si]，由此，式(1)可以重寫為:

(2)

2.2 自適應(yīng)對(duì)消法原理

本文最終的目的是從y中估計(jì)x，然而，ωi、αi、Di、k和v均為未知，直接從y中估計(jì)x極為困難，因此，需采用自適應(yīng)對(duì)消法的思想，通過估計(jì)ωi、αi、Di、k得到干擾，然后從y中將其消除，從而得到x的估計(jì)值。傳統(tǒng)的自適應(yīng)對(duì)消法如圖2所示[7]，自適應(yīng)對(duì)消器具有兩個(gè)輸入通道，一個(gè)為信號(hào)通道，即同時(shí)含有有用信號(hào)和干擾，一個(gè)為參考通道，主要為與干擾相關(guān)的函數(shù)。參考信號(hào)經(jīng)過自適應(yīng)算法運(yùn)算后，送入自適應(yīng)對(duì)消器，最后將其從信號(hào)源中消除。

圖2 傳統(tǒng)自適應(yīng)對(duì)消法框圖

盡管傳統(tǒng)的自適應(yīng)對(duì)消法理論較為完備，但是參考信號(hào)一般需要額外測(cè)量。在單晶爐工作過程中，這些周期性干擾源難以確定，且爐膛內(nèi)部負(fù)壓、高溫的環(huán)境，較難布置測(cè)量參考信號(hào)的傳感器，因此，采用傳統(tǒng)的自適應(yīng)對(duì)消法抑制干擾較為困難。

圖3 本文自適應(yīng)對(duì)消法框圖

由圖2和圖3可以得出，與傳統(tǒng)的自適應(yīng)對(duì)消法相比，本文方法具有以下優(yōu)點(diǎn)：1)無需參考信號(hào)，可以減小測(cè)量參考信號(hào)的硬件代價(jià)；2)在估計(jì)每個(gè)干擾分量的過程中，可以避免各分量之間相互影響，具有較高的估計(jì)精度和魯棒性。

2.3 基于自適應(yīng)迭代的對(duì)消法

設(shè)定干擾頻率范圍為[ωmin,ωmax]，生成一個(gè)虛擬的頻率網(wǎng)格[11]。頻率網(wǎng)格間距為Δω，則生成網(wǎng)格數(shù)量為:

(3)

通常y中含有多個(gè)頻率分量，為避免其它分量的影響，迭代自適應(yīng)法采用對(duì)每個(gè)頻率分量逐個(gè)求解的思想，即求解某個(gè)頻率分量時(shí)，其它分量都被看作噪聲。因此，迭代過程中選用加權(quán)最小二乘法對(duì)目標(biāo)函數(shù)進(jìn)行求解。根據(jù)加權(quán)最小二乘法，目標(biāo)函數(shù)可寫為:

式中，Qj是求解第j個(gè)頻率分量時(shí)所對(duì)應(yīng)的協(xié)方差矩陣[16]。

(4)

(5)

(6)

根據(jù)加權(quán)最小二乘法，可得:

(7)

根據(jù)矩陣求逆引理[17]得:

(8)

3.2 康復(fù)醫(yī)學(xué)促進(jìn)社會(huì)融合 殘疾人、老年人和其他功能障礙者的社會(huì)參與能力不同程度受限，通過康復(fù)醫(yī)學(xué)的有效干預(yù)，可以促進(jìn)殘疾人的社會(huì)融合。康復(fù)醫(yī)療強(qiáng)調(diào)的是通過功能訓(xùn)練和必要的輔助措施，改善其功能。通過康復(fù)治療改善患者的自理能力，減輕家庭和社會(huì)負(fù)擔(dān)，提高生活質(zhì)量，并能使患者早日回歸社會(huì)?？祻?fù)醫(yī)學(xué)近年也在秉承國(guó)際社會(huì)倡導(dǎo)的全人發(fā)展的理念，在提高個(gè)體功能水平的同時(shí)，越來越多地關(guān)注康復(fù)服務(wù)對(duì)象獲得醫(yī)療康復(fù)服務(wù)后的社會(huì)參與、機(jī)會(huì)均等、權(quán)利保障等問題。

將上式代入式(7)得:

(9)

根據(jù)上述推導(dǎo)過程，式(9)相對(duì)式(7)省略了求Qj的步驟，由式(4)和式(5)及算法復(fù)雜度計(jì)算可得，在一次迭代過程中復(fù)雜度上限降低了O(N3)。

第j個(gè)頻率對(duì)應(yīng)的功率為[16]:

(10)

依次計(jì)算m個(gè)頻率對(duì)應(yīng)的功率，然后迭代10～20次[16]，最后在P(ωj)(j=1,2,…,m)中選取大于設(shè)定閾值為PTH的所有峰值，則峰值個(gè)數(shù)為干擾個(gè)數(shù)k，對(duì)應(yīng)的頻率ωj和幅值αj=[a(ωj),b(ωj)]T為干擾的頻率和幅值，最后從y中將干擾消除。

2.4 算法步驟詳述

本文所提算法步驟如下。

1)設(shè)置算法初始參數(shù)：頻率范圍為[ωmin,ωmax]，頻率網(wǎng)格間距為Δω，最大迭代次數(shù)為L(zhǎng)，功率閾值為PTH；

4)判斷迭代次數(shù)是否大于L，如果大于L，轉(zhuǎn)到步驟 5)，否則轉(zhuǎn)到步驟 3)；

3 實(shí)驗(yàn)結(jié)果與分析

為驗(yàn)證本文方法的有效性，本文應(yīng)用上述算法處理以下數(shù)據(jù)：1)仿真數(shù)據(jù)；2)單晶爐熱場(chǎng)溫度檢測(cè)工程實(shí)驗(yàn)測(cè)量數(shù)據(jù)。在對(duì)單個(gè)觀察序列進(jìn)行實(shí)驗(yàn)時(shí)與經(jīng)典的匹配濾波和傅里葉變換進(jìn)行比較，對(duì)多個(gè)觀察序列進(jìn)行滑動(dòng)估計(jì)實(shí)驗(yàn)時(shí)與經(jīng)典的中值濾波和基于粒子群搜索幅值的自適應(yīng)抵消法(文獻(xiàn)[3]方法)進(jìn)行比較。

3.1 仿真實(shí)驗(yàn)

本次實(shí)驗(yàn)利用下式產(chǎn)生三組仿真數(shù)據(jù):

設(shè)置IAA算法初始參數(shù)，觀察序列數(shù)據(jù)長(zhǎng)度N=64，頻率范圍為fmin=0.01 Hz，fmax=0.49 Hz,頻率網(wǎng)格間距為Δf=0.000 1 Hz，最大迭代次數(shù)L=18，功率閾值為PTH=1。

在信噪比為23、29和31的情況下，對(duì)某個(gè)觀察序列進(jìn)行500次蒙特卡洛實(shí)驗(yàn)，并且與匹配濾波(匹配頻率選擇0.01～0.49 Hz)、傅里葉變換進(jìn)行了比較，圖4～6分別為500次實(shí)驗(yàn)平均之后的歸一化非參數(shù)化估計(jì)譜，由圖4～6可得，與匹配濾波和傅里葉變換相比，本文方法在估計(jì)每個(gè)干擾分量的時(shí)候并不會(huì)受到其它分量的影響，具有較強(qiáng)的旁瓣抑制能力，峰值較尖，對(duì)頻率估計(jì)的分辨率較高。而匹配濾波與傅里葉變換的結(jié)果由于各個(gè)分量之間影響較大，峰值較多。由上述可得，本文方法在不同的信噪比下能夠準(zhǔn)確地搜索到信號(hào)的非參數(shù)譜的峰值，進(jìn)而得到低頻干擾的頻率和幅值，并且500次實(shí)驗(yàn)結(jié)果基本一致，具有較強(qiáng)的魯棒性、穩(wěn)定性和旁瓣抑制能力。

圖4 實(shí)驗(yàn)一μ=23時(shí)非參數(shù)化譜估計(jì)結(jié)果

圖5 實(shí)驗(yàn)一μ=29時(shí)非參數(shù)化譜估計(jì)結(jié)果

圖6 實(shí)驗(yàn)一μ=31時(shí)非參數(shù)化譜估計(jì)結(jié)果

為更清晰地觀察本文方法對(duì)頻率、幅值以及對(duì)消后所得有效信號(hào)的估計(jì)誤差，現(xiàn)定義如下均方誤差(mean-square error，MSE):

其中ξ分別取f1,f2,f3,Ψ1,Ψ2,Ψ3,x。

表1為在三種信噪比下對(duì)頻率、幅值估計(jì)值以及對(duì)消后所得信號(hào)500次實(shí)驗(yàn)的MSE。由表1可得，本文方法對(duì)頻率和幅值估計(jì)具有較高的精度，并且估計(jì)精度隨著信噪比的提高而提高。同時(shí)，本文方法能夠有效地抑制周期性干擾(穩(wěn)定在±0.1內(nèi))。

表1 實(shí)驗(yàn)一本文算法檢測(cè)結(jié)果的MSE

圖7～9分別為中值濾波、文獻(xiàn)[3]方法和本文方法在信噪比為23、29和31時(shí)對(duì)200個(gè)數(shù)據(jù)點(diǎn)進(jìn)行滑動(dòng)濾波的結(jié)果。由圖可得，在不同信噪比下，本文方法的濾波效果明顯優(yōu)于均值濾波，并且在無需參考信號(hào)的條件下，略優(yōu)于文獻(xiàn)[3]方法(需要已知參考信號(hào)的頻率)。

圖7 實(shí)驗(yàn)一μ=23時(shí)不同方法濾波結(jié)果

圖8 實(shí)驗(yàn)一μ=29時(shí)不同方法濾波結(jié)果

圖9 實(shí)驗(yàn)一μ=31時(shí)不同方法濾波結(jié)果

3.2 單晶爐熱場(chǎng)溫度工程檢測(cè)實(shí)驗(yàn)

3.2.1 實(shí)驗(yàn)準(zhǔn)備

本次實(shí)驗(yàn)依托TDR-150型單晶爐進(jìn)行，TDR-150型單晶爐為極大規(guī)模集成電路所用單晶硅的全自動(dòng)生長(zhǎng)設(shè)備，投料量為170 kg，控制直徑為300 mm，等徑長(zhǎng)度為650 mm，TDR-150型單晶爐和溫度傳感器如圖10所示。

圖10 TDR-150型單晶爐及溫度傳感器

根據(jù)單晶爐熱場(chǎng)溫度先驗(yàn)信息設(shè)置參數(shù)：頻率范圍為fmin=0.1 Hz，fmax=4.9 Hz,頻率網(wǎng)格間距為Δf=0.001 Hz，采樣時(shí)間為Δt=0.1 s，功率閾值PTH=1，其余參數(shù)設(shè)置均與實(shí)驗(yàn)一相同。

3.2.2 實(shí)驗(yàn)結(jié)果及分析

圖11為對(duì)250個(gè)數(shù)據(jù)點(diǎn)進(jìn)行滑動(dòng)估計(jì)的平均非參數(shù)化譜(進(jìn)行歸一化之后與匹配濾波、傅里葉變換進(jìn)行了比較)，由圖11可知，本文方法能夠準(zhǔn)確地搜索到非參數(shù)化估計(jì)譜的峰值所對(duì)應(yīng)的干擾頻率，并且在估計(jì)某個(gè)干擾分量時(shí)能夠減小其它分量對(duì)其的影響，對(duì)旁瓣具有較強(qiáng)的抑制能力，且峰值比較尖，具有較高的分辨率。而匹配濾波和傅里葉變換由于估計(jì)每個(gè)分量時(shí)受其它分量干擾較大，不能準(zhǔn)確搜索到干擾頻率所對(duì)應(yīng)的峰值。因此，這兩種方法不能提供有效的干擾進(jìn)行對(duì)消。

圖11 實(shí)驗(yàn)二非參數(shù)化譜估計(jì)結(jié)果

圖12為中值濾波、文獻(xiàn)[3]方法和本文方法對(duì)250個(gè)數(shù)據(jù)點(diǎn)進(jìn)行滑動(dòng)濾波的結(jié)果。

根據(jù)圖12可得，本文方法在沒有參考信號(hào)的條件下，能夠有效抑制單晶爐熱場(chǎng)溫度測(cè)量值中的干擾(在±0.3℃以內(nèi))，并且結(jié)果優(yōu)于中值濾波和基于粒子群搜索幅值的對(duì)消法(干擾頻率已知，濾波結(jié)果在±1℃內(nèi)波動(dòng))，符合TDR-150型單晶爐熱場(chǎng)溫度控制系統(tǒng)的需求。

圖12 實(shí)驗(yàn)二不同方法濾波結(jié)果

4 結(jié) 論

本文提出了一種無需參考信號(hào)的自適應(yīng)對(duì)消法。該方法將迭代自適應(yīng)法引入自適應(yīng)對(duì)消法中，直接利用迭代自適應(yīng)法從輸入信號(hào)中逐個(gè)估計(jì)低頻干擾分量，然后將各個(gè)干擾分量消除。仿真和工程實(shí)驗(yàn)表明，該方法在無需參考信號(hào)的條件下，能夠準(zhǔn)確估計(jì)各個(gè)干擾分量，從而有效地抑制單晶爐熱場(chǎng)溫度測(cè)量值中的干擾，優(yōu)于現(xiàn)有方法，符合單晶爐熱場(chǎng)溫度控制系統(tǒng)的工藝要求。

參考文獻(xiàn)：

[1] Fang H S, Jin Z L, Huang X M. Study and optimization of gas flow and temperature distribution in a Czochralski configuration[J]. Journal of Crystal Growth, 2012, 361: 114-120.

[2] 楊潤(rùn)，李淑娟，彭曦，等. 直拉式單晶爐動(dòng)態(tài)性能分析與預(yù)測(cè)[J]. 西安理工大學(xué)學(xué)報(bào), 2012,28(3):298-302.

Yang R, Li Shujuan, Peng X, et al. Dynamic characteristic analysis and prediction for CZ single crystal furnace[J].Journal of Xi'an University of Technology, 2012,28(3):298-302.

[3] 梁炎明, 劉丁, 趙躍. 基于自適應(yīng)噪聲抵消的 CZ 單晶爐爐膛溫度信號(hào)處理[J]. 控制理論與應(yīng)用, 2011, 28(1): 94-100.

Liang Y M, Liu D, Zhao Y. Temperature signal processing based on adaptive noise cancellation for CZ crystal growing furnace[J].Control Theory & Applications, 2011,28(1): 94-100.

[4] 西蒙·赫金. 自適應(yīng)濾波器原理[M]. 第四版. 鄭寶玉,譯. 北京：電子工業(yè)出版社，2010: 15-17.

[5] McNamara D M, Ziarani A K, Ortmeyer T H. A new technique of measurement of nonstationary harmonics[J]. Power Delivery, IEEE Transactions on, 2007, 22(1): 387-395.

[6] 梁炎明, 劉丁, 任海鵬, 等. 基于自適應(yīng)陷波的直拉單晶直徑跟蹤[J].儀器儀表學(xué)報(bào), 2010, 31(7): 1441-1447.

Liang Y M, Liu D, Ren H P, et al. Diameter tracking in pulling of crystals based on adaptive notch filtering [J]. Chinese Journal of Scientific Instrument, 2010, 31(7): 1441-1447.

[7] 沈福民. 自適應(yīng)信號(hào)處理[M].西安：西安電子科技大學(xué)出版社，2001: 136-147.

[8] Meng J E, Li Z R, Cai H N, et al. Adaptive noise cancellation using enhanced dynamic fuzzy neural networks[J].IEEE Transactions on Fuzzy Systems, 2005, 13(3): 331-341.

[9] Gorriz J M, Ramírez J, Cruces-Alvarez S, et al. A novel LMS algorithm applied to adaptive noise cancellation[J]. IEEE Signal Processing Letters, 2009, 16(1): 34-37.

[10] Allik B, Ilg M, Zurakowski R. Ballistic Roll Estimation using EKF Frequency Tracking and Adaptive Noise Cancellation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2013, 49(4): 2546-2553.

[11] Yardibi T, Li J, Stoica P, et al. Source localization and sensing: A nonparametric iterative adaptive approach based on weighted least squares[J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(1): 425-443.

[12] 湯小為, 朱偉, 湯俊. 基于自適應(yīng)迭代的非參數(shù)化虛擬孔徑合成方法[J].清華大學(xué)學(xué)報(bào): 自然科學(xué)版, 2012, 52(9): 1265-1269.

Tang XW, Zhu W, Tang J. Non-parametric algorithm for virtual aperture synthesis based on the iterative adaptive approach[J].Tsinghua University(Sci & Tech), 2012, 52 (9): 1265-1269.

[13] Ling J, Yardibi T, Su X, et al. Enhanced channel estimation and symbol detection for high speed multi-input multi-output underwater acoustic communications[J].The Journal of the Acoustical Society of America, 2009, 125: 3067-3078.

[14] Roberts W, Stoica P, Li J, et al. Iterative adaptive approaches to MIMO radar imaging[J].IEEE Journal of Selected Topics in Signal Processing, 2010, 4(1): 5-20.

[15] 周杏鵬. 傳感器與檢測(cè)技術(shù)[M].北京：清華大學(xué)出版社, 2010: 24-30.

[16] Stoica P, Li J, He H. Spectral analysis of nonuniformly sampled data: a new approach versus the periodogram[J]. IEEE Transactions on Signal Processing, 2009, 57(3): 843-858.

[17] 張賢達(dá). 矩陣分析與應(yīng)用[M].北京：清華大學(xué)出版社, 2004: 68-69.