四参数正弦波组合式拟合算法

梁志国

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

1 引 言

正弦现象是自然界中广泛存在的一种物理现象,因而四参数正弦拟合算法是一种非常重要的参数估计方法,在计量测试行业中广为应用，并有大量先期研究成果[1～13]。时至今日,能够获得最高拟合精度的,仍然是几种四参数迭代搜索法,其它方法的拟合精度要略微差些[14～26]。由于正弦曲线的周期性特征,导致其拟合参数变化时,会出现众多周期性局域极值点,使得拟合容易收敛到局域极值点而非总体最优点上。因此,四参数迭代搜索法存在两点主要不足,其一是需要对4个拟合参数进行预先估计,且估计精度要足够高,才能保证迭代过程收敛。客观上需要用于拟合的正弦采样序列含有一个以上的波形周期。当波形少于一个完整周期时,将不能确保给出4个参数的预估计值,导致无法进行后续拟合。另外,四参数迭代搜索法的收敛性一直存在问题，它何时收敛并无明确结论和判据,也没有明确的收敛区间。只是在实用中,存在初始预估值影响收敛性的现象。初始值距离最优值越近,越能保证迭代过程收敛，否则,将越容易导致迭代过程发散。

在频率已知的情况下,三参数正弦曲线拟合是一种解析算法,无须迭代过程,因而,不存在收敛问题。在此基础上演化而来的四参数拟合算法是近年来新出现的算法。它将4个参数的非线性搜索过程转化为一个参数的单调搜索过程,可称为频率搜索法。其优点是无须对4个参数进行预估计,有明确的收敛域，可用来估计不足一个波形周期的正弦波参数；其缺点是拟合精度不如四参数迭代搜索法。

本文后续内容,将使用两种方法相结合,构造一种新的组合式四参数正弦拟合法,以取长补短,获得优良的拟合效果,并具有明确的收敛域；既适应不足一个波形周期的残周期状况,也适应含多个波形周期的多周期状况，具有四参数迭代搜索法的拟合精度。

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.5个周期的正弦波形序列,按上述三参数正弦拟合法,获得如图1所示归一化误差ρ/E与频率比w/ω关系曲线波形。图2为图1的局部细化。

其中：f=1 Hz；v=1 kSa/s；ω=2 π /1 000；E=1 V；Q=0；φ分别取0°、35°、70°、105°、140°、175°、210°、245°、280°、315°、350°；采样点数n=1 500,以w代替ω进行三参数拟合。

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

图2 1.5个周期波形拟合时归一化误差ρ/E与频率比w/ω关系曲线(局部展开)Fig.2 Variety of error ρ/E vs. frequency ratio w/ω when 1.5 periods (partial)

由图1及图2可见,对于信号频率为ω的正弦曲线,三参数最小二乘拟合法的最优频率在频率区间(0.5ω, 1.5ω]范围内存在且唯一,在该区间内,拟合残差有效值ρ的极值存在且唯一；幅度和相位的变化不影响ρ/E的变化规律以及ρ各个极值点出现的位置,可在区间(0, 1.5ω]范围内通过对ω的一维搜索找出ρ最优值对应的正弦频率极值点。

仿真实验表明,1.5个周期以下的波形序列所遵循的规律和特征与图1相同。在全频率范围内,幅度E的变化基本上不影响ρ/E随w/ω的变化时ρ/E的极值点出现的位置和幅度；相位的变化也不影响ρ/E随w/ω变化时ρ/E极值点出现的位置,而只改变非ω频率处各个ρ/E极值的大小。

使用7个周期的正弦波形序列,按上述三参数正弦拟合法,获得如图3所示归一化误差ρ/E与频率比w/ω关系曲线波形。图4为图3的局部细化。

其中：f=1 Hz；v=200 Sa/s；ω=2 π /200；E=1 V；Q=0；φ分别取0°，35°，70°，105°，140°，175°，210°，245°，280°，315°，350°；采样点数n=1 400,以w代替ω进行三参数拟合。

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

