High-precision frequency estimation algorithm based on improved Chrip-Z transformation

单一频率复正弦信号可以表示为

其中，$ A $是正弦信号的振幅，$ {f_c} $是正弦信号的频率，$ \varphi $是正弦信号的相位。

对于信号以$ {f_s} $的频率进行采样，采样点数为$ N $。采样后的信号表示为

其中，$ {m_{\rm{c}}} = \dfrac{{{f_c} \cdot N}}{{{f_s}}} $。

对于信号$ s(n) $进行FFT运算，信号的频谱为

对于获得的频谱的幅值进行搜索，幅值最大处即为频率的估计值。

其中，$ {m_0} = {\mathop{\rm int}} \left[\dfrac{{{f_c} \cdot N}}{{{f_s}}}\right] $，$ {\mathop{\rm int}} [x] $表示x最接近的正数，$ {m_0} $是最大谱线的索引值，$ \Delta {f_0} = \dfrac{{{f_s}}}{N} $是频率分辨率。

由于存在栅栏效应，$ {m_c} $与$ {m_0} $之间存在一定偏差$ \delta $，这限制了频率估计的精度。

其中，$ \delta \in [ - 0.5,0.5] $。

为了提高频率估计精度，CZT算法被提出，CZT算法对于部分频率区间进行了细化，CZT变换表示为

其中，$ z_k=A W^{-k} $，$ A=A_0 e^{j \theta_0} $，$ W=W_0 e^{-j \phi_0} $，$ \theta_0= 2 \pi \dfrac{m_0-q}{N} $是起始采样角度，$ q $是整数用于控制细化频率区间的大小，$ \phi_0=2 \pi \dfrac{2 q}{M N} $是两相邻采样点之间的角度，$ M $是频率细化倍数。

当频率的细化倍数为$ M $时，CZT算法的频率分辨率$ \Delta {f_1} $比FFT算法的频率分辨率$ \Delta {f_0} $高了$ M $倍。

其中，$ q $用于控制细化区间的大小。

使用CZT算法对于频率进行估计，虽然能够得到较高精度的频率估计值，但是栅栏效应也依旧存在，这使得估计的频率$ \widehat {{f_{CZT}}} $与信号频率$ {f_c} $之间依旧存在一定的偏差。

其中，$ \delta \in [ - 0.5,0.5] $。

针对现有算法存在的问题，需要研究一种新的频率估计算法，能够解决栅栏效应的问题从而提升频率估计精度，以满足实际工程应用的要求。

为解决栅栏效应对于频率估计精度的影响，提出在CZT算法的基础上，利用细化频谱最大谱线及其左右谱线对频率偏移值进行估计，从而提升频率估计的精度。所提算法将频率估计分成粗估计和细估计两个过程，下文将介绍算法的详细原理和整体流程，并分析了所提算法的计算量。

首先，使用FFT算法获得信号$ s(n) $的频谱，并计算出频谱的幅值$ \left| {S(k)} \right| $。

对于频谱幅值进行搜索，幅值最大处即为频率的粗估计值$ \widehat {{f_{CZT}}} $和该点的频率索引值$ {m_0} $。根据频率的粗估计值$ \widehat {{f_{CZT}}} $，以此为中心设置频率细化区间$ ({f_1},{f_2}) $。

在确定的频率细化区间对于信号做CZT变换。在进行CZT变换时，设置$ {A_0} = 1 $、$ {W_0} = 1 $，求得原信号$ s(n) $的细化频谱$ {S_{CZT}}(k) $：

然后，计算细化频谱的幅值$ {S_{CZT}}(k) $

对于$ {S_{CZT}}(k) $进行搜索，幅值最大处即为频率的估计值$ \widehat {{f_{CZT}}} $。

其中，$ {m_1} $是细化后频谱最大谱线的索引值。

信号的真实频率与细化后频谱的采样点之间存在一定的偏差$ \delta $，信号频率$ {f_c} $表示为

使用细化频谱中最大谱线及其左右谱线的幅值求出误差$ \delta $，对于频率估计值$ \widehat {{f_{CZT}}} $进行修正。最大谱线的幅值为

最大谱线左侧谱线的幅值为

最大谱线右侧谱线的幅值为

$ \left|S_{C Z T}\left(m_1+1\right)\right| $和$ \left|S_{C Z T}\left(m_1\right)\right| $的比值记为$ {a_1} $, $ {a_1} $表示为

对于较大的$ N $，$ {a_1} $可以表示为

$ \left|S_{C Z T}\left(m_1-1\right)\right| $和$ \left|S_{C Z T}\left(m_1\right)\right| $的比值记为$ {a_2} $, $ {a_2} $表示为

