基于比時法的晶振頻率測量建模與分析

摘 要:受器件老化、隨機噪聲等因素影響，晶振頻率變化較復雜。以GPS秒脈沖作為測量標準，構建了晶振頻率隨時間變化的測量系統，通過對測量數據進行一元回歸統計處理，分離出了晶振實際頻率與其標稱頻率的相對偏差及晶振的各種隨機誤差，并分析了這兩種誤差對晶振準確度及穩定度的影響。該方法可為頻率源誤差測量分析提供借鑒作用。

關鍵詞:GPS;晶振;頻率;誤差;回歸分析

中圖分類號:TM930 文獻標識碼:A

文章編號:1004-373X(2010)03-119-03

Modeling and Analysis of Crystal-oscillator Frequency Based on Measuring Time

LI Erpeng，WEN Kaizhang，WANG Weiwei

(Northwest Institute of Nuclear Technology，Xi′an，710024，China)

Abstract:The characteristics of oscillator frequency error are very complex owing to the aging and stochastic noise.The system measure model of crystal oscillator frequency error is designed with GPS as frequency standard.A method of regression analysis is used for disposing the oscillator frequency error data and separating the linearity excursion error from the stochastic error.Both two errors are analyzed in contribution of the crystal oscillator stabilization and precision.This mean is provided with reference in frequency criterion error of measure and analysis.

Keywords:GPS;crystal-oscillator;frequency;error;regression analysis

0 引 言

石英晶體振蕩器受制造工藝、器件老化以及外部溫度等因素影響，其實際頻率值與標稱頻率值存在偏差。此外，晶體振蕩器內部存在著各種噪聲也會使頻率值產生隨機起伏，從而導致晶振頻率的準確度和穩定度降低。當前，國內外在研究晶體振蕩器的老化和隨機噪聲、分析晶振頻率誤差特性方面提出了很多方法，比較常見的有時間對數線性模型法、自適應濾波法及非線性時變預測法。前兩種方法的缺點是參數較多，選擇合適的參數較難[1，2]，非線性時變預測法的模型描述能力強，但由于沒有函數參數的顯式表達式，不能將模型求解歸結為參數求解問題，一般通過學習來逼近該函數，主要用于頻率變化的預測[3]。

在上述方法的基礎上，本文利用CPLD設計了一種電路，該電路采用比時法[4，5]來測量晶振頻率變化，并根據其頻率隨時間變化擬合曲線的特點，用線性回歸法分析其頻率誤差特性。該方法模型簡單，參數易于估計，可通過簡潔的補償方法消除晶振相對頻偏，具有實際的應用價值。

1 晶振頻率測量系統組成

測量系統由GPS接收機、晶振、時差測量模塊、時鐘產生模塊、計算機數據采集處理組成。組成原理如圖1所示。

圖1 測量系統組成

GPS接收機每秒輸出1路TTL電平的標準秒脈沖(1 PPS)，晶振是時差測量和時鐘產生的頻率源。時鐘產生電路產生本地秒脈沖。時差測量電路測量GPS秒脈沖與本地秒脈沖的相位差值。

采用比時法測量晶振頻率的系統工作原理如下:首先由晶振分頻得到本地的1 Hz頻率源，將GPS秒脈沖與本地晶振秒脈沖送入時差測量模塊進行相位比較，得到兩者的相位差信號，設計時間間隔計數器[6]對此相差閘門信號計數，計數值即為晶振頻率相對于標準頻率的偏差，反映了晶振頻率的誤差特征。時間間隔計數器每秒測量一次兩者的偏差值，需要測量的頻差數據量由計算機設定，測量結果傳送到計算機進行數據統計，并對晶振頻率誤差特性進行分析。本地時鐘產生、時差測量及數據采集電路模塊等都由CPLD設計實現[7，8]。

2 模型建立

實驗以GPS秒脈沖作為標準參考信號，采用比時法對晶振的輸出頻率進行測量。按照建立的測量系統，實際對某一晶振采集到30個計數值，如表1所示。

表1 實驗數據

x/s12345678910

y/個79156235313392470548626704783

x/s11121314151617181920

y/個8619401 0181 0971 1751 2541 3321 4101 4881 566

x/s21222324252627282930

y/個1 6451 7221 8011 8791 9562 0342 1132 1932 2702 349