图4 7个周期波形拟合时归一化误差ρ/E与频率比w/ω关系曲线(局部展开)Fig.4 Variety of error ρ/E vs. frequency ratio w/ω when 7 periods(partial)

由图3及图4可见,对于含有7个波形周期的信号频率ω的正弦曲线,三参数最小二乘拟合法的最优频率在区间ω±ω/7(即[6ω/7, 8ω/7])范围内存在且唯一,拟合残差有效值ρ的极值存在且唯一,幅度和相位的变化不影响其ρ/E的变化规律,可在ω±ω/7范围内通过对ω的一维搜索,找出该ρ的最优极值对应的最优频率点。

仿真实验表明,m个周期的波形序列所遵循的规律和特征与图3相同。在全频率范围内,幅度和相位的变化对ρ/E的影响规律与7个周期序列时相同。

总结众多仿真曲线波形可见,对于1.5个波形周期以下的正弦曲线拟合,有如下规律：

1) 频率ω变化时,在全频率范围内,三参数最小二乘拟合法获得的拟合残差有效值ρ/E拥有众多等间隔极值点。

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

3) 幅度和相位的变化基本不影响ρ/E的极值出现的位置及其变化规律；但初始相位φ的变化将影响ρ/E的极值幅度值。

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

5) 可在(0, 1.5ω]范围内通过对拟合频率w的一维搜索找出ρ/E的最小值点,从而获得四个拟合参数,完成四参数正弦波拟合；在该区间范围内,拟合过程绝对收敛。

对于m个波形周期的正弦曲线拟合,有如下规律：

1) 拟合频率w变化时,在全频率范围内,三参数最小二乘拟合法获得的拟合残差有效值ρ/E拥有众多等间隔极值点。

2) 三参数最小二乘拟合法最优拟合频率w在区间[(1-1/m)·ω, (1+1/m)·ω]范围内存在且唯一,在该区间范围内,拟合残差有效值ρ/E的极小值存在且唯一,其极小值点(就是对应的最小二乘最优频率点；因而,该区间[(1-1/m)·ω, (1+1/m)·ω]边界可以作为一维频率搜索的区间上下界。

3) 幅度和相位的变化基本不影响ρ/E的极值点出现的位置及其变化规律；但初始相位φ的变化将影响ρ/E的幅度值。

4) 可在[(1-1/m)·ω, (1+1/m)·ω]范围内通过对拟合频率w的一维搜索找出ρ/E的极小值点,从而获得四个拟合参数,完成四参数正弦波拟合；在该区间范围内,拟合过程绝对收敛。

3 四参数正弦拟合法——频率搜索迭代法

3.1 原理过程

3.1.1 多周期情况

对上述三参数正弦曲线拟合方法的改造,可获得一种绝对收敛的四参数正弦曲线拟合方法。

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

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

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

(5) 若ρ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)。

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

3.1.2 残周期情况

当正弦波采样序列包含不足一个波形周期时,使用上述第3.1.1节所述的频率区间估计方法无法获得有效的搜索区间。此时,使用过零检测法获得的采集序列包含的过零点少于3个,实际上最多可能包含的波形数p将少于1.5个周期,多数情况下均将少于1个波形周期,即0

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

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

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

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

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

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

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

