CN103941087B - 欠采样速率下的高频余弦信号的频率测量方法及其装置 - Google Patents

Abstract

本发明公开了一种欠采样速率下的高频余弦信号的频率测量方法及其装置，涉及数字信号处理领域，本发明结合过零点类型和apFFT相位谱分布特征，从2L个峰值幅度谱位置中挑选出L个谱位置索引值；根据余数筛选结果，用全相位比值内插法对各路apFFT谱分别进行谱校正，得到每路的频率估计；利用各路谱校正得到的频率估计值，作为余数，再按照闭合解析形式的中国余数定理对这些余数处理，通过频率估计值重构出原始高频信号的频率。本发明实现了多路低速率欠采样下的高频余弦信号的频率测量，大大提高了对高频信号频率测量的范围；采用新型频谱校正方法提高了高频余弦信号的频率测量精度，保证了新型低速率欠采样下的高频余弦信号的抗噪性能。

Description

欠采样速率下的高频余弦信号的频率测量方法及其装置

技术领域

本发明涉及数字信号处理领域，尤其涉及一种欠采样速率下的高频余弦信号的频率测量方法及其装置。

背景技术

高频信号的频率估计与检测是雷达通信、声呐、地震监测、故障诊断乃至医学医疗等领域中的至关重要的问题，引起国内外相关学者的普遍关注。随着信号频率升高，密集谱识别工程的实现难度会进一步增大。而在无线电等工程领域，却越来越频繁涉及频率值处于百兆(MHz)、千兆(GHz)、万兆(10GHz)数量级的信号采样、分析与处理问题。如在IEEE802.15.3a、IEEE802.15.4a等标准方案中，物理层信号采用了超宽带（UWB）技术，为实现室内精确定位及高速信息传输，要求其脉冲信号频率最高值达到3.1GHz以上[1]。毫无疑问，高频信号处理首先要求进行高速采样，这样在工程上必然对模数转化设备及其后期的FPGA等数字处理器件的性能及成本提出高的要求。

因而，在无法满足密集均匀采样的情况下，如何实现高分辨率地识别并分离出密集谱成成分是学术界和工程界迫切需要解决的问题。仅靠改进硬件设备的数据采集性能，其作用是非常有限的（如提高A/D采样速率就必须付出更高的功耗与硬件成本作为代价），只有在信号处理领域提出新的谱分析理论方法，才能根本上解决这类问题。

为解决低速欠采样下(令信号频率为f₀，要求采样速率f_s<<2f₀)的高频信号的频率估计问题，我国古老的中国余数定理[2][3](Chinese Remainder Theorem,CRT)被引入该领域中。中国余数定理研究的是这样一个问题：为重构某一未知整数N，给定的只有一组相互之间满足互素关系的整数模值：M₁,M₂,...,M_L，及其未知整数N对各模值M_i的模除结果r_i(即余数r_i,满足r_i=N mod M_i),i=1,...,L，从这些个余数r_i重构未知整数N的问题。CRT具有许多应用，如密码学[3]，信道编码[4][5]，信号处理[2][6][7][8]，以及雷达系统[9]-[20]等。

近年来，各种CRT重构算法出现很多新的成果。如文献[4]-[6]提出了余数数目冗余方法(Remainder Number Redundancy method)，该方法可以解决给定L个余数重构ρ个未知整数的问题(要求L>>ρ)；文献[7][10][18][20]提出了余数冗余方法(RemainderRedundancy method)。其中，对于余数冗余方法非常适合用于低速率采样下的信号频率估计。特别是文献[21]提出了一种闭合解析形式的CRT重构方法，该方法放宽了经典中国余数定理要求的L个模值互素的要求，仅仅要求L个模值具有某一最大公约数(greatest commondivider,gcd)为M即可，即满足M=gcd{M₁,M₂,...,M_L}，且文献[21]将CRT重构对象从整数域推广到整个实数领域，且该方法对余数误差具有很高的鲁棒性，因而具有很高的实用价值。

但是现有的用中国余数定理的信号频率估计，仅仅限于复指数信号的频率估计，如文献[7]将CRT用于欠采样下的复指数频率估计中，如文献[10][9]将CRT用于合成孔径雷达系统的解相(phase unwrapping)中。文献[7][9][10]中，所得到的L路余数都是从各路低速率样本的FFT谱峰得到。

但是目前还未见有文献研究欠采样下的余弦信号的频率估计问题。该问题所涉及的波形如图1所示(L=3路重构，高速采样频率f₀=127.2Hz，低速率采样频率f_s1=16Hz，f_s2=24Hz，f_s3=40Hz，M=8)

