梁志國(guó)， 朱振宇

（1.北京長(zhǎng)城計(jì)量測(cè)試技術(shù)研究所計(jì)量與校準(zhǔn)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京100095；2.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072）

非均勻采樣正弦波形的最小二乘擬合算法

梁志國(guó)1， 朱振宇2

（1.北京長(zhǎng)城計(jì)量測(cè)試技術(shù)研究所計(jì)量與校準(zhǔn)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京100095；
2.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072）

針對(duì)非均勻采樣序列提出了一種四參數(shù)正弦波形曲線(xiàn)擬合方法。它借助于頻率已知的三參數(shù)正弦曲線(xiàn)擬合算法實(shí)現(xiàn)，既可以在非均勻采樣條件下獲得正弦波采樣序列的4個(gè)參數(shù)擬合值，也可用于均勻采樣條件下正弦波形采樣序列的曲線(xiàn)參數(shù)擬合。其特點(diǎn)是將幅度、頻率、相位和直流分量4個(gè)參數(shù)的四維非線(xiàn)性最優(yōu)化問(wèn)題轉(zhuǎn)化成頻率參數(shù)的一維線(xiàn)性搜索和最優(yōu)化問(wèn)題。其優(yōu)點(diǎn)是無(wú)需先驗(yàn)初值估計(jì)，具有良好的收斂性、魯棒性，以及明確的收斂區(qū)間。仿真及實(shí)驗(yàn)均證明了該方法的有效性和可行性。該方法可用于解決非均勻采樣和均勻采樣下正弦波測(cè)量序列的參數(shù)擬合問(wèn)題。

計(jì)量學(xué)；非均勻采樣；正弦波；曲線(xiàn)擬合；參數(shù)估計(jì)；校準(zhǔn)

1 引 言

四參數(shù)正弦波形最小二乘擬合算法是具有廣泛應(yīng)用價(jià)值和前景的基本算法，因而吸引了眾多學(xué)者對(duì)該問(wèn)題進(jìn)行了廣泛而深入的研究。這些研究中，多數(shù)是針對(duì)具有恒定采樣間隔的均勻采樣方式進(jìn)行的，這也是多數(shù)數(shù)字信號(hào)處理理論和方法中的基本假設(shè)和前提［1～13］。然而，伴隨著高速寬帶信號(hào)的測(cè)量處理需求，非均勻采樣技術(shù)開(kāi)始出現(xiàn)。例如，以多A／D采樣序列合成更高速采樣序列的高速采樣技術(shù)，以針對(duì)周期信號(hào)波形的取樣示波器的延遲采樣技術(shù)和隨機(jī)采樣技術(shù)等，盡管它們最終目的都是合成等效的均勻采樣序列，但由于技術(shù)的不完善等因素，產(chǎn)生了非均勻采樣效應(yīng)。若仍然按照均勻采樣方式進(jìn)行處理，則將造成時(shí)基失真和頻譜混疊。多年以來(lái)，有眾多研究是專(zhuān)門(mén)針對(duì)這些非均勻采樣導(dǎo)致的時(shí)基失真開(kāi)始的。

另有一類(lèi)現(xiàn)象，是在均勻采樣序列中，由于瞬態(tài)干擾或短期粗大誤差造成了測(cè)量序列中有部分波形嚴(yán)重失真和錯(cuò)誤，以往人們處理的方式是重新測(cè)量或直接切除部分波形。人們也希望找到一種方法，在直接剔除嚴(yán)重錯(cuò)誤部分后形成的“非均勻采樣”序列波形仍然能夠進(jìn)行整體有效的參數(shù)擬合。若將它們完全按照均勻采樣方式處理，將會(huì)造成額外的誤差或波形失真，根本無(wú)法進(jìn)行有效的參數(shù)擬合與估計(jì)。

