最小平方頻譜分析法

最小平方頻譜分析法（英語：Least-squares spectral analysis）是一種利用最小平方法尋找適配於資料點之最佳正弦曲線，以估算頻譜的方法。其數學原理與科學界中最常用的傅立葉分析相似^[1][2]。一般而言，傅立葉分析會將間隔較長之訊號的長周期雜訊放大，而最小平方頻譜分析法則解決了這個問題^[3]。

最小平方頻譜分析法也稱為凡尼切克法（Vaníček method）^[4]、隆布法（Lomb method）^[3][5]或隆布—史卡構法（Lomb–Scargle method）^[2][6][7]，分別取名自對其有所貢獻的佩特·凡尼切克（捷克語：Petr Vaníček）、尼可拉斯·隆布（Nicholas R. Lomb）^[8]以及傑佛瑞·史卡構（Jeffrey D. Scargle）^[9]。此外，麥可·科恩伯格（Michael Korenberg）、史考特·陳（Scott Chen）以及大衛·多諾霍（英語：David Donoho）等人也曾開發出與之關係密切的其它方法。

歷史背景

長久以來，科學界對於傅立葉分析、週期圖法與正弦曲線之最小平方擬合間的密切關聯早有所知^[10]。然而，大多數以上述理論為基礎開發的方法僅適用於取樣間距相等的訊號。1963年，荷蘭數學和計算機科學研究學會的弗里克·巴寧（Freek J. M. Barning）提出類似的方法處理了取樣間距不同的訊號^[11]：先以週期圖法分析（現又稱為隆布法），再從所得的週期圖中選取特定頻率的正弦曲線並以最小平方法擬合之；這兩段過程間則是透過匹配追蹤連接，並以反擬合後處理（post-backfitting）以避免過適^[12]、或使用正交匹配追蹤法^[13]。

新不倫瑞克大學的加拿大大地測量學家佩特·凡尼切克（捷克語：Petr Vaníček）於1969年也提出了匹配追蹤法，他稱之為「連續頻譜分析」，並稱其所得的結果為「最小平方週期圖」，且同時適用於取樣間距相同或不同的訊號^[14]。他更於1971年將該方法一般化，使其能夠分析出訊號中單一平均值以外的各種系統成分，包含可預測但大小未知的線性、二次方或指數趨勢^[15]。

雪梨大學的尼可拉斯·隆布於1976年簡化了凡尼切克法，並指出該方法與週期圖法有相當密切的關連^[8]。後來，傑佛瑞·史卡構也定義了非等距取樣訊號的週期圖，並對其做了分析^[9]；他並發現，只要將該週期圖稍作修改，便能夠得到與隆布的最小平方擬合法完全相同的結果。

史卡構曾表示他的論文「並未提出一個新的測量方法，而是對現有最常用的方法——週期圖法——在取樣間距不相等時的可靠性以及效率進行研究」。他在論文中同時引用了最小平方弦波擬合以及週期圖法分析，並指出該兩種方法（在他所提出的修改之下）是完全相等的^[9]。

媒體對相關研究的發展則有如下的摘要：

「隆布提出了一項適用於非等距取樣訊號，且完全不同的頻譜分析法。它不僅減少了舊有方法的困難，更具備一些非常理想的性質。此方法乃是基於巴寧和凡尼切克早前的研究，並後為史卡構做了更詳盡的闡述^[3]。」

1989年，皇后大學的麥可·科恩伯格開發了所謂的「快速正交搜尋法」（fast orthogonal search），以更迅速地找到一組頻譜成分分析的近似最佳解^[16]；該方法與後來出現的正交匹配追蹤法有許多相似之處。1994年，史丹佛大學的史考特·陳與大衛·多諾霍提出了「基底追蹤法」（basis pursuit）；該方法藉由將目標轉化為線性規劃問題，並使各項係數的L₁-距離最小化，以更快速地得到解答^[17]。

方法

凡尼切克法

凡尼切克法乃是利用標準的線性回歸或最小平方法，將一個離散訊號近似為許多不同頻率的弦波之加權總和^[18]。這些弦波的頻率則是以類似巴寧所提出的方法來決定，並另外藉由減少最小平方擬合後的殘差以使每次決定的新頻率得到最佳化（這與匹配追蹤法加上反擬合前處理所得到的結果相同^[12]）。並且，決定後的弦波數量必須不超過整個訊號的取樣點數量（同一頻率的正弦和餘弦波需視為不同的弦波）。