表1中x表示測量的時間(單位:s);y表示晶振脈沖計數個數。假設晶振在某秒計數值為M，它的計數周期為T，則MT為晶振秒脈沖與GPS秒脈沖的時間差值。例如在第6 s時測得計數值為470，則表示在第6 s時晶振與GPS秒脈沖的時間差為470T。為了研究時間x與計數個數y之間的關系，用ORIGIN工具軟件對數據進行擬合處理，得到的x，y關系曲線如圖2所示。

圖2 數據擬合曲線

從散點圖可以看出，測量計數值和測量時間大致呈線性關系。據此假設這兩個變量之間的內在關系是一條直線，這些點與直線的偏離是由于測量過程中其他一些隨機因素的影響而引起的[9]，這樣可以假設這組測量數據有如下結構形式:

yt=β0+βxt+εt，t=1，2，…，N

(1)

式中:ε1，ε2，…，εN分別表示其他隨機因素對變量y1，y2，…，yN影響的總和，一般假設它們是一組相互獨立，并服從同一正態分布N(0，δ)的隨機變量。變量x在實驗中為自然數，表示具體的秒脈沖數值。這樣，變量y表示實際所測得的晶振與標準頻率的計數差值，它是服從N(β0+βxt，δ)的隨機變量。

用最小二乘法[10]來估計參數β0，β。設b0，b分別是參數β0，β的最小二乘估計，于是得到一元線性回歸的回歸方程:

=b0+b

(2)

式中:b0，b是回歸方程的回歸系數，分別表示晶振相對于標準頻率的初始誤差和累積誤差。應用最小二乘法可求得回歸系數b，b0為:

b=N∑Nt=1xtyt-(∑Nt=1xt)(∑Nt=1yt)N∑Nt=1x2t-(∑Nt=1xt)2=lxylxx

(3)

b0=(∑Nt=1x2t)(∑Nt=1yt)-(∑Nt=1xt)(∑Nt=1xtyt)N∑Nt=1x2t-(∑Nt=1xt)2=-b

(4)

3 數據分析與處理

3.1 回歸系數估計

為了定量分析數據，從而確定晶振頻率誤差的組成，首先對上述測得的數據進行歸一化處理。實際測量中得到的是晶振脈沖的計數個數，設測量系統所用晶振頻率為10 MHz，可將計數數據轉化為晶振相對于標準時間每秒的時間之差。例如在x=30 s時，y=2 349，表示在第30 s時，晶振頻率相對于標準頻率的計數值為2 349，若晶振頻率f=10 MHz，則可得到在第30 s時晶振相對于GPS時間的誤差為t=y/f=2349 μs。

用Matlab對歸一化數據進行處理，依照最小二乘原理，得到計數時間x與時間差值y的均值，以及x的自相關、x和y的互相關、y的自相關及回歸方程如表2所示。

表2 數據分析表

=∑30t=1xtN=15.5 s=∑30t=1ytN=1.213 6 μsN=30

lxx=∑Nt=1(xt-)2

=2 247.5 s

lyy=∑Nt=1(yt-)2

=1.377 2×105 μs2

lxy=∑Nt=1(xt-)#8226;

(yt-)

=1.759 3×104 μs#8226;s

b=lxylxx=7.828 μs/s

b0=-b=0.030 μs

=b0+bx=0.030+7.828x

從表2得到回歸方程為:

=b0+bx=0.030+7.828x

(5)

從回歸方程可以看出，給定一個時間x值，就可得到相對應的時間差值。例如，在x=15 s時，=117.72 μs，表示在第15 s時晶振相對于標準時間的誤差為t=117.72 μs。從回歸方程可預測此晶振在24 h相對于標準時間的時間誤差值為:

7.828×3 600×24=676 339.2 μs676 ms

以上求得了回歸方程，但是該方程是否基本上符合y與x之間的客觀規律，是否符合晶振頻率誤差變化的實際特點，還需要對回歸方程做進一步的分析。在回歸分析法中，通常采用方差分析法[10]對回歸方程的顯著性進行檢驗，其實質是將N個測量值的影響從數量上區分開，然后用F檢驗法對所求回歸方程進行顯著性檢驗。

3.2 晶振頻率誤差數據方差分析及顯著性檢驗