对于图1所示的高频余弦信号的频率估计问题，其测频难度在于：余弦信号包含两个复指数边带，故对于某一频率的余弦信号，从L路低速采样做谱峰搜索，必然得到2L路余数，怎样从2L余数中挑选出CRT重构所需的有效的L路余数是个非常棘手的问题。

发明内容

本发明提供了一种欠采样速率下的高频余弦信号的频率测量方法及其装置，本发明实现了在欠采样速率下对高频信号的估计，详见下文描述：

一种欠采样速率下的高频余弦信号的频率测量方法，所述方法包括以下步骤：

（1）对高频模拟余弦信号进行过零点检测，取任一过零点作为高频模拟余弦信号中心采样位置；

（2）以过零点为中心，分别以频率f_s1～f_sL对高频模拟余弦信号进行L路低速率采样，每路均采集2M-1个样点，并存储；

（3）对低速率采样得到的L路信号，分别进行加汉宁双窗全相位快速傅里叶变换；

（4）结合过零点类型和apFFT相位谱分布特征，从2L个峰值幅度谱位置中挑选出L个谱位置索引值；

（5）根据余数筛选结果，用全相位比值内插法对各路apFFT谱分别进行谱校正，得到每路的频率估计

（6）利用各路谱校正得到的频率估计值，作为余数，再按照闭合解析形式的中国余数定理对这些余数处理，通过频率估计值重构出原始高频信号的频率f₀。

所述结合过零点类型和apFFT相位谱分布特征，从2L个峰值幅度谱位置中挑选出L个谱位置索引值的步骤具体为：

(a)由过零点过渡情况，确定过零点瞬间相位(+90°或-90°)；

(b)从每路相位谱分布图中，确定与过零点瞬间相位一致的左半边带或右半边带；

(c)从每路确定的半边带中，找出对应的幅度谱峰值位置，作为余数筛选所需的位置索引{k_i ^*,i=1,...,L}。

所述根据余数筛选结果，用全相位比值内插法对各路apFFT谱分别进行谱校正，得到每路的频率估计的步骤具体为：

(a)对输入的2N-1样本做阶数N的apFFT谱分析，得到谱估计结果Y(k)，k=0,...,N-1；

(b)找出Y(k)的峰值谱位置k=k^*，在振幅谱线中选取相邻最大的两根进行比值，将该比值的平方根记为v，Y(k^*)为谱峰幅值，Y(k^*-1)和Y(k^*+1)分别表示峰值两侧的谱线的幅值，即

$v=\sqrt{\frac{\left|Y\left(k^{*}\right)\right|}{\max \left(\right|Y(k^{*}-1)|,|Y(k^{*}+1)\left|\right)}}---\left(2\right)$

(c)根据v求取比例偏差因子△k，其中，

$\mathrm{\&Delta;k}=\frac{2-v}{1+v}---\left(3\right)$

(d)根据比例偏差因子△k进行频率校正

若峰值谱线处于k=k^*的位置上，则校正后的值为

即通过上述步骤(a)至步骤(d)可以得到一路频率估计重复执行可以得到每路的频率估计

所述方法所能达到的频率测量范围为f_max=lcm(f_s1，f_s2，...，f_sL)。

一种欠采样速率下的高频余弦信号的频率测量装置，所述频率测量装置包括：触发电路、模数转化器、DSP器件、输出驱动及其显示电路，

待测信号经过所述触发电路的过零点检测处理来决定其初始相位，然后经过多路采样频率分别为f_s1，f_s2，...，f_sL的所述模数转化器采样得到样本序列{x₁(n)，x₂(n)，...，x_L(n)}，分别以并行数字输入的形式进入所述DSP器件，经过所述DSP器件的处理，得到高频信号的频率估计；最后通过所述输出驱动显示及其显示电路显示出频率值。

本发明提供的技术方案的有益效果是：

第一、实现了在低速率欠采样条件下对高频信号频率的测量，大大提高了对高频信号频率测量的范围。

对于传统的频率测量，对于每一单路的采样速率f_si，其测量范围仅为(0,f_si/2)。而本发明由于采用多路低速率欠采样方案对高频正弦信号频率联合进行频率测量，其测量范围大大增加。对于L路低速率采样频率f_s1，f_s2，...，f_sL，本发明所能精确测量的频率范围为f_max=lcm(f_s1，f_s2，...，f_sL)。