假設${\displaystyle \phi }$為一個由許多不同頻率的基底弦波加權總和所得的訊號向量，其基底函數以矩陣${\displaystyle {\textbf {A}}}$表示，而每個基底所對應的權重則以${\displaystyle x}$向量表示，也就是：

${\displaystyle \phi \approx {\textbf {A}}x}$

會希望權重向量${\displaystyle x}$使得上式在近似${\displaystyle \phi }$時的平方誤差和為最小。利用標準的線性回歸法即可求得${\displaystyle x}$的閉合形式解為^[19]：

${\displaystyle x=({\textbf {A}}^{\mathrm {T} }{\textbf {A}})^{-1}{\textbf {A}}^{\mathrm {T} }\phi }$

其中，矩陣${\displaystyle {\textbf {A}}}$可由任意彼此獨立（但不需正交）的基底函數所組成；對頻譜分析而言，這些基底函數一般為某特定頻率範圍內均勻分布的正弦或餘弦波。如果在某過窄的頻寬內選擇了過多的頻率，則這些所選的函數可能不夠獨立、矩陣${\displaystyle {\textbf {A}}}$本身較為病態，而使得最後求得的頻譜不具意義^[19]。

當矩陣${\displaystyle {\textbf {A}}}$的所有基底函數皆彼此正交（也就是彼此不相關，且矩陣中的每一行兩兩內積皆為零）時，${\displaystyle {\textbf {A}}^{\mathrm {T} }{\textbf {A}}}$為一對角矩陣；此外，當${\displaystyle {\textbf {A}}}$之每一行的能量（即向量之元素平方和）皆相等時，${\displaystyle {\textbf {A}}^{\mathrm {T} }{\textbf {A}}}$便為一單位矩陣與某常數相乘的積，而其反矩陣的計算便相當容易。假設取樣間距為${\displaystyle t}$，那麼後者的情況即是當取樣間距相同，且基底函數選擇的是頻率範圍為${\displaystyle (0,1/2t)}$內成對的正弦和餘弦波（每個基底的頻率為${\displaystyle 1/Nt}$，並忽略區間中取樣值皆為零的最小和最大頻率），而這情況同時也是離散傅立葉變換的特例之一^[19]。

${\displaystyle x={\textbf {A}}^{\mathrm {T} }\phi }$     （離散傅立葉變換中${\displaystyle N}$個等距取樣點和頻率乘上某比例常數的情況）

隆布所提出上述的簡化式可用於大多數的情況，因為一般而言不同弦波彼此之間的相關係數較小（至少當它們的間距夠大時），惟同一頻率下成對的正弦和餘弦波彼此間並非獨立而不適用。本質上，這和傳統週期圖法的公式無異，但可用於非等距取樣的訊號。其中，${\displaystyle x}$向量可為訊號的潛在頻譜做有效的估計，但由於基底之間的相關性被忽略，${\displaystyle {\textbf {A}}x}$便不再能夠良好地近似原始的訊號，這個方法也因此不再是真正的「最小平方法」，惟學界仍繼續如是稱之。

隆布—史卡構法

在標準的週期圖法公式中，一般是直接計算訊號與正弦或餘弦波之間的內積。史卡構則對其稍作修改：首先，先找到一時間延遲項${\displaystyle \tau }$使得成對的正弦與餘弦波在取樣時間為${\displaystyle t_{j}}$時將會彼此正交；同時，他也對可能具有不同能量的兩個基底之間做了標準化，以取得在某頻率下更加準確的能量估計值^[3][9]。這些修改使得史卡構的週期圖法和隆布的最小平方法即完全相同。時間延遲項${\displaystyle \tau }$的定義如下：

${\displaystyle \tan {2\omega \tau }={\frac {\sum _{j}\sin 2\omega t_{j}}{\sum _{j}\cos 2\omega t_{j}}}}$

頻率${\displaystyle \omega }$的週期圖則可估算如下：

${\displaystyle P_{x}(\omega )={\frac {1}{2}}\left({\frac {\left[\sum _{j}X_{j}\cos \omega (t_{j}-\tau )\right]^{2}}{\sum _{j}\cos ^{2}\omega (t_{j}-\tau )}}+{\frac {\left[\sum _{j}X_{j}\sin \omega (t_{j}-\tau )\right]^{2}}{\sum _{j}\sin ^{2}\omega (t_{j}-\tau )}}\right)}$

史卡構並表示，上式所得的週期圖與等距取樣情況下的結果具有相同的統計分布^[9]。

在任一單獨的頻率${\displaystyle \omega }$下，這個方法求得的頻譜能量與最小平方法使用該頻率之弦波擬合後所得的能量相同。擬合後函數的形式如下：