(6) 若ρ(fL)η·ρ(fM), 则必有fR<1.5f0, 确定迭代左边界频率为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(0

对于含p个信号周期的测量序列的情况,如图3所示,四参数正弦曲线拟合的绝对收敛区间为[ω0(1-1/p),ω0(1+1/p)]，其中,ω0为信号的真实数字角频率。

在该绝对收敛区间之外,当p≥2时,首先,可以对测量序列的有效值Erms进行预估计：

不失一般性,假设正弦序列的噪声与信号有效值幅度之比为N/S≪1,即噪声功率远小于信号功率,则可选取判据hd取值满足N/S

然后,变化ω在(0, +∞)区间内搜索从ρ(ω)/Erms>hd变化到ρ(ω)/Erms

当导数ρ′(ωd)<0,ωd<ω0,可令ωL=ωd；

当导数ρ′(ωd)>0,ωd>ω0,可令ωR=ωd。

这样,便寻找出落在绝对收敛区间内且包含ω0的迭代区间[ωL,ωR],使用前述方法可获得绝对收敛的四参数正弦拟合结果。

采取辅助判据hd这样一种措施后,可将本文所述四参数正弦拟合方法的频率绝对收敛区间拓展到(0, +∞),从而可在任何情况下都获得收敛结果。

4 组合式四参数正弦波拟合算法

4.1 算法流程

组合式四参数正弦波拟合算法的流程如下：

1) 针对正弦波形采样序列,使用上述第3节所述的一维搜索迭代正弦拟合算法,获得其幅度、频率、初始相位、直流分量和拟合残差；

2) 将该拟合结果作为四参数搜索迭代正弦拟合算法的拟合初始值,执行文献[2]所述的四参数搜索迭代,以获得最终的四参数正弦波拟合结果,并结束拟合；

3) 若文献[2]所述的四参数搜索迭代不收敛,则仍然以第3节所述的一维搜索迭代获得的正弦拟合结果作为组合式四参数正弦波拟合的最终结果,结束拟合。

4.2 仿真搜索验证

使用幅度E0=1 V；频率f0=1 Hz；直流分量Q0=0；初始相位φ0分别取0°，35°，70°，105°，140°，175°，210°，245°，280°，315°，350°；采样点数n=1000；正弦波形周期数从p=0.001～p=1,取步进Δp=0.001,按照3.1节所述方法进行搜索计算,获得如下图5～图8所示曲线波形。其中：

幅度拟合误差：ΔA/E0=(A-E0)/E0

频率拟合误差：Δf/f0=(f-f0)/f0

相位拟合误差：Δθ=θ-φ0

直流分量拟合误差：ΔC/E0=(C-Q0)/E0

图5 幅度拟合误差ΔA/E0随拟合序列包含周期个数变化情况Fig.5 Amplitude fitting error ΔA/E0 varies with the number of cycles in the fitted sequence

图6 频率拟合误差Δf/f0随拟合序列包含周期个数变化情况Fig.6 Frequency fitting error Δf/f0varies with the number of cycles in the fitted sequence

图7 相位拟合误差Δθ随拟合序列包含周期个数变化情况Fig.7 Phase fitting error Δθ varies with the number of cycles in the fitted sequence

图8 直流分量拟合误差ΔC/E0随拟合序列包含周期个数变化情况Fig.8 DC component fitting error ΔC/E0 varies with the number of cycles in the fitted sequence

从这些图所反映的规律可见,使用本文上述组合式拟合方法,可以在包含很少的正弦周期情况下(最少为千分之五个波形周期)获得5个拟合参数,但不够稳定,拟合误差也较大,初始相位对其也有较大影响。通常,较短的波形在峰值和谷值附近出现时,拟合误差较大,在过零点附近出现时,拟合误差较小。当采样序列包含1/5个周期以上波形时,拟合误差显著降低。此时,幅度拟合误差可降到±2%以内,频率拟合误差可降到±1%以内,相位拟合误差可降到±0.4°以内,直流分量拟合误差与峰值幅度比可以降到±2%以内。

5 实验验证

对上述四参数正弦波曲线拟合算法的实验验证,使用SCO232型数据采集系统采集的数据进行,其A/D位数为12 bits,测量范围为-5 000～+5 000 mV,采集速率为2 000 Sa/s,序列长度为1 800点；信号源为5 700 A型多功能校准器,信号幅度为4 500 mV,频率为11 Hz。

用本文所述方法,获得拟合结果为：信号幅度4 458.311 mV,误差-0.93%；频率11.003 18 Hz,误差2.89×10-4；相位-20.198°,直流分量为-14.668 mV,对应的有效位数6.85 bits,噪信比5.62×10-3。序列包含9.9个波形周期。其拟合曲线与原始数据序列部分值如图9所示,而原始测量序列数据与拟合回归值间的差值如图10所示。

