组合复信号的相位差法频率估计

2022-10-17罗久飞郑明轩

振动与冲击 2022年19期

罗久飞，郑明轩，李靖

(重庆邮电大学先进制造工程学院，重庆 400065)

精确而高效地测量谐波信号的幅值、频率和相位是信号处理中的一个经典课题。工程中通常采用离散频谱分析对测量信号进行频谱参数估计。然而，时域信号的截断和离散化处理，不可避免会造成能量泄露和栅栏效应，最终可能带来较大的估计误差[1]。针对这些问题，国内外学者提出了离散频谱校正技术，能够实现对频谱参数的精确估计。目前，离散频谱校正方法主要包括插值校正法[2-3]、能量重心校正法[4-5]、FFT+FT(fast Fourier transform+Fourier tronsform)连续细化校正法[6-7]和相位差校正法[8-10]。其中，相位差校正法具有算法简单，适用于各类窗函数，并且校正公式不受窗函数的限制，抗噪性能强等优点，被广泛应用于电力监测、振动分析、故障诊断、科氏力质量流量计等工程领域[11-16]。

McMahon等[17]首次提出利用相位插值精确估计随机噪声干扰下复指数信号的频率。丁康等针对两段连续的时域信号，研究了基于时移的相位差校正法的一种特殊情况。随后，Zhu等[18]针对这一特例进行了误差分析，推导出高斯白噪声下的归一化频率估计方差为

(1)

式中：δ为理论归一化频率校正量；ENBW为等效噪声带宽；SL(·)表示扇形损失[19]。丁康等[20]通过进一步研究，提出了可选择时移长度的时移相位差法。结果表明，算法的校正精度与时移长度L有关，L越大，精度越高。黄云志等[21]提出了基于窗中心平移的相位差校正法，相较于其他方法，该算法在随机噪声下具有更强的稳定性。丁康等[22]提出了基于时移和改变窗长的综合相位差校正法。杨志坚等[23]分析了高斯白噪声下，时移相位差校正法的频率估计精度，推导出矩形窗下的归一化频率均方根误差

(2)

以及汉宁窗下的归一化频率均方根误差

(3)

式中，C可由窗函数，时移长度，归一化频率校正量计算。Huang等[24]提出基于双子段相位差的频率估计方法，该方法类似于基于时移的相位差校正法。随后，Luo等[25]研究了非对称窗的相位特性，并提出一种基于非对称窗的相位差校正法，只需利用一段信号即可精确估计频率值。

基于时移的相位差校正法具有诸多优点。例如，在随机噪声的干扰下，可通过增大时移长度来提高算法的精度。此外，时移相位差法也适用于小阻尼衰减正弦信号的频率估计。然而，一旦时移长度过大，传统的时移相位差法就会出现相位绕卷问题，导致在频率偏移量δ较大时产生相当大的估计误差。同时，在强噪声干扰下，极易出现谱线定位错误的情况，这也可能导致相位绕卷的发生，从而影响最终的估计精度。因此，本文分析了时移相位差法中可能导致相位绕卷的主要因素，及其对算法的影响。提出了一种基于组合复信号相位差的频率估计算法，较大程度上避免了相位绕卷的发生，并且在长时移下依然具有良好的抗噪声性能。

1 基于时移相位差的频率估计方法

1.1 基于时移相位差的频率估计原理

不失一般性，考虑某单频正弦信号的N点时域采样序列

x0(n)=

(4)

式中：A、f0、θ0分别为信号的幅值、频率、初始相位；fs、N分别表示采样频率和采样点数。将x0(n)时移M点后表示为

x1(n)=

(5)

对x0(n)和x1(n)加长度为N的时移窗wN(n)，并分别进行离散傅里叶变换(discrete Fourier transform, DFT)得

(6)

(7)

式中，λ0为归一化频率，W(·)和τ(·)分别为对称窗w(n)的窗谱函数和相频响应函数。归一化频率λ0，采样频率fs，采样点数N之间的关系如下

(8)

式中：l为X0(k)中幅值最大处对应的谱线号；δ为归一化频率校正量。将k=l代入式(6)和式(7)，则有

X0(l)≈

(9)

X1(l)≈

(10)

