曹昌勇

皖西學(xué)院 機械與電子工程學(xué)院，安徽 六安237012

1 引言

1/f 噪聲是一種普遍存在于物理、生物、社會、醫(yī)學(xué)、天文、地理等系統(tǒng)中的自然現(xiàn)象，1925年約翰遜（Johnson）在電子管中首次觀察到1/f噪聲。在以后的研究中，1/f噪聲一直是電子元器件和電子線路低頻噪聲的主要形式[1]。20 世紀80 年代以來的研究發(fā)現(xiàn)，電子元器件中的1/f噪聲的大小與電路的可靠性有密切的關(guān)系。一個1/f 噪聲很大的元器件或電路，很可能存在較嚴重的可靠性問題，其壽命則會較短，而且抗惡劣環(huán)境的能力也就較差。因此，人們通過檢測和分析電子系統(tǒng)1/f噪聲的性質(zhì)，就可以間接地預(yù)測電路或系統(tǒng)的質(zhì)量和可靠性[1-2]。自從發(fā)現(xiàn)1/f 噪聲的近百年時間里，人們對它的產(chǎn)生機理進行了大量的研究，提出了各種不同的模型，但仍然還沒有一個統(tǒng)一的理論。人們發(fā)現(xiàn)1/f 噪聲與系統(tǒng)的缺陷有關(guān)，是一個復(fù)雜的動力系統(tǒng)[3]，滿足長期統(tǒng)計規(guī)律，并且是長程相關(guān)的。對如何生成1/f噪聲的研究，也有不少的模型，比如分數(shù)階布朗運動過程（fBm）；通過混沌模型[4]來生成1/f 噪聲，該方法用一個logistic 混沌序列輸入到hosking 分數(shù)階差分濾波器中，并能產(chǎn)生具有混沌特性的1/f噪聲序列；通過小波變換[5]，該方法根據(jù)白噪聲的小波變換后的系數(shù)仍服從白噪聲規(guī)律，然后正交歸一化小波系數(shù)找到滿足1/f信號生成定理的正交小波系數(shù)集，通過Karhunen-Loeve 展開式生成1/f 信號。對1/f 噪聲的參數(shù)估計，主要采用的辦法有小波變換、最小二乘、時頻分析等；對于跟蹤1/f噪聲，普遍采用自適應(yīng)算法，比如最小均方誤差（LMS）算法，文獻[6]提出了一種變步長LMS 算法跟蹤fBm 過程。由于LMS 算法收斂速度慢，因此本文采用改進的最小遞歸二乘法（RLS）來跟蹤1/f噪聲，但是這種算法存在計算復(fù)雜，迭代步數(shù)增加到一定程度濾波器系數(shù)必定發(fā)散等缺點。

自適應(yīng)算法被廣泛應(yīng)用在系統(tǒng)辨識、預(yù)測建模、噪聲對消以及信道估計等領(lǐng)域。它與濾波器一起構(gòu)成自適應(yīng)濾波器，可以在一定程度上對噪聲進行過濾，基本原理如圖1 所示。

圖1 自適應(yīng)濾波器原理圖

自適應(yīng)算法有很多，用的最多的要數(shù)LMS 和RLS 算法。為了更好地對非平穩(wěn)信號進行處理，需要對算法進行改進。人們在改進這兩種算法的過程中做了大量的工作。對LMS 算法的改進一般采用的辦法是變步長[6-7]，或者是考慮引入梯度算法[8]或其他算法實現(xiàn)對快速變化的信號進行跟蹤。對于RLS 算法來說，它的收斂速度比LMS 要快得多，但它是以增加計算量為代價的。影響計算復(fù)雜度的關(guān)鍵在于該濾波器算法的結(jié)構(gòu)上。文獻[9]提出了一種分級實現(xiàn)的RLS 算法，它把大的濾波器階數(shù)分解成幾個級聯(lián)的階數(shù)小的濾波器；文獻[10]通過輸入濾波器的信號的自相關(guān)函數(shù)是服從遞歸規(guī)律建立了擴展算法（ERLS）；文獻[11-12]通過誤差反饋的辦法，在文獻[9]的基礎(chǔ)上增加一個誤差判斷開關(guān)，通過開關(guān)切換去修改變化比較大的級聯(lián)濾波器的系數(shù)，達到部分更新濾波器系數(shù)的效果，進而減少運算量。對RLS 算法的收斂速度上的改進主要是實時修正RLS 算法的記憶因子[13]。對于非平穩(wěn)信號來說，主要是期望能快速根據(jù)輸入信號的變化來調(diào)整濾波器系數(shù)。