${\displaystyle \phi (t)=A\sin \omega t+B\cos \omega t}$^[20]

一般化

標準的隆布—史卡構週期圖法適用於平均值為零的模型（即擬合用弦波的集合），因此一般而言在計算週期圖之前會先減去訊號的平均值並假設其為零。然而，當模型的平均值為非零時，這項假設並不準確。一般化的隆布—史卡構週期圖法便將此假設移除，並在計算頻譜時一併求得其平均值。如此一來，擬合後函數的形式如下：

${\displaystyle \phi (t)=A\sin \omega t+B\cos \omega t+C}$^[21]

此一般化的隆布—史卡構週期圖法也被稱為「浮動平均週期圖法」^[22]。

快速正交搜尋法

科恩伯格於1989年提出了如下的方法：先從一過完備（英語：Overcompleteness）的函數集合中選出一稀疏子集（以頻譜分析而言即指弦波）用於擬合，稱為快速正交搜尋法（fast orthogonal search）。數學上，快速正交搜尋是在縮小均方誤差（MSER）時使用了一種略為修改過的楚列斯基分解法，並以稀疏矩陣之求逆方法實作^[16][23]。快速正交搜尋如同其它的最小平方頻譜分析法，也能夠避免離散傅立葉變換最主要的缺點，並且可以求得高準確度的潛在週期性，同時對於處理非等距取樣訊號也相當有效。快速正交搜尋法同時也被應用於非線性系統鑑別（nonlinear system identification）等其它的問題。

基底追蹤法

陳與多諾霍提出的基底追蹤（英語：Basis pursuit）法（basis pursuit）也是從一過完備的函數集合中選出一稀疏子集（其中可包含弦波或其它函數）用於擬合，惟此法定義能夠使各項係數的L₁-距離得到最小化者即為最佳解。如此一來，頻譜分析便可藉此轉化為已知的線性規劃問題，並可利用現有的快速演算法得到解答^[17]。

卡方法

大衛·帕爾默（David Palmer）提出的卡方法則能夠找到任意數量之諧波的最佳擬合函數，提高了尋找非弦波諧和函數的自由度^[24]。這是一項基於快速傅立葉變換、對取樣間距任意且標準誤差不均的訊號進行加權最小平方（英語：Weighted least squares）分析的快速演算法。實作此方法的原始碼也已被公開^[25]。由於離散訊號往往不具有相等的取樣間距，此方法先將訊號「柵格化」，即稀疏地在數個取樣點間填入一串時間序列。所有時間軸上互相交會的柵格點之統計權重都將被設為零，亦即取樣點之間的誤差線（英語：Error bar）為無窮大。

應用

最小平方頻譜分析法最主要的用處是對未完整記錄的訊號進行頻譜分析，且不須對其竄改或摻入不存在的數據。

利用最小平方頻譜分析法所得的頻譜中之各強度量值即代表了某個頻率或週期對於整個時間序列的變異數之貢獻^[14]。一般而言，經由上述方法所定義的頻譜強度可用於推算輸出結果的顯著水準狀態^[26]。此外，凡尼切克法求得的頻譜強度也可由分貝值表示^[27]。值得注意的是，凡尼切克頻譜在統計上具有Β分布^[28]。

凡尼切克的最小平方頻譜之反運算可由下述方法實現：若將正運算的過程表示為一矩陣的乘法，那麼（當它為非奇異矩陣時）求取它的反矩陣或偽反矩陣便可進行反運算；當所選的基底弦波在取樣點時皆彼此獨立且其總數和取樣點數量相等時，反運算所得的結果便會與原始的訊號相同^[19]。相對地，週期圖法並不具有已知的反運算方法。

實作

最小平方頻譜分析法可用一頁以內的MATLAB程式碼實作完成^[29]。簡言之：

「欲求得最小平方頻譜，我們必須計算${\displaystyle m}$個頻譜值……可以藉由${\displaystyle m}$次的最小平方法近似，每次求得一種不同頻率（的頻譜能量）^[18]。」

