正弦波信号序列的最小二乘拟合算法

理想正弦信号可用下述四参数表达式表示：
$$y(t)=A_0\cos ({\omega }_{0}t)+B_0\sin ({\omega }_{0}t)+C_0$$
或：
$$y(t)=A_0\cos ({\omega }_{0}t+\theta )+C$$
用指定参数的正弦波作为数字存储示波器的输人，得到一组数据记录。通过改变拟合正弦信号的相位、幅度、直流偏移和频率，使拟合结果和输人记录序列各点的残差平方和最小。这里提供两种正弦拟合方法：一种用于采样频率和被测信号输入频率均已知时；另一种用于信号频率未知时。第一个途径包括两种算法:一种通过矩阵运算，另一种通过迭代过程。对于已知信号频率的情况，当初始条件相同时，上述两个算法结果一致。但两者和收敛性不一样，使用矩阵算法的收敛速度比不使用矩阵的快，特别是信号周期数小于$5$个时。

注： 下列$1$中所述的三参数算法(对已知频率)是一种闭合算法。因此，总能获得一个结果。但是，如果算法中使用的频率(假设已知)和实际输入的频率不一样，或采集速率有较大误差，三参数算法的结果比稍后所述的四参数算法结果要差。但$2$中的四参数算法在初始条件偏离较多，或有一些特别不正确的数据下，迭代过程可能发散。

1 正弦波三参数(已知频率)最小二乘拟合算法

1.1 正弦信号三参数最小二乘拟合——矩阵算法

设理想正弦信号为：
$$y(t)=A_0\cos ({\omega }_{0}t)+B_0\sin ({\omega }_{0}t)+C_0$$
数据记录序列为时刻$t_1,t_2,...,t_m$的采集样本$Y_1,Y_2,...,Y_m$，拟合过程即为选取或寻找$A_0,B_0,C_0$，使下式所述残差平方和最小：
$$\begin{array}{c}\sum _{n=1}^M[Y_n-A_0\cos ({\omega }_{0}t)-B_0\sin ({\omega }_{0}t)-C_0{]}^2\end{array}$$
式中${\omega }_0$是数字存储示波器输人信号频率(假设已知)。
为了找出合适的$A_0,B_0$和$C_0$值，首先构造下列矩阵：
$$D_0=\left[\begin{array}{ccc}\cos ({\omega }_0{t}_1)& \sin ({\omega }_0{t}_1)& 1\\ \cos ({\omega }_0{t}_2)& \sin ({\omega }_0{t}_2)& 1\\ ⋮& ⋮& ⋮\\ \cos ({\omega }_0{t}_M)& \sin ({\omega }_0{t}_M)& 1\end{array}\right]y=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_M\end{array}\right]x_0=\left[\begin{array}{c}{A}_0\\ B_0\\ C_0\end{array}\right]$$
式$(1)$可用矩阵方式表示如下
$$\begin{array}{c}(y-D_0{x}_0{)}^T(y-D_0{x}_0)\end{array}$$
这里 $(\ast {)}^T$*表示矩阵 $(\ast )$*的转置，可以得出式 $(2)$最小时的最小二乘解$\widehat{x_0}$

$$\widehat{x_0}=(D_0^T{D}_0{)}^{-1}(D_0^{T}y)$$
拟合函数为：
$$y_n^{\prime }=A_0\cos ({\omega }_0{t}_n)+B_0\sin ({\omega }_0{t}_n)+C_0$$
将其转换为幅度和相位表达形式:
$$y_n^{\prime }=A\cos ({\omega }_0{t}_n+\theta )+C$$
式 中：
$$\begin{gathered} A=\sqrt[]{A_0^2+B_0^2} \\ \theta ={\tan }^{-1}[\frac{-B_0}{A_0}],A_0\ge 0 \\ \theta ={\tan }^{-1}[\frac{-B_0}{A_0}]+\pi ,A_0<0 \end{gathered}$$
拟合残差 $r_n$，如下：
$$r_n=y_n-A_0\cos ({\omega }_0{t}_n)-B_0\sin ({\omega }_0{t}_n)-C_0$$
拟合残差有效值如下：
$${ϵ}_{rms}=\sqrt[]{\frac{1}{M}\sum _{n=1}^M{r}_n^2}$$