正規(guī)或者快速RLS 算法是建立在迭代計算基礎(chǔ)上的，因此在有限的計算機精度的條件下，必定要發(fā)散。為了解決濾波器權(quán)系數(shù)發(fā)散的問題以及快速準確的跟蹤信號變化的問題，本文根據(jù)輸入濾波器的信號的自相關(guān)矩陣應(yīng)該滿足慢變化的規(guī)律，引入誤差的非線性函數(shù)來控制這個變化以實現(xiàn)實時快速跟蹤1/f信號。該算法繼承了正規(guī)RLS算法的快收斂特性，并且不增加正規(guī)的RLS 算法的計算量，誤差也比正規(guī)RLS 算法有所減小。通過該方法對混沌1/f噪聲過程與fBm 過程進行跟蹤仿真，實驗表明它具有良好的穩(wěn)定性，算法收斂速度較正規(guī)RLS 算法沒有降低，在小的濾波器階數(shù)的條件下也具有良好的跟蹤能力。該算法在遇到突變輸入時會快速調(diào)整濾波器系數(shù)，因此在短時間內(nèi)濾波器系數(shù)的改變比正規(guī)RLS 算法要快得多，這就實現(xiàn)了快速跟蹤目標，從而減小了誤差。

RLS 算法中記憶因子對算法的跟蹤速度有比較大的影響。當(dāng)記憶因子比較小時，算法響應(yīng)快變信號的能力強，但抗干擾能力弱，適合于跟蹤非平穩(wěn)信號；當(dāng)記憶因子較大時，算法響應(yīng)快變信號能力比較差，但抗干擾能力強，適合跟蹤平穩(wěn)信號。記憶因子的大小與抗干擾強弱是矛盾的。因此記憶因子的調(diào)整量太大和太小都不好，本文給出了記憶因子與輸入信號自相關(guān)矩陣的特征值（fBm 噪聲信號）的一個表達式，為動態(tài)調(diào)整記憶因子提供了理論依據(jù)。

2 1/f噪聲的兩個模型

所謂1/f過程是指隨機過程x（t）的功率譜是冪函數(shù)的形式，即S(f)=σ2/f a，其中σ2是常數(shù)，由于其產(chǎn)生機理不明，一般都把它當(dāng)做是一種隨機過程，它在低頻部分的功率特別大，特別是頻率靠近0 的時刻。

2.1 1/f混沌模型

文獻[4]通過采用logistic 方程生成混沌序列，然后把混沌序列輸入到hosking 分數(shù)階濾波器，表示如下。

其中，H(Z)為hosking 分數(shù)階濾波器的傳遞函數(shù)，其沖擊響應(yīng)為：

其中，u(n)是單位階躍函數(shù)，式（1b）的離散遞推表達為：

通過求式（1）與式（1c）的卷積即可得到混沌1/f噪聲信號。式（1）產(chǎn)生的混沌序列的均值為0，自相關(guān)函數(shù)與白噪聲的自相關(guān)函數(shù)相似。

2.2 分數(shù)維布朗運動過程

文獻[6]提出了一種變步長的LMS 算法跟蹤fBm 過程。如果一個時間連續(xù)的隨機過程的統(tǒng)計特性是尺度不變的就叫自相似過程，自相似可用公式表示如下：