也就是說，對一所欲擬合之頻率集合中的每個頻率${\displaystyle \omega _{i}}$，將其對應的正弦（${\displaystyle \sin {\omega _{i}t}}$）和餘弦函數（${\displaystyle \cos {\omega _{i}t}}$）於訊號的取樣時間點${\displaystyle t_{j}}$取值，計算訊號與弦波之取樣值的內積並給予適當的標準化；若是使用隆布—史卡構的週期圖法，則在取內積之前先計算每個頻率${\displaystyle \omega _{i}}$的時間延遲項${\displaystyle \tau _{i}}$使得該頻率下的正弦與餘弦成分能夠彼此正交^[19]；最後，從振幅的內積值便可求得該兩項弦波成分的頻譜能量。若將這項程序套用於取樣間距相等、且擬合弦波頻率選為取樣時間全長倒數之整數倍的訊號時，即相當於離散傅立葉變換。

上述即凡尼切克的原始方法，它在處理每個弦波成分時都將其視為彼此獨立，即使在取樣點時這些弦波有可能並不彼此正交。相對地，透過矩陣方程式求解、並將原始訊號的總變異數分割至指定的弦波頻率之間的方法，也能夠在一次考慮所有弦波成分（即不獨立處理）的情況下求得頻譜值^[19]。這種矩陣求解的最小平方法可直接透過MATLAB內建的反斜線運算子得到^[30]。

麥克·克雷默（Mike Craymer）解釋道，相對於獨立處理各成分的作法（以及隆布的週期圖法），一次考慮所有成分的矩陣求解法無法擬合比取樣點數量還多的基底（即正弦與餘弦波），並指出：

「……如果所選的頻率使得某些傅立葉成分（即三角函數）彼此間接近線性獨立，則對於結果可能會有嚴重的影響，進而產生一個病態或近乎奇異的矩陣${\displaystyle N}$。若想避免產生這樣的病態條件，則必須選取一組不同的頻率來估計（例如等差頻率）、或是直接忽略${\displaystyle N}$的相關性（即非對角線上的值）並分別估算各獨立頻率的最小平方反轉換……^[19]」

相對地，隆布提出的方法便如同標準的週期圖法，可使用任意多數或任意高密度的頻率成分做為擬合用的基底；也就是說，這個方法可在頻域下以任意的倍率進行過取樣^[3]。

在傅立葉分析（例如傅立葉變換或離散傅立葉變換）中，用於擬合訊號的弦波皆為彼此正交，所以前述的獨立投影內積法和一次計算的矩陣求解最小平方法並沒有差異；也就是說，使用最小平方法將變異數分割至頻率不同且彼此正交的弦波時即不須計算反矩陣^[31]。由於此類分析可用快速傅立葉變換的方法實作，當訊號完整且取樣間距相等時，便不會使用相對耗時的最小平方頻譜分析法。

參見

• 譜密度

參考資料