1.2 正弦信号三参数最小二乘拟合——代数方程算法

设数据记录序列为包含了一系列在时刻$t_M$上采集的正弦信号样本$y_M$，数字存储示波器输人信号的角频率为$\omega$，定义:

$$\begin{gathered} {\alpha }_n=\cos (\omega t_n) \\ {\beta }_n=\sin (\omega t_n) \end{gathered}$$
然后计算下面9个和：
$$\begin{gathered} \sum _{n=1}^M{y}_n\sum _{n=1}^M{\alpha }_n\sum _{n=1}^M{\beta }_n \\ \sum _{n=1}^M{\alpha }_n{\beta }_n\sum _{n=1}^M{\alpha }_n^2\sum _{n=1}^M{\beta }_n^2 \\ \sum _{n=1}^M{y}_n{\alpha }_n\sum _{n=1}^M{y}_n{\beta }_n\sum _{n=1}^M{y}_n^2 \end{gathered}$$
使用这些和计算：
$$\begin{gathered} A_1=\frac{A_N}{A_D} \\ {B}_1=\frac{B_N}{B_D} \\ C=\bar{y}-A_1\bar{\alpha }-B_1\bar{\beta } \end{gathered}$$
式中：
$$\begin{gathered} A_N=\frac{{\sum }_{n=1}^M{y}_n{\alpha }_n-\bar{y}{\sum }_{n=1}^M{\alpha }_n}{{\sum }_{n=1}^M{\alpha }_n{\beta }_n-\bar{\beta }{\sum }_{n=1}^M{\alpha }_n}-\frac{{\sum }_{n=1}^M{y}_n{\beta }_n-\bar{y}{\sum }_{n=1}^M{\beta }_n}{{\sum }_{n=1}^M{\beta }_n^2-\bar{\beta }{\sum }_{n=1}^M{\beta }_n} \\ {A}_D=\frac{{\sum }_{n=1}^M{\alpha }_n^2-\bar{\alpha }{\sum }_{n=1}^M{\alpha }_n}{{\sum }_{n=1}^M{\alpha }_n{\beta }_n-\bar{\beta }{\sum }_{n=1}^M{\alpha }_n}-\frac{{\sum }_{n=1}^M{\alpha }_n{\beta }_n-\bar{\alpha }{\sum }_{n=1}^M{\beta }_n}{{\sum }_{n=1}^M{\beta }_n^2-\bar{\beta }{\sum }_{n=1}^M{\beta }_n} \\ {B}_N=\frac{{\sum }_{n=1}^M{y}_n{\alpha }_n-\bar{y}{\sum }_{n=1}^M{\alpha }_n}{{\sum }_{n=1}^M{\alpha }_n^2-\bar{\alpha }{\sum }_{n=1}^M{\alpha }_n}-\frac{{\sum }_{n=1}^M{y}_n{\beta }_n-\bar{y}{\sum }_{n=1}^M{\beta }_n}{{\sum }_{n=1}^M{\alpha }_n{\beta }_n-\bar{\alpha }{\sum }_{n=1}^M{\beta }_n} \\ {B}_D=\frac{{\sum }_{n=1}^M{\alpha }_n{\beta }_n-\bar{\beta }{\sum }_{n=1}^M{\alpha }_n}{{\sum }_{n=1}^M{\alpha }_n^2-\bar{\alpha }{\sum }_{n=1}^M{\alpha }_n}-\frac{{\sum }_{n=1}^M{\beta }_n^2-\bar{\beta }{\sum }_{n=1}^M{\beta }_n}{{\sum }_{n=1}^M{\alpha }_n{\beta }_n-\bar{\alpha }{\sum }_{n=1}^M{\beta }_n} \\ \bar{y}=\frac{1}{M}\sum _{n=1}^M{y}_n\bar{\alpha }=\frac{1}{M}\sum _{n=1}^M{\alpha }_n\bar{\beta }=\frac{1}{M}\sum _{n=1}^M{\beta }_n \end{gathered}$$
拟合函数如下：
$$y_n^{\prime }=A_1\cos (\omega t_n)+B_1\sin (\omega t_n)+C$$
其幅度和相位表达形式：
$$y_n^{\prime }=A\cos (\omega t_n+\theta )+C$$
式中：
$$\begin{gathered} A=\sqrt[]{A_1^2+B_1^2} \\ \theta ={\tan }^{-1}[\frac{-B_1}{A_1}],A_1\ge 0 \\ \theta ={\tan }^{-1}[\frac{-B_1}{A_1}]+\pi ,A_1<0 \end{gathered}$$
拟合残差有效值为：
$${ϵ}_{rms}=\sqrt[]{\frac{ϵ}{M}}$$
式中：
$$ϵ=\sum _{n=1}^M{y}_n^2+A_1^2\sum _{n=1}^M{\alpha }_n^2+B_1^2\sum _{n=1}^M{\beta }_n^2+M{C}^2-2A_1\sum _{n=1}^M{\alpha }_n{y}_n-2B_1\sum _{n=1}^M{\beta }_n{y}_n-2C\sum _{n=1}^M{y}_n+2A_1{B}_1\sum _{n=1}^M{\alpha }_n{\beta }_n+2A_{1}C\sum _{n=1}^M{\alpha }_n+2B_{1}C\sum _{n=1}^M{\beta }_n=\sum _{n=1}^M(y_n-A_1{\alpha }_n-B_1{\beta }_n-C{)}^2=\sum _{n=1}^M(y_n-y_n^{\prime }{)}^2$$
由于这是一种闭合算法，因而收敛是肯定的。