例如，实验1中各路低速率采样频率分别为：f_s1=2048Hz，f_s2=3072Hz，f_s3=5120Hz，各路的测量范围分别为(0,1024Hz]，(0,1536Hz]，(0,2560Hz]。对于本发明，则根据中国余数定理，最大可测频率为f_max=lcm(f_s1，f_s2，f_s3)=3.072×10⁴Hz，则测量范围为(0,3.072×10⁴Hz]，测量范围提高了1个数量级。

第二、对于高频测量，本发明采用多路低速率采样，资源耗费少，大大节省了硬件成本。

相对而言，如果采用以往的高速采样，就必须提高A/D采样速率，必然会为此付出更高的功耗与硬件成本作为代价。

例如，实验1中各路低速率采样频率分别为：f_s1=2048Hz，f_s2=3072Hz，f_s3=5120Hz，则测量范围为(0,3.072×10⁴Hz]。同样地，如果进行高速率采样，要达到同样的测量范围，要求的采样频率至少为f_s=6.144×10⁴Hz。因此本发明大大降低了采样速率，节省了硬件成本。

第三、在无噪情况下，本发明的高频测量方法对高频信号的测量几乎可以做到0误差的精确测量。

例如，实验2中对高频信号的测量结果统计，其相对测量误差σ₂约处于10^-8数量级，对应的绝对测量误差仅为约10^-3Hz数量级。这在高频测量中可以认为是0误差的精确测量。

第四、提高了高频测量的抗噪声性能，在噪声条件下也能够很准确的测量出高频信号频率。

例如，实验3中在噪声环境下对高频信号的频率测量结果统计，其相对测量误差σ_noise约处于10^-7数量级，对应的绝对测量误差不超过10^-2Hz数量级，仅仅比无噪情况下多一个数量级。因此，本测量装置具有很好的抗噪声性能。

附图说明

图1对高频余弦信号进行三路低频采样的示意图；

图2低采样速率下的高精度高频测量方法的流程图；

图3对高频信号的过零点检测示例图；

图4全相位FFT谱分析框图对照(M=4)的示意图；

图5举例说明每路采样信号的apFFT分析谱图；

图6本发明的硬件实施图；

图7 DSP内部程序流图；

图8不同信噪比下频率估计的均方根误差。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

101：对高频模拟余弦信号进行过零点检测，取任一过零点作为高频模拟余弦信号中心采样位置；

其中，对输入的高频模拟信号x(t)进行过零点检测，对于余弦信号而言，过零点存在如图3所示的两种情况：

图3中，

(a)图对应为从正波形到负波形过零，这时过零点的瞬间相位为π/2。

(b)图对应为从负波形到正波形过零，这时过零点的瞬间相位为-π/2。

过零点处π/2或-π/2的瞬间相位的符号，对于后面步骤的每路信号的apFFT比值校正法的两个余数中选取一个，起到决定作用。

模拟信号经过简单的触发电路可以很容易地确定过零点时刻。

102：以过零点为中心，分别以f_s1～f_sL对高频模拟余弦信号进行L路低速率采样，每路均采集2M-1个样点，并存储；

其中，以过零点为中心采集2M-1点，一共采集L路信号。对于高频模拟余弦信号x(t)=a·cos(2πf₀t+π/2)，所测频率为f₀，采样频率分别为f_s1～f_sL，则各路采样信号为（其中n=-M+1,...,M-1）：

$\begin{array}{cc}{x}_1\left(n\right)=a.\cos (2\mathrm{\&pi;}\frac{f_0}{f_{s1}}n+\frac{\mathrm{\&pi;}}{2})& \\ x_2\left(n\right)=a.\cos (2\mathrm{\&pi;}\frac{f_0}{f_{s2}}n+\frac{\mathrm{\&pi;}}{2})& \\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ x_L\left(n\right)=a.\cos (2\mathrm{\&pi;}\frac{f_0}{f_{\mathrm{sL}}}n+\frac{\mathrm{\&pi;}}{2})& \end{array}---\left(1\right)$

采样速率要求f_s1～f_sL满足具有公约数M，且除以公约数M后是两两互素的。

103：对低速率采样得到的L路信号，分别进行加汉宁双窗的全相位快速傅里叶变换(apFFT)，记为{Y_i(k),i=0,...,L-1；k=0,...,M-1,}；(其中i表示采样路序号，k表示谱线标号)；

图4中，即需用长为(2M-1)的卷积窗w_c(由两个长度为M的汉宁窗的卷积而成)对中心样点x(0)前后(2M-1)个数据进行加权，然后将间隔为M的数据两两进行重叠相加，再对重叠相加后的数据进行FFT即得全相位谱分析结果。

例如：直接对x(t)=a·cos(2πf₀t+π/2)进行各路低速欠采样，进行apFFT操作而取峰值处的频率值，所得到的是一组整数。涉及的波形如图5所示：