由于正弦波形四參數(shù)擬合是一個(gè)對(duì)初始值準(zhǔn)確度要求較高的非線(xiàn)性迭代過(guò)程，因而以往的四參數(shù)正弦波曲線(xiàn)擬合大都要求首先給出距離參數(shù)真值足夠近的初始估計(jì)值，若初始值估計(jì)不夠精確，則會(huì)導(dǎo)致迭代不收斂，或收斂到局部最優(yōu)值點(diǎn)上，而非總體最優(yōu)值點(diǎn)，其收斂區(qū)間大都不夠明確。由于均勻采樣條件下正弦波形參數(shù)的初始值估計(jì)可以有很多方法，如FFT方法、峰值檢測(cè)法、平均值法等，且足夠準(zhǔn)確，因此，以往的四參數(shù)正弦波擬合方法主要是針對(duì)均勻采樣條件下的正弦波形序列，在非均勻采樣序列中，初始值參數(shù)估計(jì)將面臨更大困難，甚至無(wú)法進(jìn)行有效估計(jì)，導(dǎo)致它們很多時(shí)候無(wú)法適應(yīng)非均勻采樣正弦波序列的最小二乘擬合要求。

本文后續(xù)工作，將主要針對(duì)上述這些非均勻采樣正弦波序列的最小二乘擬合開(kāi)展研究，以尋找出一種適合非均勻采樣正弦序列波形參數(shù)擬合的最小二乘算法。

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

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

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

式中：E為正弦波形幅度；f為正弦波形頻率；φ為正弦波形初始相位角；Q為正弦波形信號(hào)的直流分量值。

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

2.2 問(wèn)題討論

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

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

圖1 多周期時(shí)歸一化誤差ρ／E與數(shù)字角頻率比w／ω關(guān)系曲線(xiàn)

由于數(shù)字角頻率與頻率存在式（7）所述確定關(guān)系，故擬合過(guò)程中的數(shù)字角頻率比也是擬合過(guò)程中的信號(hào)頻率比，完全可以代表其頻率擬合特征。

其中，ω＝2π／1000；E＝100；Q＝0；φ分別取0°、40°、80°、120°、160°、200°、240°、280°、320°、360°；以w代替ω進(jìn)行三參數(shù)擬合。

從圖1及其局部細(xì)化曲線(xiàn)波形圖2可見(jiàn)，對(duì)于信號(hào)為數(shù)字角頻率ω的正弦曲線(xiàn)，三參數(shù)最小二乘擬合法在ω±ω／10（即［0.9ω，1.1ω］）范圍內(nèi)極值存在且唯一，幅度和相位的變化不影響其ρ／E的變化規(guī)律，可在ω±ω／10范圍內(nèi)通過(guò)對(duì)ω的一維搜索找出該極值點(diǎn)。

圖2 多周期時(shí)歸一化誤差ρ／E與數(shù)字角頻率比w／ω關(guān)系曲線(xiàn)（局部）

在全頻率范圍內(nèi)，幅度E的變化基本上不影響ρ／E隨w／ω的變化時(shí)極值點(diǎn)出現(xiàn)的位置和幅度；相位變化也不影響ρ／E隨w／ω變化時(shí)極值點(diǎn)出現(xiàn)的位置，而只改變非ω頻率處各個(gè)ρ／E極小值的大小，對(duì)各極大值則無(wú)影響。

從變化p值等參數(shù)的多組仿真運(yùn)算不難看出，一般情況下，對(duì)于含有p個(gè)周期的頻率為f的正弦波形序列，按上述三參數(shù)正弦曲線(xiàn)擬合有如下規(guī)律：

（1）在f±f／p（即［f－f／p，f＋f／p］）范圍內(nèi)，ρ／E的極值存在且唯一，極值點(diǎn)就是f頻率點(diǎn)；

（2）幅度和相位的變化不影響ρ／E的變化規(guī)律，但初始相位φ的變化將影響ρ／E的幅度值；

（3）可以通過(guò)f的初始估計(jì)值以及p值估計(jì)和判斷四參數(shù)正弦波曲線(xiàn)擬合算法的收斂區(qū)間［f－f／p，f＋f／p］；

（4）可在f±f／p范圍內(nèi)通過(guò)對(duì)f的一維搜索找出ρ／E的極值點(diǎn)！并且，該極值點(diǎn)處的最小二乘擬合結(jié)果就是四參數(shù)最小二乘正弦波擬合結(jié)果。

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

3.1 原理過(guò)程

