基于移动正弦拟合的快速解包络算法*

王选择, 洪 潭, 胡晓晖, 翟中生, 王淞略

(1.湖北工业大学 机械工程学院,湖北 武汉 430068；2.湖北省现代制造质量工程重点实验室,湖北 武汉 430068；3.武汉市康达电气有限公司,湖北 武汉 430014)

0 引 言

调幅信号是指载波频率固定，幅值随被测对象状态发生改变的信号，其包络曲线实时反映被测对象的幅值变化。通常包络曲线的获取方法有硬件处理与软件处理两类方法，硬件处理一般采用二极管包络检波电路处理的方式，主要是运用二极管的单向导通性，构造绝对值处理电路来获取调幅信号的绝对大小。这种方法适合信号较大的单边带调幅信号，受噪声影响较大，且元件参数选择不当，容易出现惰性失真和负峰切割失真，一般不应用于数字采样处理的高频载波的调幅信号之中[1,2]。软件处理方法主要包括Hilbert幅值解调法和小波变换法，Hilbert幅 值解调法是通过Hilbert 变换获取信号的幅值和相位，由于随机噪声的干扰，解调获得的包络轮廓较为粗糙[3]。基于线性相位小波具有奇对称实部与偶对称虚部的性质，小波变换法[4,5]利用小波分解并重构一对幅值与频率均相同的正交信号，再对调幅信号的幅值变化实现解调制，小波变换克服了Hilbert变换的不足，提高了包络提取的精度，但需要根据不同的信号选择相应的小波基函数，且不适合嵌入式系统。

正弦拟合的幅相位算法具有很高的精度，但包络计算的逐段拟合耗时长，且要求系统存储空间大，不适合高速实时处理。为此，本文以最小二乘线性正弦拟合算法[6,7]为基础，提出了一种数字式移动正弦拟合快速求解包络曲线

的方法。

具体思路如下：对直接AD采样的数字调幅信号在逐段正弦移动拟合过程中，首先对以起始点数据段进行完整的正弦拟合，然后结合前向递推思想，利用上一次拟合的中间矩阵参数，在不需要重新拟合的条件下，计算下一点开始的数据段拟合结果，并得到新的中间矩阵参数，直到获取移动过程中每一段的拟合结果，最后根据拟合参数与调制信号的关系，设计适合的计算公式，准确求取调幅信号的包络曲线。由于算法充分利用递推思想，极大降低了计算量与包络曲线求解时间，且具有正弦拟合算法抗干扰能力强的特点，同时算法灵活，不仅适合整周期，也适合非整周期采样数据的处理。另外，算法适应性广，对采样信号是否带直流偏置无要求。

1 线性正弦拟合算法

1.1 调幅信号的包络解析与正弦拟合算法

调幅信号是频率固定幅值随时间变化的曲线，假设调幅信号的公式表达为

x(t)=u(t)cos(φ0+2πfct)+c

(1)

式中fc为载波频率，u(t)为调制信号，代表被测对象的幅值，φ0为信号的初始相位，c为信号的直流偏置。调幅信号直接经过满足采样定理的A/D转换后，得到等间隔的离散采样信号

xi=Aicos(φ0+2πifc/fs)+ci

(2)

式中fs为采样频率，一般要求fs>10fc，i为采样序列，Ai为调幅信号幅值。为了得到调幅信号的包络曲线，需求得式(2)中每个采样序列对应的正弦信号的幅值Ai，由于单个数据点不能得到信号幅值和相位，需要选取一定区间内的数据集进行正弦拟合。

根据正弦信号周期性的特点，对整周期采样点数的拟合具有更高的精度，因此，包络信号的提取一般选择一个周期的采样长度进行逐段拟合。事实上，若载波周期不是采样周期的整数倍，为减小计算的复杂度且保证包络信号的准确性，则尽量选取与载波周期理论点数最相近的整数点数。

1.2 线性最小二乘正弦拟合算法思想

考虑到信号u(t)的变化频率远远低于载波频率，调幅信号的幅值变化缓慢，可认为fc在决定周期Tc内幅值与直流偏置均为常数，故在一个较小的采样时间段内，则式(2)可以写成

xi=Acos(φ0+iδ)+c+εi

(3)

式中δ=2πfc/fs，为采样相移量，εi为采样信号的误差。实际逐段拟合计算过程中，可假设每小段采样区间内的数据i从0开始，即i=0,…,nT-1。为得到参数A，φ0，c的大小，按照最小二乘拟合的思想，满足误差平方和最小原则，令有