L=3路采样，高速采样频率f₀=601.1520Hz，低速率采样频率f_s1=128Hz，f_s2=192Hz，f_s3=320Hz，易推出M=gcd{128,192,320}=64)。

图5的apFFT幅度谱Y_i(k)和相位谱具有如下规律：

(a)因实信号缘故，每路幅度谱均有两个谱峰，其位置关于频率轴中心对称。

(b)每路信号两个谱峰位置对应的相位谱值是大小相等，正负符号相反。

(c)不同路数间的左右频率半轴的相位谱值的正负符号出现的顺序有差异。

(d)从幅度谱中很明显看出存在谱泄漏；因此需要对频谱的峰值位置进行校正。

以上图5展现的峰值谱和相位谱分布规律，可为后续CRT处理提供如下依据：

1)谱泄漏分布提供提升余数精度依据：

由于所测信号频率常常是任意的，很难保证图5的理想峰值谱恰好落在整数倍的谱线位置，而是常常分布在以峰值谱线为中心的几根谱线上(即形成谱泄漏)。可以对这些泄漏出来的谱线，做进一步插值处理，估计出理性谱位置。从而提高CRT所需的余数精度。

2)相位谱分布为余数索引筛选提供分类依据：

由于中国余数定理所需的L路余数只能从峰值谱位置去确定，而图5中，每路的apFFT谱存在两个谱峰，总共有2L个谱峰。故需要从中筛选出L个谱峰索引给CRT提供余数。而图5的相位谱分布规律，给余数索引筛选提供分类依据。

104：结合过零点类型和apFFT相位谱分布特征，从2L个峰值幅度谱位置中挑选出L个谱位置索引值；

具体步骤如下：

(d)由过零点过渡情况，确定过零点瞬间相位(+90°或-90°)；

(e)从每路相位谱分布图中，确定与过零点瞬间相位一致的左半边带或右半边带；

(f)从每路确定的半边带中，找出对应的幅度谱峰值位置，作为余数筛选所需的位置索

引{k_i ^*,i=1,...,L}。

例如：过零点是从正值到负值的过零点，其瞬间相位为+90°，从图5的相位谱图中可看出：第1路到第3路与+90o相一致的半边带分别为：右边带、左边带和右边带，从这些边带提取出的峰值谱位置{k_i ^*,i=1,...,L}分别为k₁ ^*=45,k₂ ^*=8,k₃ ^*=56。

105：根据余数筛选结果，用全相位比值内插法对各路apFFT谱分别进行谱校正，得到每路的频率估计

如前所述：由于各路真实欠采样信号的频率不一定恰恰落在整数倍的谱线位置上(该位置可以由上个步骤的余数索引筛选获得)，故若取该位置的频率值作为CRT的余数，必然会引入测量误差而降低测量精度。

所以本发明提出“全相位比值内插法”用以提高CRT余数精度，该方法需对峰值谱和泄漏出的旁谱进行校正，从而可将每路的频率估计精确到小数。

“全相位比值内插法”的处理步骤介绍：

(a)对输入的2N-1样本做阶数N的apFFT谱分析(选用的窗为汉宁双窗)，得到谱估计结果Y(k)，k=0,...,N-1；

(b)找出Y(k)的峰值谱位置k=k^*，在振幅谱线中选取相邻最大的两根进行比值(即最高谱幅值除以旁边的次高谱幅值)，将该比值的平方根记为v，Y(k^*)为谱峰幅值，Y(k^*-1)和Y(k^*+1)分别表示峰值两侧的谱线的幅值，即

$v=\sqrt{\frac{\left|Y\left(k^{*}\right)\right|}{\max \left(\right|Y(k^{*}-1)|,|Y(k^{*}+1)\left|\right)}}---\left(2\right)$

(c)根据v求取比例偏差因子△k，其中，

$\mathrm{\&Delta;k}=\frac{2-v}{1+v}---\left(3\right)$

(d)根据比例偏差因子△k进行频率校正。具体为

若峰值谱线处于k=k^*的位置上，则校正后的值为

即通过上述步骤(a)至步骤(d)可以得到一路频率估计重复执行可以得到每路的频率估计

106：利用各路谱校正得到的频率估计值，作为余数，再按照闭合解析形式的中国余数定理对这些余数处理，通过频率估计值重构出原始高频信号的频率f₀。

上述所得出的各路信号的频率值即为中国余数定理中所需的余数。将各路采样频率f_s1，f_s2，...，f_sL作为CRT的各路模值，结合最大公约数M值，按照如下闭合解析形式的中国余数定理的算法步骤，估计出高频信号频率值。