假設(shè)待估計(jì)的正弦波頻率目標(biāo)值為f0，待估計(jì)的正弦波采樣序列所含信號(hào)周期個(gè)數(shù)為p；則有，最大頻率差值Δfmax＝f0／p＝1／（tn－t1），在區(qū)間［f0－Δfmax，f0＋Δfmax］內(nèi)的任意頻率f下，殘差平方和ε（f）的極值存在且唯一！這樣，便將四參數(shù)正弦波曲線(xiàn)擬合中，對(duì)幅度、頻率、相位、直流分量4個(gè)參數(shù)的四維非線(xiàn)性搜索，變成了對(duì)頻率分量f造成的ε（f）的一維線(xiàn)性搜索，可保證在區(qū)間［f0－Δfmax，f0＋Δfmax］內(nèi)，用三參數(shù)擬合法實(shí)現(xiàn)的四參數(shù)正弦曲線(xiàn)擬合過(guò)程絕對(duì)收斂。該四參數(shù)擬合過(guò)程如下：

（1）設(shè)定擬合迭代停止條件為he；he為用于判定迭代殘差增量的一個(gè)非常小的小數(shù)值，例如可令he＝10－15或更小。

否則，重復(fù)（4）～（6）的過(guò)程。

3.2 收斂性問(wèn)題

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

在該絕對(duì)收斂區(qū)間之外，當(dāng)p≥2時(shí)，首先，可以對(duì)測(cè)量序列的有效值Erms進(jìn)行預(yù)估計(jì)，

不失一般性，假設(shè)正弦序列的噪聲與信號(hào)有效值幅度之比為N／S?1，即噪聲功率遠(yuǎn)小于信號(hào)功率，則可選取判據(jù)hd取值滿(mǎn)足N／S＜hd＜1。

然后，變化f在（0，＋∞）區(qū)間內(nèi)搜索從ρ（f）／Erms＞hd變化到ρ（f）／Erms＜hd處的頻率點(diǎn)fd，根據(jù)在fd附近ρ（fd）的變化規(guī)律可以判斷：

當(dāng)導(dǎo)數(shù)ρ′（fd）＜0，fd＜f0，可令fL＝fd；

當(dāng)導(dǎo)數(shù)ρ′（fd）＞0，fd＞f0，可令fR＝fd。

這樣，便尋找出落在絕對(duì)收斂區(qū)間內(nèi)且包含f0的迭代區(qū)間［fL，fR］，使用前述方法可獲得絕對(duì)收斂的四參數(shù)正弦擬合結(jié)果。

采取輔助判據(jù)hd這樣一種措施后，可將本文所述四參數(shù)正弦擬合方法的頻率絕對(duì)收斂區(qū)間拓展到（0，＋∞），從而可在任何情況下都獲得收斂結(jié)果。

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

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

圖3 正弦擬合曲線(xiàn)y＾（i）與均勻采樣曲線(xiàn)yi

令正弦波模型參數(shù)不變，使用隨機(jī)采樣間隔進(jìn)行序列采樣，獲得采樣序列并按照本文上述方法進(jìn)行正弦波曲線(xiàn)擬合，得采樣序列曲線(xiàn)及其擬合曲線(xiàn)分別如圖5～圖8所示。

圖4 正弦擬合曲線(xiàn)與均勻采樣曲線(xiàn)差異yi－y＾（i）

表1 正弦波采樣序列波形參數(shù)擬合結(jié)果

其中，圖5為波形長(zhǎng)度約為6個(gè)周期的非均勻采樣序列曲線(xiàn)及其擬合曲線(xiàn)，其采樣間隔為（RND－0.4）×20τ0。τ0＝0.125 ms。其中，RDN為在［0，1］內(nèi)均勻分布的隨機(jī)數(shù)?？梢?jiàn)由于采樣間隔呈非均勻隨機(jī)變化狀態(tài)，按照等間隔均勻采樣模式繪制的曲線(xiàn)圖5，其中的“正弦波形”已經(jīng)變化非常大。若仍然按照均勻采樣間隔處理，將無(wú)法獲得有效擬合結(jié)果，甚至?xí)M合失敗。圖6為擬合曲線(xiàn)與采樣序列之差異曲線(xiàn)。可見(jiàn)兩者擬合得非常好。相應(yīng)擬合參數(shù)見(jiàn)表1。由表1可見(jiàn)，非均勻采樣條件下正弦曲線(xiàn)擬合誤差要略大于均勻采樣條件下，但差異不太大，造成差異的原因有待于進(jìn)一步研究。