对于较大的$ N $，$ {a_2} $可以表示为

将式(20)和式(22)相减，利用三角函数的和差化积公式进行化简求解，可得$ \delta $的表达式为

将$ \delta $的估计值代入式(15)即可获得频率的估计值$ \widehat {{f_{c}}} $。

算法具体流程如下：

步骤1：对于差频信号进行FFT变换，获得最大谱线频率索引$ {m_0} $。

步骤2：根据$ {m_0} $确定频率细化区间，对于信号进行CZT变换。

步骤3：对于CZT变换后的频谱进行搜索，获得最大谱线的幅值$ \left|S_{C Z T}\left(m_1\right)\right| $及其左右两根谱线的幅值$ \left|S_{C Z T}\left(m_1-1\right)\right| $、$ \left|S_{C Z T}\left(m_1+1\right)\right| $。

步骤4：计算$ {a_1} = \dfrac{{\left| {{S_{CZT}}\left( {{m_1} + 1} \right)} \right|}}{{\left| {{S_{CZT}}\left( {{m_1}} \right)} \right|}} $。

步骤5：计算$ {a_2} = \dfrac{{\left| {{S_{CZT}}\left( {{m_1} - 1} \right)} \right|}}{{\left| {{S_{CZT}}\left( {{m_1}} \right)} \right|}} $。

步骤:6：计算$ \delta = \dfrac{{{a_1} - {a_2}}}{{{a_1} + {a_2} - 2\cos \left( {\dfrac{{2q\pi }}{M}} \right)}} $。

步骤:7：计算$ {\hat f_c} = \left( {{m_0} - q} \right)\Delta {f_0} + \left( {{m_1} + \delta } \right)\Delta {f_1} $。

对于工程应用而言，算法的计算量是一个重要的参考因素。复数乘法的计算量远大于复数加法，因此在讨论时忽略了复数加法的计算量，下表中列出了本文比较算法的计算量。

从表1中可以看出，FFT算法和Rife算法的计算量是最小的，补零FFT算法的计算量是最大的。Zoom-FFT算法的计算量是FFT算法的两倍。CZT算法的计算量与细化倍数相关，在细化倍数相同的情况下，本文所提算法的计算量与CZT算法基本一致，略大于FFT算法。实际上，本文所提算法可以以较低的细化倍数达到较高的频率估计精度，这能够减少算法需要的计算量。

为了验证算法的性能，使用MATLAB软件对于算法性能进行仿真，将所提算法与FFT算法、补零FFT算法、CZT算法、Rife算法、Zoom-FFT算法和克拉美罗下限进行比较。FFT补零法补了$ 3N $长度的零。Zoom-FFT算法的细化倍数为64。CZT算法和所提算法设置了$ q = 1 $和$ q = 2 $两种细化区间。当q=1时，CZT算法的细化倍数设置为64倍和32倍，本文所提算法的细化倍数设置为32倍。当$ q = 2 $时，CZT算法和本文所提算法的细化倍数均设置32倍。单一频率复正弦信号的克拉美罗下限^[21]为

使用RMSE来衡量频率估计的精度，使用SD来衡量算法的稳定性。

其中，$ L $为仿真次数。

在仿真时，向信号$ s(n) $加入高斯白噪声，并进行采样后得：

其中，$ w(n) $表示的是高斯白噪声。仿真中，采样频率为$ {f_s}=92.7835e3Hz $，初相位φ随机取值，振幅$ A=1 $，进行10000次的独立蒙特卡洛仿真实验。

在信噪比为0 dB时，设置信号$ s(n) $的频率$ {f_c}=5100Hz $，比较不同采样点数下各算法实际的运行时长，结果如表2所示。

从表2中可以看出，FFT算法和Rife算法的运行时长是最短的，补零FFT算法的运算时间是最长的。Zoom-FFT算法在采样点数小时，运行时长小于CZT算法，当采样点数大时运行时长大于CZT算法。CZT算法的计算时间和细化倍数相关，细化倍数越高，运行时间越长，CZT算法的运行时长略大于FFT算法。本文提出的改进算法与CZT算法的运行时间基本一致，不会影响计算的速度。

设置信号$ s(n) $的频率$ {f_c}=5100Hz $，采样点数$ N = 1024 $，向信号中加入高斯白噪声，比较信噪比在$ [-12dB,12dB] $的情况下各算法的性能，结果如图3和图4所示。