z(t)即為一個具有自相似指數(shù)H（Hurst 指數(shù)）的自相似隨機過程，它對任意的標量c＞0 都具有自相似性。式（2）中的不等號僅僅在統(tǒng)計意義上對所有有限分布滿足自相似性。具有平穩(wěn)增加的非平穩(wěn)的，且具有自相似性的高斯隨機過程就是分數(shù)維布朗運動過程，計為BH(t)，一個fBm 過程的一階導(dǎo)數(shù)稱為分數(shù)維高斯噪聲（fGn）。離散的fBm 可表示如下：

Ts為采樣周期，因為分數(shù)維布朗運動具有自相似性，不失一般性，取Ts=1。離散fBm 的均值、方差與自相關(guān)函數(shù)表示如下：

引理1 Hurst 指數(shù)為0 ＜H ＜1 的離散第一階分數(shù)維布朗運動（1-fBm）的自相關(guān)函數(shù)R(n)（M×M 階矩陣）能被對角化為R(n)=QΛ(n)QT，如果它的近似估計在n 比較大時成立，這里Q(n)≈Q 是一個常數(shù)正交矩陣，Λ(n)=diag{λ1,λ2,…,λM-1,λM(n)}。 R(n)的前面M-1個特征值是與時間無關(guān)的，特征值λM(n)是時間n 的函數(shù)。所有特征值都與Hurst指數(shù)H 相關(guān)，指數(shù)H 刻畫了離散分數(shù)維布朗運動，具體見公式（2a）。其最大的特征值是時間n 和Hurst指數(shù)H 的函數(shù)。該特征值可用下式表示：

3 RLS 算法及修正

3.1 基本RLS 算法

假如濾波器階數(shù)為M，n 時刻的權(quán)系數(shù)為W(n)=[W1(n),W2(n),…,WM(n)]T，輸入數(shù)據(jù)為向量U(n)=[u(n),u(n-1),…,u(n-M+1)]T，則輸出設(shè)為y(n)其表達如下：

期望輸出為d(n)，表達如下：

其中W0為最優(yōu)濾波器系數(shù)，假設(shè)η(n)是均值為0，方差為的高斯白噪聲，文獻[14]考慮了η(n)服從獨立同分布的噪聲，采用類A 噪聲模型進行分析，本文僅考慮高斯白噪聲的情況。n 時刻的估計誤差為：

基于最小二乘法的RLS 算法的最小化代價函數(shù)如下：

λ 為遞歸最小二乘算法的記憶因子，滿足0 ＜λ ＜1，它對當(dāng)前的誤差記憶強，對過去時刻的誤差記憶弱。把式（5）代入式（6）并兩邊對W 求導(dǎo)數(shù)，令其導(dǎo)數(shù)為0，即

對式（9）利用矩陣求逆引理可得到其逆矩陣的遞推表達式：

令P(n)=R-1(n)為輸入信號的自相關(guān)矩陣的逆矩陣：

將式（11）與式（12）代入到式（8）得到權(quán)系數(shù)的遞推公式（13b）：

上面的四個表達式給出了常規(guī)的RLS 算法的一個迭代表達，要滿足迭代，需要初始化的項為P(n)和W(n)。一般W(0)初始化為長度為M 的0 向量，P(0)=δ-1IM，IM為M 階的單位矩陣，δ 當(dāng)輸入信號的信噪比大時取小的正數(shù)，反之取大的正常數(shù)[15]。

在式（12）中，P(n)是輸入向量的時間平均相關(guān)陣的逆矩陣。

3.2 修正RLS 算法

本文對式（13c）P(n)進行修正，用一個非線性的誤差函數(shù)來修正它，并用式（15a）取代式（13c）。