圖7為圖5的測(cè)量曲線(xiàn)不經(jīng)過(guò)時(shí)間從小到大的排序直接進(jìn)行四參數(shù)擬合獲得的擬合結(jié)果與測(cè)量序列曲線(xiàn)，圖8為此情況下兩者之間的差異曲線(xiàn)，擬合參數(shù)如表1所示。由表1數(shù)據(jù)及圖7和圖8可見(jiàn)，未經(jīng)時(shí)間順序排序的采集波形擬合參數(shù)與經(jīng)過(guò)排序后再執(zhí)行擬合的結(jié)果沒(méi)有本質(zhì)差異。

圖5 正弦擬合曲線(xiàn)y＾（i）與非均勻采樣曲線(xiàn)yi（先進(jìn)行時(shí)間排序）

圖6 正弦擬合曲線(xiàn)與非均勻采樣曲線(xiàn)差異yi－（先進(jìn)行時(shí)間排序）

圖7 擬合曲線(xiàn)）與非均勻采樣曲線(xiàn)yi（未進(jìn)行時(shí)間排序）

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

用FLUKE5700多功能校準(zhǔn)器作激勵(lì)源，給出幅度8.0000V的正弦電壓信號(hào)，頻率10.000Hz，以北京阿爾泰公司ART2001型數(shù)據(jù)采集系統(tǒng)執(zhí)行采集，其A／D位數(shù)12 Bit，工作量程±10 V，采樣速率5000 Sa／s，采集數(shù)據(jù)個(gè)數(shù)2000點(diǎn)。

使用上述四參數(shù)正弦波擬合方法獲得其正弦波形幅度為8.016459 V，頻率為10.00005 Hz，相位為－17.9767°，直流分量為0.483 mV。殘差有效值ρ＝2.21 mV。其測(cè)量序列波形及擬合曲線(xiàn)如圖9所示，測(cè)量序列波形與擬合曲線(xiàn)之差如圖10所示。

圖8 擬合曲線(xiàn)與非均勻采樣曲線(xiàn)差異yi－（未進(jìn)行時(shí)間排序）

圖9 擬合曲線(xiàn)）與均勻采樣曲線(xiàn)yi

圖10 擬合曲線(xiàn)與均勻采樣曲線(xiàn)差異yi－

將上述測(cè)量曲線(xiàn)波形截取三段后，變成非均勻采樣序列，使用上述四參數(shù)正弦波擬合方法獲得其正弦波形幅度為8.016445V，頻率為10.00003Hz，相位為－17.9755°，直流分量為0.443mV。殘差有效值ρ＝2.22 mV。其測(cè)量序列波形及擬合曲線(xiàn)如圖11所示，測(cè)量序列波形與擬合曲線(xiàn)之差如圖12所示。

6 討 論

圖11 擬合曲線(xiàn)與非均勻采樣曲線(xiàn)yi

圖12 擬合曲線(xiàn)與非均勻采樣曲線(xiàn)差異yi－

從上述仿真和實(shí)際實(shí)驗(yàn)結(jié)果可見(jiàn)，本文所述正弦波擬合測(cè)量方法可以用于非均勻采樣正弦波形參數(shù)的測(cè)量估計(jì)，并可以給出幅度、頻率、相位、直流分量等基本信息參量。其最大特點(diǎn)是針對(duì)同時(shí)具有采樣幅度信息和采樣時(shí)間信息的正弦波序列，不論其采樣間隔是否均勻，均可以獲得有效收斂結(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í)際實(shí)驗(yàn)中，兩者的差異不明顯。

7 結(jié) 論

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

［1］ Kuffel J，M cComb T R，Malewski R.Comparative Evaluation of Computer Methods for Calculating the Best-Fit Sinusoid to the Digital Record of a High-Purity Sine Wave［J］.IEEETransactionsonInstrumentationand Measurement，1987，36（2）：418－422.

［2］ McComb T R，Kuffel J，Le Roux B C.A Comparative Evaluation of Some Practical Algorithms Used in the Effective Bits Test of Waveform Recorders［J］.IEEE TransactionsonInstrumentationandMeasurement，1989，38（1）：37－42.