根据式(9)和式(10)得两段信号幅值谱最高谱线处的相位分别为

φ0(l)=θ0+τ(-δ)

(11)

φ1(l)=θ1+τ(-δ)

(12)

由于通过DFT得到的相位为原始相位的主值，其取值范围为(-π,π)，则根据式(11)和式(12)可将相位差表示为

Δφ=θ1-θ0+2πh=φ1(l)-φ0(l)+2πh

(13)

Δφ=2πη+2πh

(14)

η的取值区间为(-1,1)。同时，相位差还可由时移系数S=M/N表示成

Δφ=2πλ0S=2π(l+δ)S

(15)

联立式(14)和式(15)可得

δS=η+h-lS

(16)

当S<1时，-0.5<δS<0.5。定义{·}为对1的取模运算，则有

{δS}={η+h-lS}={η-lS}

(17)

令γ={η-lS}，并进行如下调整

(18)

则可计算得到归一化频率偏移量

(19)

最后，可计算出正弦信号的频率为

f0=(l+δ)fs/N

(20)

1.2 相位绕卷问题对时移相位差法频率估计的影响

上述理论推导表明，在忽略相位绕卷问题的前提下，归一化频率校正量δ可由相位差Δφ、峰值谱线对应的频率l、时移系数S计算得到。其中，S为已知参数，Δφ和l可通过离散频谱分析获得。然而，在实际工程中，信号往往受到随机噪声的干扰，这极易造成相位的改变，同时还可能出现谱线定位错误的问题，从而导致相位绕卷，带来较大的估计误差。通过适当增大时移长度可以有效降低随机噪声导致的误差，但时移长度过大时，也可能会引入相位绕卷问题。

1.2.1 随机噪声的影响

由于δ在-0.5～0.5范围内变化，且时移系数S<1，则式(18)中γ在-0.5～0.5范围内变化。当γ接近于±0.5时，随机噪声的干扰会造成谱峰处的相位发生改变，进而导致γ的值可能超出-0.5～0.5范围，并且经过式(18)的调整后使δ的误差接近于1。文献[23]对这一误差产生的原因进行了深入分析，建议适当控制时移长度，使其不接近于采样点数，可减小此类误差。然而，时移长度的限制大大阻碍了传统算法抗噪声性能的进一步提升。

1.2.2 谱线定位错误的影响

在随机噪声的干扰下，可能出现谱线定位错误的情况，导致抽选到与谱峰相邻的左右两侧谱线。

当δ>0时，若谱线定位错误，则上述推导中式(17)应改写成

(21)

设ξ={η-(l+1)S}={γ-S}，若γ-S>-0.5，则ξ=γ-S，归一化频率校正量为

(22)

归一化频率为

(23)

因此该情况下的谱线定位错误对最终频率估计结果不产生影响。若γ-S<-0.5，则ξ=γ-S+1，归一化频率校正量为

(24)

归一化频率为

(25)

因此该情况下的谱线定位错误导致归一化频率的估计值产生1/S的误差。

同理，当δ<0时，若谱线定位错误，则式(17)应改写成

(26)

设ξ={η-(l-1)S}={γ+S}，若γ+S<0.5，则ξ=γ+S，归一化频率校正量为

(27)

归一化频率为

(28)

该情况下的谱线定位错误对最终频率估计结果不产生影响。若γ+S>0.5，则ξ=γ+S-1，归一化频率校正量为

(29)

归一化频率为

(30)

因此该情况下的谱线定位错误导致归一化频率的估计值产生-1/S的误差。

1.2.3 长时移的影响

若时移长度过大，也可能导致相位绕卷问题的发生。一旦时移系数超过1，δS的取值将超出限制区间(-0.5, 0.5)，此时需考虑相位绕卷问题，则式(19)改写为

(31)

式中，h为相位绕卷数。然而，相位绕卷数无法根据时移系数来确定，导致归一化频率校正量无法被唯一确定。例如，当时移系数S=3时，-1.5<δS<1.5。若归一化频率理论值λ0=32.35，则γ=0.05。根据式(31)，h可能的取值为-1,0,1。则此时的归一化频率校正量为

(32)

