非均勻采樣條件下殘周期正弦波形的最小二乘擬合算法

梁志國

(北京長城計量測試技術(shù)研究所 計量與校準技術(shù)重點實驗室，北京 100095)

1 引 言

正弦波形四參數(shù)曲線擬合是測量測試中的基本算法，具有廣泛的應用前景和工程價值，因而，有眾多學者對該問題進行了廣泛而深入的研究[1～25]。而殘周期正弦波曲線擬合，則在超低頻參數(shù)的快速估計、調(diào)制解調(diào)等方面具有特別的價值和意義，其最主要的困難在于殘周期條件下，正弦波擬合參數(shù)的初始值很難獲得，因而在四參數(shù)非線性迭代的傳統(tǒng)算法中不易實現(xiàn)最小二乘擬合。直到文獻[11]采用了三參數(shù)擬合結(jié)合一維搜索迭代方式解決了該問題，才最終獲得了殘周期均勻采樣條件下的最小二乘擬合結(jié)果。

本文主要內(nèi)容，是針對采樣間隔不恒定的非均勻采樣情況下，殘周期正弦波形的四參數(shù)曲線擬合開展研究，以解決非均勻采樣條件下殘周期正弦波形的曲線擬合問題。它將可以用于均勻采樣、非均勻采樣、隨機采樣等同時提供時間坐標和幅度信息的采樣波形序列的正弦波形擬合。

2 三參數(shù)正弦波形的最小二乘擬合算法

2.1 三參數(shù)正弦曲線擬合法

設(shè)理想正弦波形為：

y(t)=E1cos(2 π ft)+E2sin(2 π ft)+Q

=Ecos(2 π ft+φ)+Q

(1)

數(shù)據(jù)記錄序列為已知時刻t1,t2,…,tn的采集樣本y1,y2,…,yn。三參數(shù)正弦曲線擬合過程，即為輸入信號的頻率f已知，選取或?qū)ふ褹1、B1、C，使下式所述殘差平方和ε最?。?/p>

(2)

由式(2)的ε最小，可得擬合函數(shù)：

=Acos(2 π fti+θ)+C

(3)

(4)

(5)

擬合殘差有效值為：

由于這是一種閉合算法，因而收斂是肯定的。

2.2 問題討論

上述三參數(shù)正弦曲線擬合過程，是已知信號頻率f下進行的。

特別地，令tmin=min{ti}；tmax=max{ti}；則，平均采樣速率v=(n-1)/(tmax-tmin)；數(shù)字角頻率 ω=2 π f/v。

在均勻采樣條件下，式(1)可表示成下列離散形式：

y(i)=E1cos(ω·i)+E2sin(ω·i)+Q

=Ecos(ω·i+φ)+Q

(6)

則擬合函數(shù)將是：

=Acos(ω·i+θ)+C

(7)

不失一般性，設(shè)給定信號測量序列的數(shù)字角頻率為w，而不是ω，則，使用1/3個周期的正弦波形序列，按上述三參數(shù)正弦曲線擬合得如圖1所示歸一化誤差ρ/E與頻率比w/ω關(guān)系曲線波形。圖2為圖1的局部細化。

圖1 1/3周期波形擬合時歸一化誤差ρ/E與頻率比w/ρ關(guān)系曲線Fig.1 Variety of error ρ/E via frequency ratio w/ω when 1/3 period

圖2 p/E≤0.05時圖1曲線的局部細節(jié)Fig.2 Details of Fig.1 in the range of p/E≤0.05

其中ω=2 π/10 000；E=4；Q=0；φ分別取0°、35°、70°、105°、140°、175°、210°、245°、280°、315°、350°；以w代替ω進行三參數(shù)擬合。

從圖1及其局部細化曲線波形圖2可見，對于信號頻率ω的不足一個信號周期的殘周期正弦曲線波形擬合，有如下規(guī)律：

(1) 三參數(shù)最小二乘擬合法在(0,2ω]范圍內(nèi)ρ/E極小值存在且唯一，極值點就是ω頻率點。

(2) 幅度和相位的變化基本不影響ρ/E的變化規(guī)律；但初始相位φ的變化將影響ρ/E的幅度值。

(3) 當使用的擬合頻率w大于信號頻率的2倍時，ρ/E的量值均普遍大于擬合頻率w小于信號頻率的情況。該規(guī)律可用于判定收斂區(qū)間上界，而收斂區(qū)間的下界則可由足夠接近0頻的一個微小頻率值代替。

(4) 可在(0,2ω]范圍內(nèi)通過對ω的一維搜索找出ρ/E的極小值點。

使用其它部分周期波形曲線進行擬合可以獲得與上述規(guī)律相同的結(jié)論。

3 四參數(shù)正弦波形的最小二乘擬合算法

3.1 原理過程

對上述三參數(shù)正弦曲線擬合方法的改造，可獲得一種絕對收斂的殘周期正弦波形四參數(shù)曲線擬合方法。