2 正弦波四参数最小二乘拟合算法

2.1正弦波信号四参数最小二乘拟合——矩阵迭代算法

设正弦数据记录序列中时刻$t_1,t_2,...,t_m$上采集的样本$y_1,y_2,...,y_m$，可使用迭代过程寻找到$A_i，B_i,C_i$和${\omega }_i$值，使得下式所述残差平方和最小：
$$\begin{array}{c}\sum _{n=1}^M{\left[y_n-A_1\cos \left({\omega }_i{t}_n\right)-B_i\sin \left({\omega }_i{t}_n\right)-C_i\right]}^2\end{array}$$

式中，${\omega }_i$为数字存储示波器输人信号的频率。

其操作步骤如下：

a）设置循环指针$i=0$，对数字存储示波器输人信号的角频率${\omega }_0$作一个初始估计。可以用离散傅里叶变换($DFT$)来计算频率；或者通过数波形过零点个数计算频率；或简单地输人一个测量频率初始值。使用$1$中给定的三参数矩阵算法，进行拟合以确定$A_0,B_0$和$C_0$。
b）设置$i=i+1$作下一次迭代。
c）使用下式获得新的角频率：
$${\omega }_i={\omega }_{i-1}+\Delta {\omega }_{i-1}(当i=1时，\Delta {\omega }_{i-1}=0)$$
d）构造如下矩阵

$$\begin{gathered} y=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_M\end{array}\right] \\ \begin{array}{l}{\mathrm{D}}_{\mathrm{i}}=\left[\begin{array}{cccc}\cos \left({\omega }_i{t}_1\right)& \sin \left({\omega }_i{t}_1\right)& 1& -A_{i-1}{t}_1\sin \left({\omega }_i{t}_1\right)+B_{i-1}{t}_1\cos \left({\omega }_i{t}_1\right)\\ \cos \left({\omega }_i{t}_2\right)& \sin \left({\omega }_i{t}_2\right)& 1& -A_{i-1}{t}_2\sin \left({\omega }_i{t}_2\right)+B_{i-1}{t}_2\cos \left({\omega }_i{t}_2\right)\\ ⋮& ⋮& ⋮& ⋮\\ \cos \left({\omega }_i{t}_M\right)& \sin \left({\omega }_i{t}_M\right)& 1& -A_{i-1}{t}_M\sin \left({\omega }_i{t}_M\right)+B_{i-1}{t}_M\cos \left({\omega }_i{t}_M\right)\end{array}\right]\\ x_i=\left[\begin{array}{l}{A}_i\\ B_i\\ C_i\\ \Delta {\omega }_i\end{array}\right]\end{array} \end{gathered}$$
e）使式(3)达到最小的最小二乘解用矩阵形式表示如下:
$$\widehat{x_i}=(D_i^T{D}_i{)}^{-1}(D_i^{T}y)$$
f）按下式计算幅度$A$和相位$\theta$：
$$y_n^{\prime }=A\cos (\omega t_n+\theta )+C$$

