梁志國(guó)

(北京長(zhǎng)城計(jì)量測(cè)試技術(shù)研究所 計(jì)量與校準(zhǔn)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100095)

1 引 言

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

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

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

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

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

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

=Ecos(2 π ft+φ)+Q

(1)

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

(2)

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

=Acos(2 π fti+θ)+C

(3)

(4)

(5)

擬合殘差有效值為：

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

2.2 問(wèn)題討論

上述三參數(shù)正弦曲線(xiàn)擬合過(guò)程，是已知信號(hào)頻率f下進(jìn)行的。

特別地，令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è)給定信號(hào)測(cè)量序列的數(shù)字角頻率為w，而不是ω，則，使用1/3個(gè)周期的正弦波形序列，按上述三參數(shù)正弦曲線(xiàn)擬合得如圖1所示歸一化誤差ρ/E與頻率比w/ω關(guān)系曲線(xiàn)波形。圖2為圖1的局部細(xì)化。

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

圖2 p/E≤0.05時(shí)圖1曲線(xiàn)的局部細(xì)節(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代替ω進(jìn)行三參數(shù)擬合。

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

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

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

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

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

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

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

3.1 原理過(guò)程

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

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

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

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

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

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

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

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

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

(7)在fL上執(zhí)行三參數(shù)正弦曲線(xiàn)擬合，獲得AL、θL、CL、ρL；在fR上執(zhí)行三參數(shù)正弦曲線(xiàn)擬合，獲得AR、θR、CR、ρR；在fM上執(zhí)行三參數(shù)正弦曲線(xiàn)擬合，獲得AM、θM、CM、ρM；在fT上執(zhí)行三參數(shù)正弦曲線(xiàn)擬合，獲得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 收斂性問(wèn)題

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

4 仿真實(shí)驗(yàn)驗(yàn)證

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

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

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

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

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

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

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

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

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

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

5 實(shí)驗(yàn)驗(yàn)證

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

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

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

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

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

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

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

6 討 論

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

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

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

7 結(jié) 論

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