假設(shè)平均采集速率為v，待估計的正弦波頻率目標值為f0，待估計的正弦波采樣序列所含信號不足一個周期，個數(shù)為p(0q/τ；因而，可以肯定f0∈[q/τ,2/τ]，在區(qū)間[q/τ,2/τ]內(nèi)的任意頻率f下，殘差平方和ε(f)的極值存在且唯一！這樣，便將四參數(shù)正弦波曲線擬合中，對幅度、頻率、相位、直流分量4個參數(shù)的四維非線性搜索，變成了對頻率分量f造成的ε(f)的一維線性搜索，可保證在區(qū)間[q/τ,2/τ]內(nèi)，用三參數(shù)擬合法實現(xiàn)的四參數(shù)正弦曲線擬合過程收斂。該四參數(shù)擬合過程如下：

(1)設(shè)定擬合迭代停止條件為he。

(2)從已知時刻t1,t2,…,tn的正弦波采集樣本y1,y2,…,yn。使用計點法獲得信號波形占用時間長度為τ=tn-t1；平均采集速率v=(n-1)/(tn-t1)，選取因子q(例如q=10-5)，確定目標頻率f0的存在區(qū)間[q/τ,2/τ]。

(3)確定迭代左邊界頻率：fL=q/τ；迭代右邊界頻率：fR=2/τ。

(4)令中值頻率：fM=(fR+fL)/2；在左邊界頻率和右邊界頻率和中值頻率上分別利用3參數(shù)擬合公式計算各自的擬合殘差ρ(fL)、ρ(fM)和ρ(fR)。

(5)判斷是否ρ(fL)<η·ρ(fM)？其中，η為判據(jù)因子，取值范圍為1～1.5。

若ρ(fL)<η·ρ(fM)，則令fR=fM，fL=不變，重復執(zhí)行(4)～(5)的過程。

(6)若ρ(fL)≥η·ρ(fM)，則必有fR<2f0，確定迭代左邊界頻率為fL；迭代右邊界頻率fR；中值頻率：fM=fL+0.618×(fR-fL)；fT=fR-0.618×(fR-fL)。

(7)在fL上執(zhí)行三參數(shù)正弦曲線擬合，獲得AL、θL、CL、ρL；在fR上執(zhí)行三參數(shù)正弦曲線擬合，獲得AR、θR、CR、ρR；在fM上執(zhí)行三參數(shù)正弦曲線擬合，獲得AM、θM、CM、ρM；在fT上執(zhí)行三參數(shù)正弦曲線擬合，獲得AT、θT、CT、ρT。

(8)若ρM<ρT，則ρ=ρM，有f0∈[fT,fR]，fL=fT，fT=fM；fM=fL+0.618×(fR-fL)；若ρM>ρT，則ρ=ρT，有f0∈[fL,fM]，fR=fM，fM=fT；fT=fR-0.618×(fR-fL)。

(9)判定是否|(ρM(k)-ρT(k))/ρT(k)|

3.2 收斂性問題

四參數(shù)正弦曲線擬合是一個迭代過程，也有收斂性問題。如上所述，對于含p(p?1)個信號周期的測量序列的情況，如圖1所示，四參數(shù)正弦曲線擬合的收斂區(qū)間為(0,2ω0]；其中，ω0為信號的真實數(shù)字角頻率。

4 仿真實驗驗證

選取測量范圍-5～+5 V，幅度4 V，直流分量0 V，初始相位100°，頻率6 Hz，均勻采樣速率8 kSa/s，采樣間隔τ0=0.125 ms，采樣序列長度n=400的采樣序列，約含0.3個波形周期。使用本文上述四參數(shù)擬合算法進行正弦參數(shù)估計，獲得擬合曲線與采樣曲線波形如圖3所示，兩者之差異曲線如圖4所示。其擬合參數(shù)如表1所示。

表1 殘周期正弦序列參數(shù)擬合結(jié)果Tab.1 Results of partial period sine wave curve-fitting

圖3 殘周期擬合曲線與均勻采樣曲線yiFig.3 Uniform sampling points yi and curve-fitting

圖4 殘周期擬合曲線與均勻采樣曲線差異Fig.4 difference between uniform sampling series and

令正弦波模型參數(shù)不變，使用隨機采樣間隔進行序列采樣，獲得采樣序列并按照本文上述方法進行正弦波曲線擬合，得采樣序列曲線及其擬合曲線分別如圖5～圖8所示。其中，圖5為波形長度約為0.4個周期的非均勻采樣序列曲線及其擬合曲線，其采樣間隔為(RND-0.4)×20τ0。其中，RND為在[0,1]內(nèi)均勻分布的隨機數(shù)?？梢娪捎诓蓸娱g隔呈非均勻隨機變化狀態(tài)，按照等間隔均勻采樣模式繪制的曲線圖5，其中的“正弦波形”已經(jīng)變化非常大。

圖6為擬合曲線與采樣序列之差異曲線?？梢妰烧邤M合得非常好。相應擬合參數(shù)如表1所示。由表1數(shù)據(jù)可見，非均勻采樣條件下正弦曲線擬合誤差要略大于均勻采樣條件下，但差異不太大，造成差異的原因有待于進一步研究。圖7為波形長度約為1.2個周期的非均勻采樣序列曲線及其擬合曲線，其采樣間隔為(RND-0.3)×20τ0。

圖6 殘周期擬合曲線與非均勻采樣曲線差異Fig.6 Difference between non-uniform sampling series and

圖7 擬合曲線與非均勻采樣曲線yiFig.7 Non-uniform sampling points yi and curve-fitting

圖8為擬合曲線與采樣序列之差異曲線，擬合情況良好，相應擬合參數(shù)如表1所示。

圖8 擬合曲線與非均勻采樣曲線差異Fig.8 Difference between non-uniform sampling series and

5 實驗驗證

用超低頻振動標準裝置作激勵源，給出位移幅度36.32 mm的正弦振動，頻率0.040 000 Hz，用ASQ-1CA型位移傳感器進行測量，以NI PXI-6281型數(shù)據(jù)采集系統(tǒng)執(zhí)行采集，其A/D位數(shù)18 Bit，工作量程±2.5 V，采樣速率200 Sa/s，采集數(shù)據(jù)個數(shù) 5 000 點。為含有約1個周期波形的振動序列。

使用上述四參數(shù)正弦波擬合方法獲得其振動波形幅度為1 771.02 mV，頻率為0.040 56 Hz，相位為230.098°，直流分量為120.53 mV。殘差有效值ρ=30.157 mV。波形總失真度為2.4%。其測量序列波形及擬合曲線如圖9所示，測量序列波形與擬合曲線之差如圖10所示。

圖9 擬合曲線與均勻采樣曲線yiFig.9 Uniform sampling points yi and curve-fitting

圖10 擬合曲線與均勻采樣曲線差異Fig.10 Difference between uniform sampling series and

將上述測量曲線波形截取3段后，變成非均勻采樣序列，使用上述四參數(shù)正弦波擬合方法獲得其振動波形幅度為1 746.93 mV，頻率為0.040 42 Hz，相位為231.354°，直流分量為123.94 mV。殘差有效值ρ=34.154 mV，波形總失真度為2.8%。其測量序列波形及擬合曲線如圖11所示，測量序列波形與擬合曲線之差如圖12所示。

圖11 擬合曲線與非均勻采樣曲線yiFig.11 Non-uniform sampling points yi and curve-fitting

圖12 擬合曲線與非均勻采樣曲線差異Fig.12 Difference between non-uniform sampling series and

6 討 論

從上述仿真和實際實驗結(jié)果可見，本文所述基于非均勻采樣條件的殘周期正弦波擬合的測量方法，可以用于不足一個波形周期的殘周期非均勻采樣條件下正弦波形參數(shù)的測量估計，并可以給出幅度、頻率、相位、直流分量等基本信息參量。目標是處理殘周期條件下的正弦波擬合，實際上可以在2個波形周期以下的情況均可直接使用本文方法處理。其最大的特點是針對殘周期正弦波模型，不論其采樣間隔是否均勻，均可以獲得有效收斂結(jié)果，特別是針對均勻采樣序列剔出一段異常波形情況，在實際工作中時有發(fā)生，而本來希望獲得均勻采樣序列，但由于周期過長以及測量過程不完善等因素，使得采樣序列變成非均勻采樣序列的情況也比比皆是。對于隨機采樣情況，甚至時間排序出現(xiàn)前后顛倒的情況下，本文方法依然可以有效獲得擬合結(jié)果。體現(xiàn)出算法良好的收斂性和魯棒性，可直接用于解決上述范疇內(nèi)的工程技術(shù)問題。

仿真實驗表明，與均勻采樣條件相比，非均勻采樣序列正弦擬合誤差要略大一些，在本文實驗中，幅度差異約在0.2%，頻率差異約在0.5%，相位差異約在1.5°，直流分量差異在0.3%(相對于正弦幅度)。實際實驗情況也驗證了這些。

此前的研究已經(jīng)表明，與多周期情況下的測量參數(shù)相比，殘周期參數(shù)擬合存在更大的誤差，特別是有較大直流分量情況下。由于信息不全，波形失真、序列樣本長度、波形周期長度等眾多因素都將影響參數(shù)估計效果，因而對于波形序列長度要求也是不一樣的，研究表明，在通常的波形質(zhì)量條件下，1/5以上周期的波形長度即可以進行參數(shù)的正確估計。

7 結(jié) 論

綜上所述可見，本文主要是提出了非均勻采樣條件下殘周期正弦波的四參數(shù)擬合的一種方法，給出了詳細過程和收斂區(qū)間。分別在隨機間隔采樣和非隨機間隔采樣兩種條件下給出了比較結(jié)果，結(jié)果表明，本文所述方法具有良好的收斂性和魯棒性，僅要求已知采樣序列和每個采樣點的時刻，沒有任何其它先決條件要求，是一種表現(xiàn)良好的正弦參數(shù)擬合方法，可直接用于非均勻采樣條件下殘周期正弦波形的參數(shù)擬合，對于超低頻參數(shù)估計尤其具有特別的意義和價值。