式中：

$$\begin{gathered} A=\sqrt[]{A_i^2+B_i^2} \\ \theta ={\tan }^{-1}[\frac{-B_i}{A_i}],A_i\ge 0 \\ \theta ={\tan }^{-1}[\frac{-B_i}{A_i}]+\pi ,A_i<0 \end{gathered}$$
9）重复步骤b）到f），根据上一次迭代得到的$A_i,B_i,{\omega }_i$的值，继续迭代直到与$A,\omega ,\theta$和$C$的变化小到满足要求。
拟合残差如下:
$$r_n=y_n-A_i\cos ({\omega }_i{t}_n)-B_i\sin ({\omega }_i{t}_n)-C_i$$
误差有效值如下：
$${ϵ}_{rms}=\sqrt[]{\frac{1}{M}\sum _{n=1}^M{r}_n^2}$$

2.2 正弦信号四参数最小二乘拟合——非矩阵迭代算法

对采集记录下来的数据序列估计一个角频率初始值$\omega$和相位初始值$\theta$，这里的相位指的是记录中的第一个点对应的相位。角频率用弧度每秒来表示，可以用离散傅里叶变换($DFT$)或计算序列过零次数得到，也可以直接使用输人信号频率算得。
相位用弧度来表示，可以按$1$中所述三参数法计算得到。或者用下面的公式计算得到:
$$\left[\mathrm{sign}\left(y_2-y_1\right)\right]{\cos }^{-1}\left(\frac{y_1-C}{A}\right)$$
式中：
$$\mathrm{sign}\left(y_2-y_1\right)=\{\begin{array}{l}1,y_2⩾y_1\\ -1,y_2<y_1\end{array}$$
$y_1$是$t=0$的第一个样本点；$y_2$是紧跟着$y_1$的下一个样本点;$C,A$分别是正弦波的直流偏移和幅度。
使用上式计算时，估计波形幅度可以参照如下方法:(1)如果数字存储示波器噪声不大，用记录数据中的最大值和最小值的代数差的一半作为波形幅度。(2)用众数法找最大和最小值，计算波形幅度。估计偏移的方法:(1)可以使用记录中的最大值和最小值之和的一半。(2)偏移估计的另一方法是取整数个周期的数据的平均值。注意，如果在每个周期中点太少，符号$(y_2-y_1)$能给出不正确的情况。特别是在${\cos }^{-1}$值接近$0$或$\pi$时。
设数据记录序列中包含了一系列时刻$t_n$的取样$y_n$，用以估计$\omega$和$\theta$，计算下面的16个和:
$$\begin{array}{llll}{\sum }_{n=1}^M{\alpha }_n& {\sum }_{n=1}^M{\beta }_n& {\sum }_{n=1}^M{y}_n& {\sum }_{n=1}^M{\alpha }_n{y}_n\\ {\sum }_{n=1}^M{\beta }_n{y}_n& {\sum }_{n=1}^M{\alpha }_n{\beta }_n& {\sum }_{n=1}^M{\beta }_n^2& {\sum }_{n=1}^M{\alpha }_n^2\\ {\sum }_{n=1}^M{y}_n^2& {\sum }_{n=1}^M{\beta }_n{t}_n{y}_n& {\sum }_{n=1}^M{\alpha }_n{\beta }_n{t}_n& {\sum }_{n=1}^M{\beta }_n^2{t}_n\\ {\sum }_{n=1}^M{\beta }_n^2{t}_n^2& {\sum }_{n=1}^M{\beta }_n{t}_n& {\sum }_{n=1}^M{\alpha }_n{t}_n& {\sum }_{n=1}^M{\beta }_n{t}_n^2\end{array}$$
式中：
$$\begin{array}{l}{\alpha }_n=\cos \left(\omega t_n+\theta \right)\\ {\beta }_n=\sin \left(\omega t_n+\theta \right)\end{array}$$
现在使用$\omega$和$\theta$的估计值，计算：
$$\begin{array}{rl}\Psi & =\omega +\frac{a_{22}R-a_{12}S}{a_{11}{a}_{22}-a_{12}{a}_{21}}\\ \varphi & =\theta +\frac{a_{11}S-a_{21}R}{a_{11}{a}_{22}-a_{12}{a}_{21}}\end{array}$$
式中：
$$\begin{array}{l}{a}_{11}=\frac{{\sum }_{n=1}^M{\beta }_n{t}_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\alpha }_n{t}_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}-\frac{{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\beta }_n{t}_n^2\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}\\ a_{12}=\frac{{\sum }_{n=1}^M{\beta }_n{t}_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\alpha }_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}-\frac{{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\beta }_n{t}_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}\\ a_{21}=\frac{{\sum }_{n=1}^M{\beta }_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\alpha }_n{t}_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}-\frac{{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\beta }_n{t}_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}\\ a_{22}=\frac{{\sum }_{n=1}^M{\beta }_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\alpha }_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}-\frac{{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right){\sum }_{n=1}^M{\beta }_n\left({\beta }_n-\bar{\beta }\right)}{{\left[{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)\right]}^2}\\ R=\frac{{\sum }_{n=1}^M{\beta }_n{t}_n\left(y_n-\bar{y}\right)}{{\sum }_{n=1}^M{\alpha }_n\left(y_n-\bar{y}\right)}-\frac{{\sum }_{n=1}^M{t}_n{\beta }_n\left({\alpha }_n-\bar{\alpha }\right)}{{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)}\\ S=\frac{{\sum }_{n=1}^M{\beta }_n\left(y_n-\bar{y}\right)}{{\sum }_{n=1}^M{\alpha }_n\left(y_n-\bar{y}\right)}-\frac{{\sum }_{n=1}^M{\beta }_n\left({\alpha }_n-\bar{\alpha }\right)}{{\sum }_{n=1}^M{\alpha }_n\left({\alpha }_n-\bar{\alpha }\right)}\\ \bar{y}=\frac{1}{M}{\sum }_{n=1}^M{y}_n\\ \bar{\alpha }=\frac{1}{M}{\sum }_{n=1}^M{\alpha }_n\\ \bar{\beta }=\frac{1}{M}{\sum }_{n=1}^M{\beta }_n\\ {\alpha }_n=\cos \left(\omega t_n+\theta \right)\\ {\beta }_n=\sin \left(\omega t_n+\theta \right)\end{array}$$
使用$\Psi$和$\varphi$作为新的$\omega$和$\theta$估计值，重复上述过程，直到两者的差别小到满足要求，产生其拟合函数的形式如下:
$$y^{\prime }{}_n=A\cos \left(\Psi t_n+\varphi \right)+C$$
可按下列式子计算
$$\begin{array}{rl}A& =\frac{{\sum }_{n=1}^M{\beta }_n\left(y_n-\bar{y}\right)\left({\alpha }_n+{\beta }_n+{\beta }_n{t}_n\right)}{{\sum }_{n=1}^M\left({\alpha }_n-\ddot{\alpha }\right)\left({\alpha }_n+{\beta }_n+{\beta }_n{t}_n\right)}\\ C& =\bar{y}-A\bar{\alpha }\end{array}$$
残差有效值为$(\frac{ϵ}{M}{)}^{\frac{1}{2}}$，这里
$$\begin{array}{rl}\frac{\epsilon }{M}& =\frac{1}{M}\sum _{n=1}^M{y}_n^2+\frac{A^2}{M}\sum _{n=1}^M{\alpha }_n^2-2\frac{A}{M}\sum _{n=1}^M{\alpha }_n{y}_n+C^2-2C\bar{y}+2AC\bar{\alpha }\\ & =\frac{1}{M}\sum _{n=1}^M{\left(y_n-A{x}_n-C\right)}^2=\frac{1}{M}\sum _{n=1}^M{\left(y_n-y_n^{\prime }\right)}^2\end{array}$$
由于这是一个迭代过程，所以对一些误差很大的估计初始值$\omega$和$\theta$有可能引起发散。