(a)从所给的余数（1≤i≤L）计算参数其中：

${\hat{q}}_{i,1}=\left[\frac{{\hat{f}}_i-{\hat{f}}_1}{M}\right],2\&le;i\&le;L---\left(5\right)$

(b)计算模除Γ_i的余数：

${\widehat{\mathrm{\&xi;}}}_{i,1}={\hat{q}}_{i,1}{\stackrel{\&OverBar;}{\mathrm{\&Gamma;}}}_{i,1}\mathrm{mod}{\mathrm{\&Gamma;}}_i---\left(6\right)$

对于2≤i≤L，其中，是Γ₁关于Γ_i的模逆，可以提前算出。

(c)计算

${\hat{n}}_1=\underset{i=2}{\overset{L}{\mathrm{\&Sigma;}}}{\widehat{\mathrm{\&xi;}}}_{i,1}{b}_{i,1}\frac{{\mathrm{\&gamma;}}_1}{{\mathrm{\&Gamma;}}_i}\mathrm{mod}{\mathrm{\&gamma;}}_1---\left(7\right)$

其中b_i,1是关于Γ_i的模逆，且γ₁由定义。

(d)计算:(2≤i≤L)

${\hat{n}}_i=\frac{{\hat{n}}_1{\mathrm{\&Gamma;}}_1-{\hat{q}}_{i,1}}{{\mathrm{\&Gamma;}}_i}---\left(8\right)$

(e)计算

由上述得到的对于1≤i≤L，有：

${\hat{f}}_{0i}={\hat{n}}_{i}M{\mathrm{\&Gamma;}}_i+{\hat{f}}_i---\left(9\right)$

为了减小误差，取平均值：

${\hat{f}}_0=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\hat{f}}_{0i}---\left(10\right)$

通过以上算法即可测出频率值

文献[21]指出，基于中国余数定理的最大可测频率为f_max=lcm(f_s1，f_s2，...，f_sL)，其中lcm为最小公倍数(least common multiplier)。

本发明所能达到的频率测量范围为f_max=lcm(f_s1，f_s2，...，f_sL)，上述步骤中，要求f_s1～f_sL为整数，且其最大公约数（gcd）为M，且除以公约数M后是两两互素的。

参见图6，该欠采样速率下的高频余弦信号的频率测量装置包括：触发电路、模数转化器、DSP（Digital Signal Processor，数字信号处理器）器件、输出驱动及其显示电路，待测信号首先经过触发电路的过零点检测处理来决定其初始相位，然后经过多路采样频率分别为f_s1，f_s2，...，f_sL的A/D（模数转化器）采样得到样本序列{x₁(n)，x₂(n)，...，x_L(n)}，分别以并行数字输入的形式进入DSP器件，经过DSP器件的内部算法处理，得到高频信号的频率估计；最后借助输出驱动显示及其显示电路显示出频率值。

其中，图6的DSP为核心器件，在信号频率估计过程中，完成如下主要功能：

(1)调用核心算法，完成对各路接收信号的频率估计与校正，以及待测高频信号的频率估计处理；

(2)根据实际需要调整采样率f_s1，f_s2，...，f_sL，尽量地满足实际需要；

(3)将频率估计结果实时输出至驱动和显示模块。

需指出，由于采用了数字化的估计方法，因而决定了图6系统的复杂度、实时程度和稳定度的主要因素并不是图6中DSP器件的外围连接，而是DSP内部程序存储器所存储的核心估计算法。

DSP器件的内部程序流程如图7所示。

本发明将所提出的“欠采样速率下的高频余弦信号的高精度测量方法”这一核心估计算法植入DSP器件内，基于此完成高精度、低复杂度、高效的高频余弦信号的频率估计。

图7流程分为如下几个步骤：

(1)首先需根据具体应用要求（如医学和军事等的具体测量要求），粗略估计高频信号的频率范围，并根据具体需要设定测量范围和各路采样频率f_s1，f_s2，...，f_sL。该步骤是从工程方面提出具体需求，以使得后续流程有针对性地进行处理。

(2)然后，CPU主控器从I/O端口读采样数据，进入内部RAM。

(3)后续的“去直流处理”，是为了消除待测信号中的直流成分的影响。否则，直流成分的存在，会降低测量精度。直流成分很容易测出，仅需计算样点的平均值即可得到。

(4)按图2本发明的处理过程进行频率测量是DSP算法最核心的部分，运行该算法后，即可得到频率测量值。

(5)判断本发明方法是否满足工程需求，若不满足，程序返回，重新根据要求设定采样频率和最大可测范围。