1. Cafer Ibanoglu. Variable Stars As Essential Astrophysical Tools. Springer. 2000. ISBN 0-7923-6084-2.
2. D. Scott Birney; David Oesper; Guillermo Gonzalez. Observational Astronomy. Cambridge University Press. 2006. ISBN 0-521-85370-2.
3. Press. Numerical Recipes 3rd. Cambridge University Press. 2007. ISBN 0-521-88068-8.
4. J. Taylor; S. Hamilton. Some tests of the Vaníček Method of spectral analysis. Astrophysics and Space Science. 1972-03-20, 17 (2): 357–367. Bibcode:1972Ap&SS..17..357T. doi:10.1007/BF00642907.
5. Alistair I. Mees. Nonlinear Dynamics and Statistics. Springer. 2001. ISBN 0-8176-4163-7.
6. Frank Chambers. Climate Change: Critical Concepts in the Environment. Routledge. 2002. ISBN 0-415-27858-9.
7. Hans P. A. Van Dongen. Searching for Biological Rhythms: Peak Detection in the Periodogram of Unequally Spaced Data. Journal of Biological Rhythms. 1999, 14 (6): 617–620. PMID 10643760. doi:10.1177/074873099129000984.
8. Lomb, N. R. Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science. 1976, 39 (2): 447–462. Bibcode:1976Ap&SS..39..447L. doi:10.1007/BF00648343.
9. Scargle, J. D. Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data. Astrophysical Journal. 1982, 263: 835. Bibcode:1982ApJ...263..835S. doi:10.1086/160554.
10. David Brunt. The Combination of Observations 2nd. Cambridge University Press. 1931.
11. Barning, F. J. M. The numerical analysis of the light-curve of 12 Lacertae. Bulletin of the Astronomical Institutes of the Netherlands. 1963, 17: 22. Bibcode:1963BAN....17...22B.
12. Pascal Vincent; Yoshua Bengio. Kernel Matching Pursuit (PDF). Machine Learning. 2002, 48: 165–187 [2018-10-21]. doi:10.1023/A:1013955821559. （原始內容存檔 (PDF)於2016-03-03）.
13. Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, "Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition," in Proc. 27th Asilomar Conference on Signals, Systems and Computers, A. Singh, ed., Los Alamitos, CA, USA, IEEE Computer Society Press, 1993.
14. Vaníček, P. Approximate spectral analysis by least-squares fit. Astrophysics and Space Science. 1969, 4 (4): 387–391. Bibcode:1969Ap&SS...4..387V. doi:10.1007/BF00651344.
15. Vaníček, P. Further development and properties of the spectral analysis by least-squares. Astrophysics and Space Science. 1971, 12: 10–33. Bibcode:1971Ap&SS..12...10V. doi:10.1007/BF00656134.
16. Korenberg, M. J. A robust orthogonal algorithm for system identification and time-series analysis. Biological Cybernetics. 1989, 60 (4): 267–276. PMID 2706281. doi:10.1007/BF00204124.
17. S. Chen and D.L. Donoho (1994), "Basis Pursuit." Technical Report, Department of Statistics, Stanford University, Available at [1] （頁面存檔備份，存於網際網路檔案館）
18. Wells, D.E., P. Vaníček, S. Pagiatakis, 1985. Least-squares spectral analysis revisited. Department of Surveying Engineering Technical Report 84, University of New Brunswick, Fredericton, 68 pages, Available at [2] （頁面存檔備份，存於網際網路檔案館）.
19. Craymer, M.R., The Least Squares Spectrum, Its Inverse Transform and Autocorrelation Function: Theory and Some Applications in Geodesy^{[永久失效連結]}, Ph.D. Dissertation, University of Toronto, Canada (1998).
20. William J. Emery; Richard E. Thomson. Data Analysis Methods in Physical Oceanography. Elsevier. 2001. ISBN 0-444-50756-6.
21. M. Zechmeister; M. Kürster. The generalised Lomb–Scargle periodogram. A new formalism for the floating-mean and Keplerian periodograms. Astronomy & Astrophysics. March 2009, 496 (2): 577–584 [2018-10-21]. Bibcode:2009A&A...496..577Z. arXiv:0901.2573 . doi:10.1051/0004-6361:200811296. （原始內容存檔於2021-01-20）.
22. Andrew Cumming; Geoffrey W. Marcy; R. Paul Butler. The Lick Planet Search: Detectability and Mass Thresholds. The Astrophysical Journal. December 1999, 526 (2): 890–915 [2018-10-21]. Bibcode:1999ApJ...526..890C. arXiv:astro-ph/9906466 . doi:10.1086/308020. （原始內容存檔於2021-07-11）.
23. Korenberg, Michael J.; Brenan, Colin J. H.; Hunter, Ian W. Raman Spectral Estimation via Fast Orthogonal Search. The Analyst. 1997, 122 (9): 879–882. Bibcode:1997Ana...122..879K. doi:10.1039/a700902j.
24. Palmer, David M. A Fast Chi-squared Technique For Period Search of Irregularly Sampled Data. The Astrophysical Journal: 496–502. Bibcode:2009ApJ...695..496P. arXiv:0901.1913 . doi:10.1088/0004-637X/695/1/496.
25. David Palmer: The Fast Chi-squared Period Search. [2018-10-21]. （原始內容存檔於2021-03-18）.
26. Beard, A.G., Williams, P.J.S., Mitchell, N.J. & Muller, H.G. A special climatology of planetary waves and tidal variability, J Atm. Solar-Ter. Phys. 63 (09), p.801–811 (2001).
27. Pagiatakis, S. Stochastic significance of peaks in the least-squares spectrum, J of Geodesy 73, p.67-78 (1999).
28. Steeves, R.R. A statistical test for significance of peaks in the least squares spectrum, Collected Papers of the Geodetic Survey, Department of Energy, Mines and Resources, Surveys and Mapping, Ottawa, Canada, p.149-166 (1981)
29. Richard A. Muller; Gordon J. MacDonald. Ice Ages and Astronomical Causes: Data, Spectral Analysis, and Mechanisms. Springer. 2000. ISBN 3-540-43779-7.
30. Timothy A. Davis; Kermit Sigmon. MATLAB Primer. CRC Press. 2005. ISBN 1-58488-523-8.
31. Darrell Williamson. Discrete-Time Signal Processing: An Algebraic Approach. Springer. 1999. ISBN 1-85233-161-5.

外部連結

• LSSA software freeware download^{[永久失效連結]} (via ftp), FORTRAN, Vaníček's method, from the Natural Resources Canada.