从图3和图4中可以看出，由于存在栅栏效应，FFT的频率估计值一直维持在一个稳定值，这使得FFT算法的频率估计误差为一个定值，SD也恒为零。只有提升频率分辨率才能提升FFT算法的精度，这需要消耗大量的计算资源。补零FFT算法的频率估计精度高于FFT算法，与FFT算法类似，频率估计的精度受到了栅栏效应的影响，只有提升频率分辨率，才能提升频率估计的精度。在同样的细化倍数下，Zoom-FFT算法的精度低于CZT算法。CZT算法的精度受细化倍数和细化区间的影响，当细化倍数一定时，细化区间小的算法频率估计精度高；当细化区间一定时，细化倍数越高频率估计精度越高。虽然CZT算法具有较FFT算法、Zoom-FFT算法具有较高的精度，但是CZT算法依旧受到了栅栏效应的影响，信噪比高时，频率估计误差将为一个稳定值。CZT算法的SD一直都较小，频率估计具有较高的稳定性。Rife算法在信噪比大时，具有较好的性能，但是受噪声影响较大，SD受噪声影响较大。本文所提算法的SD一直较小，频率估计的结果稳定性较强，受噪声影响较小。本文所提算法性能与CRLB接近，频率估计的精度高于上述所有算法，在信噪比为$ [-12dB,12dB] $的情况下均接近CRLB，具有较好的抗噪性能。本文所提算法在细化倍数为32时，频率估计精度高于细化倍数为64的CZT算法，能以更低的细化倍数获得更高的频率估计精度。

为了分析本文所提算法在不同频率下的性能，在信噪比为0 dB时，采样点数$ N = 1024 $，仿真信号频率$ {f_c} $由5000 Hz以1 Hz的步长递增至5150 Hz的条件下，进行仿真比较，结果如图5和图6所示。

从图5和图6中可以看出，FFT算法、补零FFT算法、CZT算法，这几个受栅栏效应影响的算法，在信号频率改变时，频率估计RMSE呈周期性，当信号频率为频谱两个采样点中间位置时，受噪声影响频率估计值会出现抖动，导致SD较大。Rife算法受到算法本身性能的限制，RMSE在部分频段时较大，SD也一直大于本文所提算法。本文所提算法的RMSE一直趋近于CRLB优于所有比较的算法，频率估计精度不受频率变化的影响，且SD一直较小，具有较好的稳定性。

为了验证算法在实际工程应用中的性能，使用FMCW雷达来验证算法性能。FMCW雷达测量的距离与中频信号的频率成正比

实验中，雷达的采样率为$ {f_s}=123.2e3Hz $，采样点数$ N = 2048 $，带宽$ B = 999.807e6Hz $。实验场景如图7所示，使用导轨控制雷达与墙面之间的距离，使用激光测距仪来确定雷达与墙面之间的距离值。设置雷达与墙面之间的距离为1049 mm到1549 mm，每间隔50 mm测量一次距离，对于FMCW雷达采集到的数据进行分析处理，利用式(28)求出距离估计值。

通过第三部分的仿真结果可知，补零FFT算法的性能时高于FFT算法的。因此在实测数据分析时，未使用FFT算法进行比较。补零FFT算法补了$ 3N $长度的零。Zoom-FFT算法的细化倍数为32。CZT算法和所提算法的细化区间设置为$ q = 1 $，算法的细化倍数为32。

从表3中可以看出，本文所提算法的测距性能时优于补零FT算法和常规CZT算法的。Rife算法少数点的测距精度高于本文所提算法，这是因为Rife算法的精度与误差相关，在部分区间估计精度高，算法的鲁棒性较差。Zoom-FFT算法的精度与常规CZT算法相近，仅在个别点的性能优于本文所提算法，整体性能不如本文所提算法。

为解决栅栏效应对于频率估计精度影响，提升频率估计的精度和算法抗噪性能，提出了一种基于改进Chrip-Z变换的高精度频率估计算法，该算法将频率估计过程分成粗估计和细估计两个过程，通过FFT算法确定频率细化区间，使用CZT变换对于频谱进行细化，使用细化频谱最大谱线及其左右谱线的幅值对于频率估计值进行修正，获得高精度的频率估计值。

在不同条件下对于算法进行仿真，表明本文所提算法具有较好的抗噪声性能，算法具有较高的频率估计精度，RMSE与CRLB贴近，优于现有算法。本文所提算法不受信号频率的影响，在各频率下都具有较好的性能，优于现有算法。本文所提算法能以较低的细化倍数，获得较高精度的频率估计值，降低了计算量，能够满足工程应用对于算法实时性的要求。通过实验对于本文算法进行验证，将算法应用于FMCW雷达测距场景中，实验表明本文所提算法在实际应用场景中有较好的表现，优于现有算法。