(6)直至测量结果符合工程要求，然后通过DSP的输出总线输出至外部显示驱动设备，将频率测量结果进行数码显示。

需指出，由于采用了DSP实现，使得整个频率估计操作变得更为灵活，可根据信号所包含的各种分量的具体情况，通过编程灵活改变算法的内部参数设置，如谱分析的阶数M、采样频率f_s1，f_s2，...，f_sL等。

实验结果

将实验得到的频率估计的精度用相对误差来反映，定义如下

$\mathrm{\&sigma;}=\frac{|{\hat{f}}_0-f_0|}{f_0}\&times;100\%---\left(11\right)$

(1)无噪情况下不引入频谱校正的情况

实验1

选取的余弦信号为x(t)=a·cos(2πf₀n+π/2)，其中a=2。令最大公约数M=1024，f₀为待测高频信号的频率。采用L=3路低速率采样，各路低速率采样频率分别为：f_s1=2048Hz，f_s2=3072Hz，f_s3=5120Hz，则根据中国余数定理，最大可测频率为f_max=lcm(f_s1，f_s2，f_s3)=3.072×10⁴Hz。

在无噪声的情况下，在(0，f_max]范围内任意选取10个频率值f₀，利用本发明的装置进行频率测量，其中为直接经过apFFT处理直接从幅度谱峰位置并经过余数索引筛选后读取到的频率值(即余数信息)，它们是一组整数，为频率估计结果，σ₁是相对误差，得到如下表1所示的频率统计结果：

表1无校正情况下高频信号频率测量结果统计

从表1中可以看出，各路的频率信息是一组整数，对高频信号的频率估计存在微小的误差（不超过10），且相对误差σ₁在10^-5数量级。

(2)无噪情况下引入全相位比值谱校正的情况

实验2

仍采用实验1的测量信号、参数及条件，本情况区别仅在于对apFFT频谱做全相位比值内插法对各路的频率值进行校正，以提高测量精度。为高频信号的频率估计结果，σ₂是相对误差，得到如下表2所示的频率统计结果：

表2频谱校正后频率测量结果统计

从表2中可以看出，经过校正后的各路频率信息是一组实数，用这组实数作为余数计算频率在精确到个位的情况下，是不存在误差的，且相对误差在10^-8数量级。

比较表1中的未引入频谱校正的频率估计相对误差σ₁和表2中引入频谱校正的频率估计相对误差σ₂的实验数据，可以发现：对于表1的对各路频率值不做校正的情况，也可以很准确地计算出高频信号的频率信息，其相对测量误差σ₁约处于10^-5数量级，对应的绝对测量误差不超过10Hz。但是对于表2的对各路频率值用全相位比值内插法做校正的情况，其相对测量误差σ₂约处于10^-8数量级，即频率精度普遍提高3个数量级，对应的绝对测量误差仅为约10^-3Hz数量级。

因此，全相位比值内插法是测量精度提高的根本措施。

(3)同一噪声条件下，不同测量频率对象情况

为衡量有噪情况下的测频精度，引入均方根误差（root-mean-square error，RMSE）来度量，均方根误差为

$\mathrm{RMSE}=\sqrt{E\left\{{|{\hat{f}}_0-f_0|}^2\right\}}---\left(12\right)$

实验3

仍采用实验1的测量信号和参数，本情况区别仅在于测量噪声条件，实验1为无噪情况，本实验选取高斯白噪声，其信噪比（Signal to Noise Ratio，SNR）环境设为SNR=10dB。按照图2所示的结合全相位比值谱校正与CRT的频率测量流程，可以得到如下表3所示的频率统计结果：

表3加噪声情况下的频率测量及相对误差结果统计

从3表可以看出，在噪声环境下，本专利所提出的高频测量方法也可以很准确地计算出高频信号的频率信息，其相对测量误差σ_noise约处于10^-7数量级，对应的绝对测量误差不超过10^-2Hz数量级，精度仅仅比无噪声条件降低一个数量级。因此，本测量装置具有很好的抗噪声性能。

(4)不同噪声条件下，同一测量频率对象情况

实验4

仍采用实验1的测量信号和参数，本情况区别仅在于测量噪声条件，实验1为无噪情况，本实验选取高斯白噪声，其信噪比环境变化范围设为SNR=1～50dB，每次加噪计算次数为20次，实验中我们取f₀=2×10⁴Hz。按照图2所示的结合全相位比值谱校正与CRT的频率测量流程，所得到的均方根误差随信噪比（SNR）的变化如图8所示。