上述分析表明，随机噪声的干扰，谱线定位错误，以及长时移均可能导致相位绕卷，使传统时移相位差算法产生较大的估计误差。因此，需要研究一种能克服相位绕卷问题，并且在长时移下可以有效估计信号频率的算法。

2 组合复信号相位差法频率估计

2.1 组合复信号相位差法频率估计原理

利用两段时间序列构造新序列y0(n)为

y0(n)=x0(n)+jx1(n)

(33)

式中：x0(n)和x1(n)由式(4)和式(5)定义。对y0(n)加窗并进行快速傅里叶变换(fast Fourier transform, FFT)，则有

Y0(k)=X0(k)+jX1(k)

(34)

式中：X0(·)和X1(·)分别为x0(n)和x1(n)的DFT。当时移系数S>1时，为避免出现相位绕卷问题，需要增加一个初始归一化频率校正值λ1，以确保剩余归一化频率校正量δ1满足Sδ1在区间(-0.5,0.5)内。

设m为峰值谱线对应的频率值，并定义

(35)

(36)

对于矩形窗，归一化频率校正量为[27]

(37)

式中的符号由次高谱线位置决定，当次高谱线位置为m+1时，取正号，反之取负号。则计算得初始归一化频率校正值

λ1=|m+δ0|

(38)

将k=±λ1代入式(34)，可得

Y0(±λ1)=X0(±λ1)+jX1(±λ1)

(39)

其中Y0(·)可由下式计算

(40)

设X0(±λ1)=a0±jb0，X1(±λ1)=a1±jb1，则可将式(39)改写成

Y0(λ1)=X0(λ1)+jX1(λ1)=a0-b1+

j(a1+b0)

(41)

Y0(-λ1)=X0(-λ1)+jX1(-λ1)=a0+b1+

j(a1-b0)

(42)

根据式(41)和式(42)可知

(43)

φ0(λ1)=arctan(b0/a0)

(44)

φ1(λ1)=arctan(b1/a1)

(45)

根据式(44)和式(45)可知相位差为

Δφ=φ1(λ1)-φ0(λ1)

(46)

将Δφ代入式(13)～式(19)所示的时移相位差算法，计算得剩余归一化频率校正量δ1。则最终的频率校正值为

f1=(λ1+δ1)fs/N

(47)

2.2 强噪声下的算法改进

为进一步提高算法在强噪声干扰下的校正精度，对组合信号进行时域调整。将k=±m代入式(34)，并根据式(39)～式(46)相似推导计算得相位差Δφ1。利用Δφ1对x1(n)进行调整，并定义新序列y1(n)为

y1(n)=x0(n)+x1(n)e-jΔφ1

(48)

设Y1(·)为y1(·)的FFT结果，p为Y1(·)中幅值最大谱线对应的频率，根据式(35)～式(38)相似推导可得归一化频率校正值λ2。进一步地，利用迭代插值算法[28]将λ2进行1次迭代运算，得到更为精确的归一化频率校正值λ3。计算Y1(±λ3)，则有

Y1(λ3)=X0(λ3)+X1(λ3)e-jΔφ1

(49)

Y1(-λ3)=X0(-λ3)+X1(-λ3)e-jΔφ1

(50)

设X0(±λ3)=a2±jb2，X1(±λ3)=a3±jb3，并根据式(41)～式(43)的相似推导可解算得a2，b2，a3，b3，则相位差Δφ2为

Δφ2=arctan(b3/a3)-arctan(b2/a2)

(51)

将Δφ2代入式(13)～式(19)所示的时移相位差算法，计算归一化频率校正量δ2，则最终的频率校正值为

f2=(λ3+δ2)fs/N

(52)

2.3 算法复杂度

若只考虑必要的乘法和加法，并忽略所有与N无关的计算项，则本文所提组合复信号相位差法共需1次N点FFT计算，2次单点DFT计算，改进算法额外需要1次N点FFT计算和2次单点DFT计算。

对于组合复信号相位差算法，1次N点FFT计算需要(N/2)log2N次复数乘法和Nlog2N次复数加法[29]，1次单点DFT计算需要N次复数乘法和N-1次复数加法，则算法复杂度如表1所示。

表1 算法复杂度Tab.1 Calculation complexity

