最小平方频谱分析法

长久以来，科学界对于傅立叶分析、周期图法与正弦曲线之最小平方拟合间的密切关联早有所知^[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]。由于此类分析可用快速傅立叶变换的方法实作，当讯号完整且取样间距相等时，便不会使用相对耗时的最小平方频谱分析法。

• 谱密度

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