从图8可以得出，各种信噪比条件下，即便是在信噪比很低的情况下，本专利所提出的高频测量方法的均方根误差不高于10^-1Hz,相比于10⁴Hz的频率来说，其误差值很小；且在信噪比大约大于30dB以后，其测量结果和在无噪情况下所测结果几乎是吻合的。证明本发明具有很好的抗噪声性能和很高的测频精度。

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

参考文献

[1]崔海涛，黄翔东，蒋长丽.基于全相位FFT的激光测距法[J].计算机工程与应用，2011,47(8s):61-63

[2]J.H.McClellan and C.M.Rader,Number Theory in Digital SignalProcessing.Englewood Cliffs,NJ:Prentice-Hall,1979.

[3]C.Ding,D.Pei,and A.Salomaa,Chinese Remainder Theorem:Ap-plicationsin Computing,Coding,Cryptography.Singapore:World Scientific,1999.

[4]O.Goldreich,D.Ron,and M.Sudan,“Chinese remaindering with er-rors,”IEEE Trans.Inf.Theory,vol.46,no.7,pp.1330–1338,Jul.2000.

[5]V.Guruswami,A.Sahai,and M.Sudan,“Soft-decision decoding of Chineseremainder codes,”inProc.41st IEEE Symp.Foundations Computer Science,RedondoBeach,CA,2000,pp.159–168.

[6]X.-G.Xia and K.Liu,“A generalized Chinese remainder theorem forresidue sets with errors and its application in frequency determination frommultiple sensors with low sampling rates,”IEEE Signal Process.Lett.,vol.12,pp.768–771,Nov.2005.

[7]X.W.Li,H.Liang,and X.-G.Xia,“A robust Chinese remainder the-oremwith its applications in frequency estimation from undersam-pled waveforms,”IEEE Trans.Signal Process.,vol.57,no.11,pp.4314–4322,Nov.2009.

[8]C.Wang,Q.Y.Yin,and W.J.Wang,“An efficient ranging mehtod forwireless sensor networks,”inProc.Int.Conf.Acoustics,Speech,Signal Processing(IEEE ICASSP),Dallas,TX,Mar.2010,pp.2846–2849.

[9]X.-G.Xia and G.Wang,“Phase unwrapping and a robust Chineseremainder theorem,”IEEE Signal Process.Lett.,vol.14,no.4,pp.247–250,Apr.2007.

[10]X.W.Li and X.-G.Xia,“A fast robust Chinese remainder theorembased phase unwrapping algorithm,”IEEE Signal Process.Lett.,vol.15,pp.665–668,Oct.2008.

[11]W.Xu,E.C.Chang,L.K.Kwoh,H.Lim,and W.C.A.Heng,“Phase unwrapping ofSAR interferogram with multi-frequency or multi-base-line,”inProc.IGARSS,1994,pp.730–732.

[12]D.P.Jorgensen,T.R.Shepherd,and A.S.Goldstein,“A dual-pulserepetition frequency scheme for mitigating velocity ambiguities of the NOAAP-3airborne Doppler radar,”J.Atmos.Ocean.Technol.,vol.17,no.5,pp.585–594,May2000.

[13]G.Wang,X.-G.Xia,V.C.Chen,and R.L.Fiedler,“Detection,loca-tion,andimaging of fast moving targets using multifrequency antenna array SAR,”IEEETrans.Aerosp.Electron.Syst.,vol.40,no.1,pp.345–355,Jan.2004.

[14]M.Ruegg,E.Meier,and D.Nuesch,“Capabilities of dual-frequencymillimeter wave SAR with monopulse processing for ground moving targetindication,”IEEE Trans.Geosci.Remote Sens.,vol.45,no.3,pp.539–553,Mar.2007.

[15]Y.M.Zhang and M.Amin,“MIMO radar exploiting narrowband fre-quency-hopping waveforms,”presented at the16th Eur.Signal Pro-cessing Conf.(EUSIPCO),Lausanne,Switzerland,Aug.25–29,2008.

[16]J.Bioucas-Dias,V.Katkovnik,J.Astola,and K.Egiazarian,“Multi-frequency phase unwrapping from noisy data:Adaptive local max-imum likelihoodapproach,”inImage Analysis,Lecture Notes in Com-puter Science.New York:Springer,Jul.2009,vol.5575/2009,pp.310–320.

[17]W.-K.Qi,Y.-W.Dang,and W.-D.Yu,“Deblurring velocity ambiguity ofdistributed space-borne SAR based on Chinese remainder theorem,”J.Electron.Inf.Tech.,vol.31,no.10,pp.2493–2496,Oct.2009.