3 仿真分析

3.1 无噪声情况下组合复信号相位差算法的频率估计性能

本小节分析比较了无噪声情况下传统的时移相位差算法和组合复信号相位差算法的频率估计精度。根据式(4)和(5)生成测试信号，并设置采样点数N=256，采样率fs=256 Hz，幅值A=1，信号频率f0在区间[63.5,64.5]内以步长0.025变化，初相θ0在[-π,π]内以步长π/20变化。同时，为分析时移系数对算法的影响，设置时移系数以0.062 5的间隔在[0.062 5,4]的范围内变化。对于不同的相位，选择频率估计的最大绝对误差作为校正误差。矩形窗下时移相位差校正法和组合复信号的相位差校正法的最大绝对误差变化曲线，如图1所示。

(a) 基于时移的相位差法

(b) 组合复信号的相位差法图1 矩形窗下频率估计的最大绝对误差Fig.1 Maximum absolute error of frequency with rectangular window

由图1可知，在无噪声情况下，基于时移的相位差法的有效校正范围随着时移系数的增加而减小，并且由相位绕卷问题引入的最大误差约为1。本文提出的组合复信号的相位差法具有更高的频率估计精度，最大绝对误差随着时移系数的增大而逐渐减小，且不同时移系数下的最大绝对误差均低于2.3×10-3。

3.2 随机噪声干扰下组合复信号相位差法的频率估计性能

3.2.1 时移系数对估计精度的影响

为探究时移系数对组合复信号的相位差法估计精度的影响，设置采样点数N=256，采样率fs=256 Hz，幅值A=1，信号频率f0在区间[63.5,64.5]内随机分布，初相θ0在[-π,π]内随机分布，时移系数在[1,20]内以步距1变化。不同时移系数下组合复信号相位差法频率估计的RMSE变化曲线，如图2所示。

(a) 组合复信号相位差法

(b) 改进的组合复信号相位差法图2 不同时移系数S下的RMSEFig.2 RMSE versus S

仿真结果显示，对于不同信噪比，当时移系数较小时，频率估计的RMSE值随着时移系数的增大而减小。而在时移系数达到某一阈值后，RMSE随时移系数的增大而逐渐增大，并且信噪比越大，该时移系数阈值也越大。在信噪比较低的情况下，改进的组合复信号相位差法具有更强的抗噪声性能。

3.2.2 信噪比对估计精度的影响

本小节分析了不同N·SNR下，组合复信号的相位差法及其改进算法的频率估计性能。设置采样点数N=256，采样率fs=256 Hz，幅值A=1，信号频率f0在区间[63.5,64.5]内随机分布，初相θ0在[-π,π]内随机分布，时移系数S=3，结果如图3所示。

图3 不同N·SNR下频率估计的RMSEFig.3 RMSE versus N·SNR

由图3可知，在N·SNR大于510时，组合复信号的相位差法及其迭代改进算法的估计性能几乎相同，但是N·SNR小于510时，后者RMSE值明显小于前者。结果表明，在噪声较强的情况下，建议采用迭代运算来获得更高的频率估计精度。

3.2.3 不同频偏下各算法的频率估计性能比较

为分析频率偏移量对各算法估计精度的影响，设置采样点数N=256，采样率fs=256 Hz，幅值A=1，信号频率f0在区间[63.5,64.5]内以步长0.025变化，初相θ0在[-π,π]内随机，时移系数S=3。信噪比分别为10 dB，5 dB，0，-5 dB时各算法的频率估计均方根误差，如图4所示。

(a) SNR=10 dB

(b) SNR=5 dB

(d) SNR=-5 dB图4 不同频偏下各算法的RMSEFig.4 RMSE versus deviations

从图4可知，在较长时移下由于存在相位绕卷问题，导致基于时移的相位差法在频率偏移量|δ|>0.1时具有较大的均方根误差。而基于窗中心平移的相位差法由于进行了补零处理，相较于时移相位差法，具有更高的估计精度，但仍然可能发生相位绕卷，导致某些特定频率偏移处具有较大的估计误差。相比之下，本文提出的组合复信号的相位差法具有良好的频率估计性能及算法稳定性，其均方根误差与频率偏移量近似无关。