梁志國， 朱振宇

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

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

梁志國1， 朱振宇2

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

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

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

1 引 言

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

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

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

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

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

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

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

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

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

2.2 問題討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 原理過程

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

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

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

3.2 收斂性問題

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

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

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

然后，變化f在（0，＋∞）區(qū)間內(nèi)搜索從ρ（f）／Erms＞hd變化到ρ（f）／Erms＜hd處的頻率點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。

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

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

4 仿真實驗驗證

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

圖3 正弦擬合曲線y＾（i）與均勻采樣曲線yi

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

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

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

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

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

圖5 正弦擬合曲線y＾（i）與非均勻采樣曲線yi（先進行時間排序）

圖6 正弦擬合曲線與非均勻采樣曲線差異yi－（先進行時間排序）

圖7 擬合曲線）與非均勻采樣曲線yi（未進行時間排序）

5 實驗驗證

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

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

圖8 擬合曲線與非均勻采樣曲線差異yi－（未進行時間排序）

圖9 擬合曲線）與均勻采樣曲線yi

圖10 擬合曲線與均勻采樣曲線差異yi－

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

6 討 論

圖11 擬合曲線與非均勻采樣曲線yi

圖12 擬合曲線與非均勻采樣曲線差異yi－

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

仿真實驗表明，與均勻采樣條件相比，非均勻采樣序列正弦擬合誤差要略大一些。實際實驗中，兩者的差異不明顯。

7 結(jié) 論

綜上所述，本文提出了非均勻采樣條件下正弦波測量序列四參數(shù)擬合的一種方法，給出了詳細過程和收斂區(qū)間，并分別在隨機間隔采樣和非隨機間隔采樣兩種條件下給出了比較結(jié)果。結(jié)果表明，本文所述方法具有良好收斂性和魯棒性，僅要求已知采樣序列和每個采樣點的時刻，沒有任何其它先決條件要求，是一種表現(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］ 田社平，王堅，顏德田，等.基于遺傳算法的正弦波信號參數(shù)提取方法［J］.計量技術(shù)，2005，（5）：3－5.

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

［9］ 梁志國，孟曉風(fēng).殘周期正弦波形的四參數(shù)擬合［J］.計量學(xué)報，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

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