非均匀采样条件下残周期正弦波形的最小二乘拟合算法

梁志国

(北京长城计量测试技术研究所 计量与校准技术重点实验室，北京 100095)

1 引 言

正弦波形四参数曲线拟合是测量测试中的基本算法，具有广泛的应用前景和工程价值，因而，有众多学者对该问题进行了广泛而深入的研究[1～25]。而残周期正弦波曲线拟合，则在超低频参数的快速估计、调制解调等方面具有特别的价值和意义，其最主要的困难在于残周期条件下，正弦波拟合参数的初始值很难获得，因而在四参数非线性迭代的传统算法中不易实现最小二乘拟合。直到文献[11]采用了三参数拟合结合一维搜索迭代方式解决了该问题，才最终获得了残周期均匀采样条件下的最小二乘拟合结果。

本文主要内容，是针对采样间隔不恒定的非均匀采样情况下，残周期正弦波形的四参数曲线拟合开展研究，以解决非均匀采样条件下残周期正弦波形的曲线拟合问题。它将可以用于均匀采样、非均匀采样、随机采样等同时提供时间坐标和幅度信息的采样波形序列的正弦波形拟合。

2 三参数正弦波形的最小二乘拟合算法

2.1 三参数正弦曲线拟合法

设理想正弦波形为：

y(t)=E1cos(2 π ft)+E2sin(2 π ft)+Q

=Ecos(2 π ft+φ)+Q

(1)

数据记录序列为已知时刻t1,t2,…,tn的采集样本y1,y2,…,yn。三参数正弦曲线拟合过程，即为输入信号的频率f已知，选取或寻找A1、B1、C，使下式所述残差平方和ε最小：

(2)

由式(2)的ε最小，可得拟合函数：

=Acos(2 π fti+θ)+C

(3)

(4)

(5)

拟合残差有效值为：

由于这是一种闭合算法，因而收敛是肯定的。

2.2 问题讨论

上述三参数正弦曲线拟合过程，是已知信号频率f下进行的。

特别地，令tmin=min{ti}；tmax=max{ti}；则，平均采样速率v=(n-1)/(tmax-tmin)；数字角频率 ω=2 π f/v。

在均匀采样条件下，式(1)可表示成下列离散形式：

y(i)=E1cos(ω·i)+E2sin(ω·i)+Q

=Ecos(ω·i+φ)+Q

(6)

则拟合函数将是：

=Acos(ω·i+θ)+C

(7)

不失一般性，设给定信号测量序列的数字角频率为w，而不是ω，则，使用1/3个周期的正弦波形序列，按上述三参数正弦曲线拟合得如图1所示归一化误差ρ/E与频率比w/ω关系曲线波形。图2为图1的局部细化。

图1 1/3周期波形拟合时归一化误差ρ/E与频率比w/ρ关系曲线Fig.1 Variety of error ρ/E via frequency ratio w/ω when 1/3 period

图2 p/E≤0.05时图1曲线的局部细节Fig.2 Details of Fig.1 in the range of p/E≤0.05

其中ω=2 π/10 000；E=4；Q=0；φ分别取0°、35°、70°、105°、140°、175°、210°、245°、280°、315°、350°；以w代替ω进行三参数拟合。

从图1及其局部细化曲线波形图2可见，对于信号频率ω的不足一个信号周期的残周期正弦曲线波形拟合，有如下规律：

(1) 三参数最小二乘拟合法在(0,2ω]范围内ρ/E极小值存在且唯一，极值点就是ω频率点。

(2) 幅度和相位的变化基本不影响ρ/E的变化规律；但初始相位φ的变化将影响ρ/E的幅度值。

(3) 当使用的拟合频率w大于信号频率的2倍时，ρ/E的量值均普遍大于拟合频率w小于信号频率的情况。该规律可用于判定收敛区间上界，而收敛区间的下界则可由足够接近0频的一个微小频率值代替。

(4) 可在(0,2ω]范围内通过对ω的一维搜索找出ρ/E的极小值点。

使用其它部分周期波形曲线进行拟合可以获得与上述规律相同的结论。

3 四参数正弦波形的最小二乘拟合算法

3.1 原理过程

对上述三参数正弦曲线拟合方法的改造，可获得一种绝对收敛的残周期正弦波形四参数曲线拟合方法。