測量值y1，y2，…，yN之間的差異(稱為變差)是由兩方面的原因引起的。一是自變量取值的不同，二是其他因素(包括試驗誤差、隨機誤差等)的影響。為了對回歸方程進行檢驗，把兩者所引起的變差從y的總變差中分解出來。根據上述數據可得:

U=∑Nt=1(t-)2=blxy=1.377 18×105

Q=∑Nt=1(yt-t)2=lyy-blxy=1.996

式中:

U稱為回歸平方和，它反映了在y的總變差中由于x和y的線性關系而引起的y變化的部分;Q稱為殘余平方和，即所有測量點距回歸直線的殘余誤差平方和。

若總的平方和由N項組成，其自由度就為N-1，總的離差平方和的自由度可分為回歸平方和的自由度vU和殘余平方和的自由度vQ之和，即:

vS=vU+vQ

(6)

在一元線性回歸中，vU=1，vS=N-1，則Q的自由度vQ=N-2。由回歸平方和與殘余平方和的意義可知，一個回歸方程是否顯著，也就是y與x的線性關系是否密切，取決于U及Q的大小，U愈大Q愈小，說明y與x的線性關系愈密切。通常可以采用F檢驗法來對方程進行顯著性檢驗。對于一元線性回歸，將U及Q的值代入上式得到統計量F:

F=U/1Q/(N-2)=1.931 528×106



由F分布表可以查出，F≥F0.01(1，28)=7.64。可認為回歸是高度顯著的。

殘余平方和Q除以它的自由度vQ所得商:

σ2=Q/(N-2)=0.071 3 μs2

稱為殘余方差，它可以看作排除了x對y的線性影響后，衡量y隨機波動大小的一個估計量。殘余方差的平方根:

σ=Q/(N-2)=0.267 μs

稱為殘余標準差，它可用來衡量所有隨機因素對y一次性測量平均變差的大小，σ愈小，則回歸直線的精度愈高。

把平方和及自由度進行分解的方差分析數據結果歸納在一個表格中，如表3所示。

表3 線性模型數據的方差分析表

來源平方和 /μs2自由度方差 /μs2F顯著性

回歸

殘余1.377 18×105

1.996128…

0.071 31.931 528×106α=0.01

…

總計137 719.99629………

從表3可以看出，在30 s時間內，晶振實際頻率與其標稱頻率的相對偏差引起的誤差平方和為1.377 18×105 μs2，其他各種隨機因素引起的誤差的平方和為1.996 μs2。晶振相對頻偏引起的誤差遠遠大于其他因素引起的誤差，表明了晶振頻率誤差貢獻主要來源于實際頻率與標稱頻率之間的頻偏。通過頻率修正即可得到一種高穩定度的頻標源。

4 結 語

以上分析結果表明，該晶振頻率誤差主要來源于實際頻率與標稱頻率的相對偏差，隨機誤差對晶振整體誤差的貢獻很小。通過對晶振標稱頻率值進行在線補償可以消除其相對頻偏，從而獲得一種具有較高穩定度的頻率源，可為需要時間顯示的場所提供高精度的時間服務。

參考文獻

[1]Su W，FILLER R L.Application of Kalman Filtering Techniques to the Precision Clock with Non-constant Aging [A].Proc.IEEE Frequency Control Symposium[C].1992:231-237.

[2]SU W.A Novel Method for Aging Estimation of Crystal Oscillators[A].Proc.IEEE Frequency Symposium[C].1996:890-896.

[3]柳麗，陳之純.石英晶振的建模與模型求解[J].上海航天，2004(3):1-6.

[4]曾祥君，尹項根，林干，等.晶振信號同步GPS信號產生高精度時鐘的方法及實現[J].電力系統自動化，2003，27(8):49-51.

[5]葛先軍.時間頻率計量技術[D].煙臺:海軍航空工程學院，2005.

[6]魏西峰.全同步數字頻率測量方法的研究[J].現代電子技術，2005，28(12):101-102.

[7]周立功.單片機與CPLD綜合應用技術[M].北京:北京航空航天大學出版社，2003.

[8]馬忠梅.單片機的C語言應用程序設計[M].北京:北京航空航天大學出版社，2003.

[9]寇艷紅，張其善.GPS接收機中晶振誤差模擬方法[J].電子與信息學報，2004，26(8):1 320-1 324.

[10]費業泰.誤差理論與數據處理[M].4版.北京:機械工業出版社，2003.