公式（14）所對應(yīng)的關(guān)系曲線如圖2 所示，其中a 取0.8，b 取3，可以根據(jù)程序設(shè)置不同的值來討論曲線的變化情況；公式（15）不需要贅畫，因為IM是一個單位矩陣。 f(e(n))是關(guān)于誤差e(n)的非負非線性函數(shù)，把絕對值與a 去掉即為參考文獻[3]中函數(shù)的形式。式（14）中的a 與b 是正數(shù)，一般與被跟蹤的信號有關(guān)，可以通過圖2 曲線圖討論來確定。誤差一般不為0，當(dāng)誤差增大時P(n)增大，即表示信號在快速變化，信號發(fā)生了突然的變化，這需要在下一次比較大的調(diào)整權(quán)系數(shù)W(n)，通過觀察式（13），由于P(n)增加，導(dǎo)致k(n)增加比較大，從而W(n)得到了大的調(diào)整，說明信號是快變的；相反，當(dāng)誤差很小時，表示濾波器權(quán)系數(shù)不需要進行大的調(diào)整，從式（14）可見小的誤差不會對原來的P(n) 產(chǎn)生大的改變，因此k(n) 的變化也比較小，從而W(n)的調(diào)整很小，這表明輸入信號在這一小段時間內(nèi)是平穩(wěn)的。另一方面，由于P(n)的改變，可以說是對它進行了一個變化的初始化工作。

圖2 F（n）同e（n）的關(guān)系曲線

對于RLS 算法，P(n)是集平均輸入向量自相關(guān)陣R(n)在n 時刻的逆矩陣，在正規(guī)的RLS 算法中，該項是一直迭代下去的，由于計算機精度問題，必定會導(dǎo)致P(n)不可逆，導(dǎo)致R(n)的特征值發(fā)生突變。特征值發(fā)生突變的根本原因是有部分非平穩(wěn)數(shù)據(jù)的輸入，造成誤差增大，導(dǎo)致權(quán)系數(shù)調(diào)整過大，最終算法發(fā)散。因此一種辦法是通過關(guān)于誤差的非線性函數(shù)來修正權(quán)系數(shù)的調(diào)整[14]，即非線性RLS 算法。本文用非線性函數(shù)來微量調(diào)整P(n)矩陣，這樣可以避免特征值的突變，但算法仍然是線性的RLS 算法。

3.3 與λ 相關(guān)的一個表達式

下面推導(dǎo)記憶因子λ 與輸入信號（1-fBm）的自相關(guān)矩陣的特征值的一個重要關(guān)系。

把式（4）、式（16）代入式（13a）得到：

由式（12）、式（13b）、式（16）和（17）推出：

對式（18）兩邊分別取期望得到：

根據(jù)P(n)滿足慢變規(guī)律，可知它相對于輸入U(n)與誤差e(n)來說，它的變化顯得很小，因此可以把它看做是與輸入和誤差相互獨立的，故可表達如下：

由式（16）和式（13b）知，V(n)與U(n)獨立，假設(shè)觀測噪聲η(n)與U(n)相互統(tǒng)計獨立，且η(n)的均值為0，將式（17）代入到式（20）得到：

其中，RU(n)=E[U(n)UT(n)]為n 時刻輸入濾波器數(shù)據(jù)的自相關(guān)函數(shù)，在式（21）中，根據(jù)式（12）得到，且RLS 算法需要初始化P(0)=δ-1IM，即R(0)=δIM，故：

把式（22）、式（23）代入式（21），并考慮到輸入數(shù)據(jù)的自相關(guān)矩陣是緩變的，把當(dāng)前時刻自相關(guān)陣與之前時刻的相關(guān)陣看做是可以近似統(tǒng)一的（記憶因子的存在也為這種近似統(tǒng)一提供了依據(jù)），因此得到：

由引理1 可知，式（24）RU(n)=QΛ(n)QT，其中Q 為正交的矩陣，Λ(n)是由RU(n)的特征值構(gòu)成的對角陣，設(shè)其特征值為λ1,λ2,…,λM-1,λM(n)，其中λM(n)在式（2c）中給出。因此式（24）轉(zhuǎn)化為：