假设平均采集速率为v，待估计的正弦波频率目标值为f0，待估计的正弦波采样序列所含信号不足一个周期，个数为p(0q/τ；因而，可以肯定f0∈[q/τ,2/τ]，在区间[q/τ,2/τ]内的任意频率f下，残差平方和ε(f)的极值存在且唯一！这样，便将四参数正弦波曲线拟合中，对幅度、频率、相位、直流分量4个参数的四维非线性搜索，变成了对频率分量f造成的ε(f)的一维线性搜索，可保证在区间[q/τ,2/τ]内，用三参数拟合法实现的四参数正弦曲线拟合过程收敛。该四参数拟合过程如下：

(1)设定拟合迭代停止条件为he。

(2)从已知时刻t1,t2,…,tn的正弦波采集样本y1,y2,…,yn。使用计点法获得信号波形占用时间长度为τ=tn-t1；平均采集速率v=(n-1)/(tn-t1)，选取因子q(例如q=10-5)，确定目标频率f0的存在区间[q/τ,2/τ]。

(3)确定迭代左边界频率：fL=q/τ；迭代右边界频率：fR=2/τ。

(4)令中值频率：fM=(fR+fL)/2；在左边界频率和右边界频率和中值频率上分别利用3参数拟合公式计算各自的拟合残差ρ(fL)、ρ(fM)和ρ(fR)。

(5)判断是否ρ(fL)<η·ρ(fM)？其中，η为判据因子，取值范围为1～1.5。

若ρ(fL)<η·ρ(fM)，则令fR=fM，fL=不变，重复执行(4)～(5)的过程。

(6)若ρ(fL)≥η·ρ(fM)，则必有fR<2f0，确定迭代左边界频率为fL；迭代右边界频率fR；中值频率：fM=fL+0.618×(fR-fL)；fT=fR-0.618×(fR-fL)。

(7)在fL上执行三参数正弦曲线拟合，获得AL、θL、CL、ρL；在fR上执行三参数正弦曲线拟合，获得AR、θR、CR、ρR；在fM上执行三参数正弦曲线拟合，获得AM、θM、CM、ρM；在fT上执行三参数正弦曲线拟合，获得AT、θT、CT、ρT。

(8)若ρM<ρT，则ρ=ρM，有f0∈[fT,fR]，fL=fT，fT=fM；fM=fL+0.618×(fR-fL)；若ρM>ρT，则ρ=ρT，有f0∈[fL,fM]，fR=fM，fM=fT；fT=fR-0.618×(fR-fL)。

(9)判定是否|(ρM(k)-ρT(k))/ρT(k)|

3.2 收敛性问题

四参数正弦曲线拟合是一个迭代过程，也有收敛性问题。如上所述，对于含p(p⦤1)个信号周期的测量序列的情况，如图1所示，四参数正弦曲线拟合的收敛区间为(0,2ω0]；其中，ω0为信号的真实数字角频率。

4 仿真实验验证

选取测量范围-5～+5 V，幅度4 V，直流分量0 V，初始相位100°，频率6 Hz，均匀采样速率8 kSa/s，采样间隔τ0=0.125 ms，采样序列长度n=400的采样序列，约含0.3个波形周期。使用本文上述四参数拟合算法进行正弦参数估计，获得拟合曲线与采样曲线波形如图3所示，两者之差异曲线如图4所示。其拟合参数如表1所示。

表1 残周期正弦序列参数拟合结果Tab.1 Results of partial period sine wave curve-fitting

图3 残周期拟合曲线与均匀采样曲线yiFig.3 Uniform sampling points yi and curve-fitting

图4 残周期拟合曲线与均匀采样曲线差异Fig.4 difference between uniform sampling series and

令正弦波模型参数不变，使用随机采样间隔进行序列采样，获得采样序列并按照本文上述方法进行正弦波曲线拟合，得采样序列曲线及其拟合曲线分别如图5～图8所示。其中，图5为波形长度约为0.4个周期的非均匀采样序列曲线及其拟合曲线，其采样间隔为(RND-0.4)×20τ0。其中，RND为在[0,1]内均匀分布的随机数。可见由于采样间隔呈非均匀随机变化状态，按照等间隔均匀采样模式绘制的曲线图5，其中的“正弦波形”已经变化非常大。

图6为拟合曲线与采样序列之差异曲线。可见两者拟合得非常好。相应拟合参数如表1所示。由表1数据可见，非均匀采样条件下正弦曲线拟合误差要略大于均匀采样条件下，但差异不太大，造成差异的原因有待于进一步研究。图7为波形长度约为1.2个周期的非均匀采样序列曲线及其拟合曲线，其采样间隔为(RND-0.3)×20τ0。

图6 残周期拟合曲线与非均匀采样曲线差异Fig.6 Difference between non-uniform sampling series and

图7 拟合曲线与非均匀采样曲线yiFig.7 Non-uniform sampling points yi and curve-fitting

图8为拟合曲线与采样序列之差异曲线，拟合情况良好，相应拟合参数如表1所示。

图8 拟合曲线与非均匀采样曲线差异Fig.8 Difference between non-uniform sampling series and

5 实验验证

用超低频振动标准装置作激励源，给出位移幅度36.32 mm的正弦振动，频率0.040 000 Hz，用ASQ-1CA型位移传感器进行测量，以NI PXI-6281型数据采集系统执行采集，其A/D位数18 Bit，工作量程±2.5 V，采样速率200 Sa/s，采集数据个数 5 000 点。为含有约1个周期波形的振动序列。

使用上述四参数正弦波拟合方法获得其振动波形幅度为1 771.02 mV，频率为0.040 56 Hz，相位为230.098°，直流分量为120.53 mV。残差有效值ρ=30.157 mV。波形总失真度为2.4%。其测量序列波形及拟合曲线如图9所示，测量序列波形与拟合曲线之差如图10所示。

图9 拟合曲线与均匀采样曲线yiFig.9 Uniform sampling points yi and curve-fitting

图10 拟合曲线与均匀采样曲线差异Fig.10 Difference between uniform sampling series and

将上述测量曲线波形截取3段后，变成非均匀采样序列，使用上述四参数正弦波拟合方法获得其振动波形幅度为1 746.93 mV，频率为0.040 42 Hz，相位为231.354°，直流分量为123.94 mV。残差有效值ρ=34.154 mV，波形总失真度为2.8%。其测量序列波形及拟合曲线如图11所示，测量序列波形与拟合曲线之差如图12所示。

图11 拟合曲线与非均匀采样曲线yiFig.11 Non-uniform sampling points yi and curve-fitting

图12 拟合曲线与非均匀采样曲线差异Fig.12 Difference between non-uniform sampling series and

6 讨 论

从上述仿真和实际实验结果可见，本文所述基于非均匀采样条件的残周期正弦波拟合的测量方法，可以用于不足一个波形周期的残周期非均匀采样条件下正弦波形参数的测量估计，并可以给出幅度、频率、相位、直流分量等基本信息参量。目标是处理残周期条件下的正弦波拟合，实际上可以在2个波形周期以下的情况均可直接使用本文方法处理。其最大的特点是针对残周期正弦波模型，不论其采样间隔是否均匀，均可以获得有效收敛结果，特别是针对均匀采样序列剔出一段异常波形情况，在实际工作中时有发生，而本来希望获得均匀采样序列，但由于周期过长以及测量过程不完善等因素，使得采样序列变成非均匀采样序列的情况也比比皆是。对于随机采样情况，甚至时间排序出现前后颠倒的情况下，本文方法依然可以有效获得拟合结果。体现出算法良好的收敛性和鲁棒性，可直接用于解决上述范畴内的工程技术问题。

仿真实验表明，与均匀采样条件相比，非均匀采样序列正弦拟合误差要略大一些，在本文实验中，幅度差异约在0.2%，频率差异约在0.5%，相位差异约在1.5°，直流分量差异在0.3%(相对于正弦幅度)。实际实验情况也验证了这些。

此前的研究已经表明，与多周期情况下的测量参数相比，残周期参数拟合存在更大的误差，特别是有较大直流分量情况下。由于信息不全，波形失真、序列样本长度、波形周期长度等众多因素都将影响参数估计效果，因而对于波形序列长度要求也是不一样的，研究表明，在通常的波形质量条件下，1/5以上周期的波形长度即可以进行参数的正确估计。

7 结 论

综上所述可见，本文主要是提出了非均匀采样条件下残周期正弦波的四参数拟合的一种方法，给出了详细过程和收敛区间。分别在随机间隔采样和非随机间隔采样两种条件下给出了比较结果，结果表明，本文所述方法具有良好的收敛性和鲁棒性，仅要求已知采样序列和每个采样点的时刻，没有任何其它先决条件要求，是一种表现良好的正弦参数拟合方法，可直接用于非均匀采样条件下残周期正弦波形的参数拟合，对于超低频参数估计尤其具有特别的意义和价值。