a=Acosφ0,b=-Asinφ0

(4)

则有线性化最小二乘拟合为

f(a,b,c)=∑(xi-acos (iδ)+bsin(iδ)-c)2

=min

(5)

因此，求解参数A，φ0，c转变成求解最小二乘拟合值参数a，b，c的问题。为使最小二乘拟合误差最小，满足式(2)的条件为

(6)

将式(6)写成矩阵形式，得到等式方程

(7)

令

(8)

则根据式(8)可以得到

(9)

求解出最小二乘拟合值参数a，b，c后，即可得到正弦信号幅值A和初始相位φ0。

一次正弦拟合只能得到调幅信号的一个包络点，为了得到所有包络点，形成一个包络曲线，需要求解出包络曲线上每一点对应的包络幅值，即需要对调幅信号进行逐段正弦拟合，但这会导致运算量较大，极大地增加了包络曲线的求解时间,因而，需要一种能够进行快速拟合的算法。

2 移动正弦拟合快速算法

2.1 调幅信号的包络解析与正弦拟合算法

从线性最小二乘正弦拟合算法的计算过程可以看出，为了得到采样数据某一点的信号幅值，需要对以该点为中心的一段区间数据进行正弦拟合。且随着被拟合数据点的移动，拟合数据区间也对应移动。移动正弦拟合的数据区间移动过程如图1所示,当求解移动数据点T(m)，T(m+1)，T(m+2), …对应的正弦幅值时，需要分别对区间m，m+1，m+2,…内的数据进行正弦拟合。从图中可以看出，在相邻数据区间移动过程中，实际参与拟合的数据只是去掉首端值，尾端增加一个值，总数据量不变。移动拟合过程与移动平均[8]一样，上一个数据拟合区间与下一个数据拟合区间存在大量重复数据点的处理。根据拟合算法的特点，可对重复点的数据进行适当的区分处理，实现移动正弦拟合的快速算法。为此，这里引入向前递推的算法，来降低拟合计算量，提高移动正弦拟合的计算速度。

图1 调幅信号移动正弦拟合区间

2.2 向前递推算法与求解公式

假设当前移动拟合数据区间的数据组Xk=[xkxk+1…xk+nT-1]，则下一组为Xk+1=[xk+1xk+2…xk+nT]。根据式(9)可知，在区间拟合过程中，每一个区间都作为一组新的数据进行正弦拟合，而每组数据正弦拟合的关键在于中间矩阵B与Y的求解。由拟合参数的推导过程可以看出，参与计算的矩阵B仅与拟合数据量以及采样相移量相关，因此在移动正弦拟合过程中，矩阵B与逆矩阵B-1是与采样数据无关的常矩阵，可对其预先初始化。在实际移动拟合过程中，为得到区间长度为nT的拟合区间中间点的幅值与相位，式(8)中的矩阵B改造为如下形式

B=

其中,r=(nT-1)/2，并且由式(8)可知矩阵Y是随着区间数据组的变化而发生改变的。假设数据组Xk对应的矩阵为Yk，数据组Xk+1对应的矩阵为Yk+1，那么分别有

应用三角函数公式的数学推导，可以发现Yk+1与Yk之间存在如下的递推关系

(10)

显然在Yk已知的条件下，此递推公式通过预存cosδ，sinδ，cos((nT-1)δ)与sin((nT-1)δ)四个常量，仅需要少许的计算量即可求得Yk+1的值，极大地降低了矩阵Yk+1的计算工作量，从而大大降低了整体正弦移动拟合的计算工作量。在得到矩阵B-1与Yk的基础上，利用式(9)即可以快速得到拟合值ak，bk，ck。

值得指出的是，由于正弦拟合需要一定的数据量，因此在移动正弦拟合算法同移动平均一样，出现边沿损失现象，所获得的包络曲线数据量比采样数据少一个正弦拟合数据量nT。

2.3 递推算法的运算复杂度分析

结合Yk递推计算量与正弦拟合公式(9)，计算得到递推法与直接拟合的运算量对比图如图2所示。

图2 递推法与直接拟合的运算量对比