設(shè)B=I-Λ2(n)-1，由于，因此B 的特征值可以用Λ(n)的特征值來表示。式（26）給出了一個迭代過程。RLS 算法收斂意味著當(dāng)n 趨于無窮時：

即當(dāng)n 趨于無窮時，W(n+1)→W0達到最優(yōu)值。

其中λΛ2(j,i)為Λ2在第i 次迭代過程中第j 個特征值。由可以把Λ2表示出來：

最后由式（28）得到收斂的條件為：

公式（30）給出了記憶因子λ 與輸入信號的自相關(guān)矩陣的特征值λΛ(j,n)在n 時刻的一個表達式，這為在跟蹤過程中動態(tài)調(diào)整RLS 算法的記憶因子λ 給出了理論上的依據(jù)。

公式（30）說明調(diào)整過程中算法的記憶因子是受到自相關(guān)矩陣的特征值λΛ(j,n)控制的。該公式是記憶因子滿足的一個關(guān)系表達式，可以根據(jù)這個公式來設(shè)置記憶因子，當(dāng)然，如果記憶因子在每次迭代過程中的設(shè)置不滿足這個條件時，算法肯定會發(fā)散，從這個角度看，屬于收斂條件范疇的討論。

4 仿真結(jié)果

4.1 對混沌1/f信號跟蹤的仿真

根據(jù)第2章中的第一個模型產(chǎn)生1/f時間序列（初始值取0.3，hosking 分數(shù)階濾波器d=0.4），取3 000 個離散數(shù)據(jù)。設(shè)修正的RLS 算法濾波器階數(shù)為M=10，觀測噪聲方差為0.4，得到如圖3 的跟蹤結(jié)果（誤差均值與方差分別為0.046 5、0.467 4）。

圖3 修正RLS 算法對混沌1/f噪聲序列的跟蹤（a=5，b=0.005）

標準RLS 算法的跟蹤結(jié)果如圖4 所示（誤差均值為0.047 7，方差為0.480 4）。

圖4 標準RLS 算法跟蹤混沌1/f噪聲信號

從對混沌1/f噪聲信號的跟蹤情況看，修正RLS 算法比標準RLS 算法弱占優(yōu)勢，其均值和方差均比標準算法小。

4.2 對1-fBm 噪聲信號的跟蹤

取H=0.7，濾波器階數(shù)為10，修正算法與標準算法結(jié)果如圖5（誤差均值與方差分別為0.001 7、0.104 0）和圖6 所示（均值為-0.005 6，方差為0.274 2）。

圖5 修正RLS 算法對混沌1/f噪聲序列的跟蹤（a=5，b=0.5）

圖6 標準算法跟蹤1-fBm 信號

從對混沌1/f 信號或1-fBm 的結(jié)果看，修正的算法比標準算法要好。從圖7 與圖8 在對跟蹤1-fBm 過程中濾波器系數(shù)的調(diào)整過程可以看出，修正算法濾波器的系數(shù)可以發(fā)生快速改變，而標準算法只是很慢地變化，這也導(dǎo)致它的誤差比較大，跟蹤速度達不到的原因。

圖7 修正算法對濾波器系數(shù)的調(diào)整過程

圖8 標準RLS 算法對濾波器系數(shù)的調(diào)整過程

圖7、圖8 中不同的顏色曲線表示的是RLS 算法中各個濾波器系數(shù)的動態(tài)曲線，由于濾波器系數(shù)一般都不止一個，因此為了區(qū)別每一個系數(shù)的動態(tài)變化，作圖時采用了不同的顏色。

從圖7 可見，迭代時會出現(xiàn)濾波器系數(shù)的調(diào)整過程；對比圖8 的系數(shù)調(diào)整可見，算法不會影響收斂速度。

5 結(jié)束語