图9 原始数据序列{yi}与拟合曲线Fig.9 Both sampling series {yi} and fitting curve

图10 原始数据序列与拟合曲线之差Fig.10 Difference between sampling series and

在上述采集序列中截取60个采集点,变成不足一个波形周期的残周期状态。使用上述方法进行四参数拟合,获得拟合结果为：信号幅度4 878.354 mV,误差8.41%；频率10.430 19 Hz,误差-5.18%；相位-19.883°,直流分量为-474.981 mV,对应的有效位数7.78 bits,噪信比2.69×10-3。序列包含0.313个波形周期。其拟合曲线与原始数据序列如图11所示,而原始测量序列数据与拟合回归值间的差值如图12所示。

图11 原始数据序列{yi}与拟合曲线Fig.11 Both sampling series {yi} and

图12 原始数据序列与拟合曲线之差Fig.12 Difference between sampling series and

由图9、图11可见,拟合曲线与测量序列基本重合,而由图10、图12的差值曲线可以看出,拟合残差正负方向分布基本均匀,体现了拟合过程的有效性,验证了本文上述方法的实用性。

其中,多周期序列的幅度与频率拟合参量和输入信号标准值的差异很小,而仅仅截取部分周期的残周期序列的幅度与频率拟合参量和输入信号标准值的差异有所增大,符合拟合运算的原理趋势,并且也见证了本文所述方法在实际工作中,可以同时适用于多周期正弦序列和残周期正弦序列的参数拟合。

6 讨 论

通过上述仿真及实测实验可见,使用本文所述的组合式正弦拟合方法,可以将四参数搜索迭代正弦拟合方法拓展到小于一个周期的残周期情况,扩展了其应用范围和使用空间,使得那些由于周期极长,因而数据不全的正弦拟合问题得以有解决之道。另外,由于单参数搜索迭代正弦拟合算法获得的四个参数已经足够精确,并总能给出收敛拟合结果。因此以它们为初始值,再进行四参数搜索迭代时,可以很容易获得收敛结果。若有不收敛的情况出现时,直接以单参数搜索迭代正弦拟合算法获得的四个参数给出拟合结果,从而能够保证本文所述的组合式正弦拟合方法总能获得收敛结果。避免了单纯的四参数搜索算法的收敛区间不明确、并有时会出现不明原因的发散问题。

本文所述方法是在等间隔均匀采样正弦波序列上实现的。实际上,单参数搜索迭代正弦拟合算法与四参数搜索迭代正弦拟合方法均能够适用于非均匀采样序列情况,使得上述方法也能同时推广应用于非均匀采样序列的参数拟合中,从而使本文所述方法有望成为一种性能优良的普适型方法。可用于0.2个波形周期以上的正弦曲线序列的四参数拟合,给出收敛的最优结果。当然,其先决条件是被测波形中不含有谐波；若残周期波形本身含有过高的谐波成分,将使问题复杂化,导致拟合结果的含义存在争议。

本文使用组合式方法获得了收敛性和高精度的双重效果,精度可以达到文献[2]所述四参数迭代的情况,高于频率搜索方法,具体量值可参见文献[22]；而付出的代价是算法的运算量大幅提升,近似逼近频率搜索法和四参数迭代法的时间之和。这是其明显的弱点,有待后续进行快速算法的研究,以降低算法的复杂程度。

7 结 语

综上所述,通过本文提供的方法和流程,可以在任何情况下给出正弦波形采样序列拟合参数。无论是多周期情况,还是不足一个波形周期的残周期情况,拓展了四参数搜索正弦拟合的应用范围；并且,由于组合方法的预估计方法具有明确的收敛域,可保证组合后的方法在该区间内获得收敛的拟合结果。因此,这使得它也可以应用到调制信号的数字化精确解调上[27],而多数拟合方法,因无法保证过程的收敛性,不适于调制信号解调应用。