［3］ Jeng Y C.High precision sinusoidal frequency estimator based on weighted least square method［J］.IEEE TransactionsonInstrumentationandMeasurement，1987，36（1）：124－127.

［4］ Deyst J P，Souders T M，Solomon O.Bounds on Least-Squares Four-Parameter Sine-Fit Errors Due to Harmonic Distortion and Noise［J］.IEEETransactionson InstrumentationandMeasurement，1995，44（3）：637－642.

［5］ Zhang JQ，Zhao Xinmin，Hu Xiao，etal.Sine wave fit algorithm based on total least-squares method with application to ADC effective bitsmeasurement［J］.IEEE transactionsonInstrumentationandMeasurement，1997，46（4）：1026－1030.

［6］ Handel P.Properties of the IEEE-STD-1057 Four-Parameter Sine Wave Fit Algorithm［J］.IEEE TransactionsonInstrumentationandMeasurement，2000，49（6）：1 189－1193.

［7］ 田社平，王堅(jiān)，顏德田，等.基于遺傳算法的正弦波信號(hào)參數(shù)提取方法［J］.計(jì)量技術(shù)，2005，（5）：3－5.

［8］ 梁志國(guó)，朱濟(jì)杰，孟曉風(fēng).四參數(shù)正弦曲線(xiàn)擬合的一種收斂算法［J］，儀器儀表學(xué)報(bào)，2006，27（11）：1513－1519.

［9］ 梁志國(guó)，孟曉風(fēng).殘周期正弦波形的四參數(shù)擬合［J］.計(jì)量學(xué)報(bào)，2009，30（3）：245－249.

［10］ Bilau T Z，Megyeri T，Sárhegyi A，etal.Fourparameter fitting of sinewave testing result：Iteration and convergence［J］.ComputerStandards＆Interfaces，2004，26（1）：51－56.

［11］ Palfi V，Kollar I.Acceleration of the ADC Test With Sine-Wave Fit［J］.IEEETransactionsonInstrumentation andMeasurement，2013，62（5）：880－888.

［12］ Chen K F.Estimating Parameters of a Sine Wave by Separable Nonlinear Least Squares Fitting［J］.IEEE TransactionsonInstrumentationandMeasurement，2010，59（12）：3214－3217.

［13］ Sárhegyi A，Kollár I.Robust Sine Wave Fitting in ADC Testing［C］／／Proceedings of the IEEE Instrumentation and Measurement Technology Conference，2006.Sorrento，Italy，2006，24－27.

［14］ IEEE Std 1057-1994.IEEE standard for digitizing waveform recorders［S］.1994.

The Sinewave Fit Algorithm Based on Total Least-square Method w ith Non-uniform Samp ling

LIANG Zhi-guo1， ZHU Zhen-yu2
（1.Changcheng Institute of Metrology and Measurement，Beijing 100095，China；
2.College of Precision Instrument and Opto-electronics Engineering，Tianjin University，Tianjin 300072，China）

Aiming at the non-uniform sampling series，a four-parameter sine wave curve-fitmethod is introduced.It is based on the three-parameter sine wave curve-fit method with frequency known.It can attain the four-parameter of sinusoidal curve-fit results with both the uniform sampling series and the non-uniform sampling series.The speciality of the arithmetic is that it turns the optimization of four parameters（amplitude，frequency，phase and offset）into the optimization of one parameter（only frequency），and without any original parameter estimation.Themethod has excellent convergence，good robustness，and clearly convergence interval.The validity and feasibility are proved by both simulation and experiments，thismethod can be applied to the four-parameter sinewave curve-fitwith non-uniform sampling series.

Metrology；Non-uniform sampling；Sinusoidal；Curve-fit；Parameter estimation；Calibration

TB973

A

1000-1158（2014）05-0494-06

10．3969／j．issn．1000-1158．2014．05．18

2013-05-17；

2013-09-23

梁志國(guó)（1962－），男，黑龍江巴彥人，博士，北京長(zhǎng)城計(jì)量測(cè)試技術(shù)研究所研究員，主要研究方向?yàn)閿?shù)字化測(cè)量與校準(zhǔn)、模式識(shí)別、動(dòng)態(tài)校準(zhǔn)、精確測(cè)量。Lzg304＠sina.com