RLS 算法是一種快速收斂的算法，但是在跟蹤非平穩(wěn)信號時無法快速調(diào)整濾波器系數(shù)以便實時跟蹤，并且在迭代次數(shù)達到很大時算法必將發(fā)散，本文通過引入一個非線性函數(shù)來對輸入信號相關(guān)逆矩陣進行修正，這樣可以快速調(diào)整濾波器系數(shù)，以適應(yīng)信號的快速變化。對兩個1/f噪聲模型產(chǎn)生的序列進行仿真實驗，結(jié)果表明修正方法比正規(guī)的RLS 算法要好。此外，本文還推導(dǎo)了RLS 算法的記憶因子與輸入信號自相關(guān)矩陣特征值的一個關(guān)系，這為動態(tài)調(diào)整記憶因子提供了理論上的依據(jù)。

[1] 莊奕琪，馬中發(fā)，杜磊.1/f 噪聲之謎[J].世界科技研究與發(fā)展，1999，21（4）：69-72.

[2] 楊冬云，王司.1/f 分形噪聲理論及其在信號處理中的應(yīng)用研究綜述[J].黑龍江工程學(xué)院學(xué)報：自然科學(xué)版，2004（3）：31-35.

[3] 許生龍.1/f噪聲的動力學(xué)統(tǒng)計理論[J].紅外技術(shù)，2003（5）：63-67.

[4] 張強，馬潤年，許進.1/f過程的混沌模型[J].信號處理，2001（5）：392-394.

[5] 李光，于盛林.基于小波變換的擬1/f 信號生成[J].數(shù)據(jù)采集與處理，2003（1）：119-122.

[6] Gupta A，Joshi S D.A novel least mean squares algorithm for tracking a discrete-time fBm process[C]//Proceedings of 2006 Annual India Conference，2006.

[7] 覃景繁，歐陽景正.一種新的變步長LMS 自應(yīng)用濾波算法[J].數(shù)據(jù)采集與處理，1997（3）：171-174.

[8] Ang W P，F(xiàn)arhang-Boroujeny B.A new class of gradient adaptive step-size LMS algorithms[J].IEEE Trans on Signal Processing，2001，49：805-809.

[9] Wu An-Yeu，Liu K J R.Split recursive least-squares-algorithms，architectures，and applications[J].IEEE Trans on Circuits and Systems Part II：Analog and Digital Signal Processing，1996，43：645-658.

[10] Haykin S，Sayed A H，Zeidler J R，et al.Adaptive tracking of linear time-variant systems by extended RLS algorithms[J].IEEE Transactions on Signal Processing，1997，45（5）.

[11] Jin Minglia.Partial updating RLS algorithm[C]//Proceedings of the 7th International Conference on Signal Processing（ICSP’04），Beijing，China，2004：392-395.

[12] Niu Wei，Wang Minxi，Chen Kaya.A novel reduced order RLS predistortion[C]//Proceedings of Microwave Conference，Asia-Pacific，2005.

[13] Park D，Jun B E，Kim J H.Fast tracking RLS algorithm using novel variable forgetting factor with unity zone[J].Electronics Letters，1991，27：2150-2151.

[14] Weng J F，Leung S H.Nonlinear RLS algorithm for amplitude estimation in class A noise[J].IEEE Proceedings Communications，2000，147：81-86.

[15] Wei P C，Han Jun，Zeidler J R，et al.Comparative tracking performance of the LMS and RLS algorithms for chirped narrowband signal recovery[J].IEEE Transactions on Signal Processing，2002，50：1602-1609.

[16] Chun B，Kim B，Lee Y H.Generalization of exponentially weighted RLS algorithm based on a state-space model[C]//Proceedings of the IEEE International Symposium on Circuits and Systems.[S.l.]：IEEE Press，1998：198-201.

[17] Yeo H K，Sharif B S，Hinton O R，et al.Improved RLS algorithm for time-variant underwater acoustic communications[J].Electronics Letters，2000，36：191-192.