[18]G.Li,H.Meng,X.-G.Xia,and Y.-N.Peng,“Range and velocity estimationof moving targets using multiple stepped-frequency pulse trains,”Sensors,vol.8,pp.1343–1350,2008.

[19]G.Li,J.Xu,Y.-N.Peng,and X.-G.Xia,“Location and imaging of movingtargets using non-uniform linear antenna array,”IEEETrans.Aerosp.Electron.Syst.,vol.43,no.3,pp.1214–1220,Jul.2007.

[20]X.W.Li and X.-G.Xia,“Multiple-frequency interferomentric velocitySAR location and imaging of elevated moving target,”inProc.Int.Conf.Acoustics,Speech,Signal Processing (IEEE ICASSP),Dallas,TX,Mar.2010,pp.2810–2813.

[21]Wang,W.and X.-G.Xia(2010)."A closed-form robust Chinese remaindertheorem and its performance analysis."Signal Processing,IEEE Transactionson58(11):5655-5666.

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种欠采样速率下的高频余弦信号的频率测量方法，其特征在于，所述方法包括以下步骤：

(1)对高频模拟余弦信号进行过零点检测，取任一过零点作为高频模拟余弦信号中心采样位置；

(2)以过零点为中心，分别以频率f_s1～f_sL对高频模拟余弦信号进行L路低速率采样，每路均采集2M-1个样点，并存储；

(3)对低速率采样得到的L路信号，分别进行加汉宁双窗全相位快速傅里叶变换；

(4)结合过零点类型和apFFT相位谱分布特征，从2L个峰值幅度谱位置中挑选出L个谱位置索引值；

(5)根据余数筛选结果，用全相位比值内插法对各路apFFT谱分别进行谱校正，得到每路的频率估计

(6)利用各路谱校正得到的频率估计值，作为余数，再按照闭合解析形式的中国余数定理对这些余数处理，通过频率估计值重构出原始高频信号的频率f₀；

其中，所述结合过零点类型和apFFT相位谱分布特征，从2L个峰值幅度谱位置中挑选出L个谱位置索引值的步骤具体为：

(a)由过零点过渡情况，确定过零点瞬间相位(+90°或-90°)；

(b)从每路相位谱分布图中，确定与过零点瞬间相位一致的左半边带或右半边带；

(c)从每路确定的半边带中，找出对应的幅度谱峰值位置，作为余数筛选所需的位置索引{k_i ^*,i＝1,…,L}。

2.根据权利要求1所述的一种欠采样速率下的高频余弦信号的频率测量方法，其特征在于，所述根据余数筛选结果，用全相位比值内插法对各路apFFT谱分别进行谱校正，得到每路的频率估计的步骤具体为：

(a)对输入的2N-1样本做阶数N的apFFT谱分析，得到谱估计结果Y(k)，k＝0,…,N-1；

(b)找出Y(k)的峰值谱位置k＝k^*，在振幅谱线中选取相邻最大的两根进行比值，将该比值的平方根记为v，Y(k^*)为谱峰幅值，Y(k^*-1)和Y(k^*+1)分别表示峰值两侧的谱线的幅值，即

$v=\sqrt{\frac{\left|Y\left(k^{*}\right)\right|}{max\left(\right|Y(k^{*}-1)|,|Y(k^{*}+1)\left|\right)}}---\left(2\right)$

(c)根据v求取比例偏差因子Δk，其中，

$\mathrm{\&Delta;}k=\frac{2-v}{1+v}---\left(3\right)$

(d)根据比例偏差因子Δk进行频率校正

若峰值谱线处于k＝k^*的位置上，则校正后的值为

即通过上述步骤(a)至步骤(d)可以得到一路频率估计重复执行可以得到每路的频率估计

3.根据权利要求1所述的一种欠采样速率下的高频余弦信号的频率测量方法，其特征在于，所述方法所能达到的频率测量范围为f_max＝lcm(f_s1，f_s2，…，f_sL)。

4.一种实施权利要求1所述的欠采样速率下的高频余弦信号的频率测量方法的频率测量装置，其特征在于，所述频率测量装置包括：触发电路、模数转化器、DSP器件、输出驱动及其显示电路，

待测信号经过所述触发电路的过零点检测处理来决定其初始相位，然后经过多路采样频率分别为f_s1，f_s2，…，f_sL的所述模数转化器采样得到样本序列{x₁(n)，x₂(n)，…，x_L(n)}，分别以并行数字输入的形式进入所述DSP器件，经过所述DSP器件的处理，得到高频信号的频率估计；最后通过所述输出驱动显示及其显示电路显示出频率值。