图2中横纵坐标均为对数坐标显示，没有采用递推算法的每一次正弦拟合，需要2nT+9次浮点乘法运算，3(nT-1)+6次加法运算；而采用递推算法后仅要求15次乘法运算与14次加法运算，与区间数据量nT大小无关。从预存函数值空间来看，非递推算法要求预存2nT个正余弦函数值，而递推算法仅需4个函数值，所占内存空间少，易于单片机嵌入式系统的实现。由于非递推拟合算法的运算复杂度正比于拟合数据量，递推算法随着nT的增加，速度优势更加明显。假设1 000个采样数据，不同周期采样点的计算工作量比较。如每个周期达到100点以上的条件下，递推移动正弦拟合运算相对于直接拟合算法乘法计算次数缩小1/15，加法次数缩小到1/20，这有利于高频率采样信号的处理，且总体计算量与采样频率无关。

3 包络曲线的快速求解及实例

3.1 包络曲线的求解方法

理论上,包络曲线所求结果是调幅信号的幅值曲线。在得到正弦移动拟合结果ak，bk，ck后，可以按照如下公式(11)，直接计算幅值Ak

(11)

嵌入式系统中是通过迭代运算求取平方根，为了减小计算量，可以采用对拟合值直接求绝对值的方法来计算包络曲线。根据公式(4)，可知正弦参数拟合值得绝对值满足

在载波采样周期nT内，幅值Ak近似满足如下关系

(12)

因此，可直接通过求解拟合参数的绝对值，来求取包络幅值曲线，将平方运算转换成加法与乘法运算，可以显著降低移动拟合后包络曲线的计算量。

如图3所示为某包络信号正弦拟合参数拟合值ak与bk的绝对值及移动平均曲线图。

图3 拟合参数的绝对值及其移动平均结果

图3中所显示的是采样信号第750～850个采样点的处理结果，故结果中不包括移动拟合和移动平均所导致的边沿效应。从图中局部放大的绝对值移动平均曲线来看，两参数的ak,bk绝对值移动平均存在较大的周期波动起伏，进一步观察可以看出，两者的波动大小正好相反，可以互相弥补，两者的绝对值均值的移动平均更能反映实际的包络曲线。故包络幅值的计算公式修正如下

(13)

3.2 超声测距解包络实例

在超声解调测距系统中，控制输出激励信号，返回的超声波回波放大信号如图4所示。

图4 超声波回波采集信号

其中超声波发射探头输出频率fc为40 kHz的方波激励信号，回波放大信号直接采用STM32单片机内部的12位A/D转换器，并设置为其最高采样频率6/7 MHz，即采样相位间隔δ为7π/75，则一个信号周期内的采样点数为21.43，采样序列总长度为2 000。图中的回波曲线不同于常规的回波曲线，原因在于系统的超声发射采用了正反向反复激励的输出方式[9,10]。

在超声回波移动正弦拟合过程中，首先选取21个采样点作为固定拟合长度进行移动正弦拟合(非整周期拟合)并计算出常数矩阵B与B-1，再根据采样信号获得3个初始正弦拟合参数，运用递推算法得到所有区间段的拟合参数值，最后根据式(13)得到移动正弦拟合过程中的第500～1 000数据点的拟合相位、拟合信号以及拟合误差的如图5所示.从图5 (a)可以看出采样信号的波谷对应的拟合相位为0，波峰对应的拟合相位为π (或-π)，从图5(b)可以看出:原始信号与拟合信号基本重合，图5 (c)显示出拟合误差较小且为白噪声，与预期结果一致。

图5 移动正弦拟合结果

快速移动正弦拟合求解后，对拟合求解参数ak，bk的绝对值和进行移动平均，利用包络曲线的求解公式(13)就能够快速计算出调幅信号的上包络。图6显示了有效部分采样点的原始信号及根据移动正弦拟合计算得到的包络曲线图。

图6 超声测距输出信号解包络

图6中为了更加清晰显示包络曲线的求解效果，包络曲线由包络幅值Ak与直流偏置ck之和替代，可以看出移动正弦拟合能够较为准确地获取调幅信号的包络曲线。

4 结 论

分析对比了直接正弦拟合与递推正弦拟合求解包络曲线的算法复杂度，表明该方法在高频采样中运算量与存储

空间量都具有明显的优势。提出利用拟合参数绝对值均值的移动平均方法代替幅值公式求解调幅信号的包络线，进一步降低了计算量，且适宜嵌入式系统的算法设计。同时该算法适用于非整周期采样拟合，应用于超声波信号的包络曲线解调实验中，能快速有效获取采样信号的包络曲线。