CN104007316A - 一种欠采样速率下的高精度频率测量方法及其测量仪 - Google Patents

Abstract

本发明公开了一种欠采样速率下的高精度频率测量方法及其测量仪，方法包括：用三谱线内插法进行谱校正，得到谱校正后的L对频率估计值，以及峰值位置校正之后的L对相位信息；结合过零点类型和校正之后的L对相位信息，从L对频率估计值中挑选出L个频率估计值；利用各路谱校正和L个频率估计值，作为余数，再按照闭合解析形式的中国余数定理重构出原始高频信号的频率。测量仪包括：待测信号首先经过触发电路的过零点检测处理来决定其初始相位，然后经过多路采样频率的模数转化器采样得到样本序列，分别以并行数字输入的形式进入DSP器件，经过DSP器件的内部处理，得到高频信号的频率估计；最后借助输出驱动及其显示电路显示出频率值。

Description

一种欠采样速率下的高精度频率测量方法及其测量仪

技术领域

本发明涉及数字信号处理领域，尤其涉及一种欠采样速率下的高精度频率测量方法及其测量仪，本发明涉及对高频余弦信号在多路低采样速率下得到的样本，进行高精度测量。

背景技术

欠采样速率下的高精度频率测量仪，一直是频率仪设计的难题。例如：10MHz的正弦波，要想在几路几百kHz甚至几十kHz的采样速率下，测量出频率，在传统频率仪设计中，是很难实现的。因为，为获得信号的频率信息，根据奈奎斯特采样定理，至少要求一个信号周期内采样到2个以上的样点，这样必然要求采样速率高于20MHz以上。仅仅在kHz数量级的采样率下进行采样，会引入信号失真，而无法测出准确频率值。

然而，随着信号频率升高，若机械地应用奈奎斯特采样定理做信号采样，必然会对模数转换器(Analog to Digital Converter，ADC)的转换速率、功耗、硬件成本等方面提出更高的要求。在某些特定场合，甚至是不可实现的。如在软件无线电中，为提高系统的可靠性和灵活性，其基本思想就是将数模转换器(D/A)尽可能地向射频RF(Radio Frequency)靠近，以便于对射频段采集的样本做数字解调、解码等[1]。如果采用高速采样，就必须提高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]提出了余数冗余方法(Remainder Redundancymethod)。其中，对于余数冗余方法非常适合用于低速率采样下的信号频率估计。特别是文献[21]提出了一种闭合形式的CRT重构方法，该方法放宽了经典中国余数定理要求的L个模值互素的要求，仅仅要求L个模值具有某一最大公约数(greatest common divider,gcd)为M即可，即满足M＝gcd{M₁,M₂,…,M_L}，且文献[21]将CRT重构对象从整数域推广到整个实数领域，且该方法对余数误差具有很高的鲁棒性，因而具有很高的实用价值。

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

但是目前还未见有文献解决多路欠采样下的余弦信号的频率估计问题，在当前学术界和工程界，仍属技术空白。例如，为测出f₀＝10MHz的余弦信号频率值，仅仅选择L＝3路低速率采样频率f_s1＝59.392kHz,f_s2＝34.816kHz,f_s3＝38.912kHz,如何从这些采样结果中设计算法方案，准确测量频率。

发明内容

本发明提供了一种欠采样速率下的高精度频率测量方法及其测量仪，本发明实现了多路低速率欠采样下的高频余弦信号的频率测量，详见下文描述：

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

(1)对低速率采样得到的L路信号，分别进行快速傅里叶变换得到幅度谱以及相位谱；

(2)找出各路幅度谱中的峰值位置，用三谱线内插法进行谱校正，得到谱校正后的L对频率估计值，以及峰值位置校正之后的L对相位信息；

(3)结合过零点类型和校正之后的L对相位信息，从L对频率估计值中挑选出L个频率估计值；

(4)利用各路谱校正和L个频率估计值，作为余数，再按照闭合解析形式的中国余数定理对这些余数处理，以重构出原始高频信号的频率f₀。

所述用三谱线内插法进行谱校正，得到谱校正后的L对频率估计值，以及峰值位置校正之后的L对相位信息的步骤具体为：

1)利用幅度谱中峰值位置以及与峰值位置相邻的左右旁瓣的两根谱线的幅度信息求取频率偏差值△k，若峰值谱线处于k＝k^*的位置上，k^*表示峰值位置，F(k^*)为谱峰幅值，F(k^*-1)和F(k^*+1)分别表示峰值两侧的谱线的幅值，Real表示取实部，则

$\mathrm{\&Delta;k}=\frac{\tan (\mathrm{\&pi;}/M)}{\mathrm{\&pi;}/M}\mathrm{Real}\left\{\frac{F(k^{*}-1)-F(k^{*}+1)}{2F\left(k^{*}\right)-F(k^{*}-1)-F(k^{*}+1)}\right\}$

2)根据频率偏差值△k进行频率校正和相位校正，获得校正后的频率估计值和相位信息值分别为

$\hat{f}=(k^{*}+\mathrm{\&Delta;k})f_s/M$

所述结合过零点类型和校正之后的L对相位信息，从L对频率估计值中挑选出L个频率估计值的步骤具体为：

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

2)从每路提供的两个校正之后的相位值中，筛选出其中与过零点瞬间相位正负符号一致的相位值；

3)将每路筛选出的相位值所对应的半边频带的频率估计值作为CRT所需的余数。

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

待测信号首先经过所述触发电路的过零点检测处理来决定其初始相位，然后经过多路采样频率分别为f_s1，f_s2，…，f_sL的所述模数转化器采样得到样本序列，分别以并行数字输入的形式进入所述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＝4096Hz，f_s2＝10240Hz，f_s3＝14336Hz，各路的测量范围分别为(0,2048Hz]，(0,5120Hz]，(0,7168Hz]。对于本发明，则根据中国余数定理，最大可测频率为f_max＝lcm(f_s1，f_s2，f_s3)＝1.4336×10⁵Hz，则测量范围为(0,1.4336×10⁵Hz]，测量范围提高了将近两个数量级。

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

相对而言，如果采用以往的高速采样，就必须提高A/D采样速率，必然会为此付出更高的功耗与硬件成本作为代价。例如，实验1中各路低速率采样频率分别为：f_s1＝4096Hz，f_s2＝10240Hz，f_s3＝14336Hz，则测量范围为(0,1.4336×10⁵Hz]。同样地，如果进行高速率采样，要达到同样的测量范围，要求的采样频率至少为f_s＝2.8672×10⁵Hz。因此本发明大大降低了采样速率，节省了硬件成本。

第三、在无噪情况下，本发明方法对高频信号的测量误差几乎完全可忽略。

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

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

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

附图说明

图1为欠采样速率下的高精度频率测量方法的流程图；

图2为过零点检测示例图；

图3为举例说明每路采样信号的FFT分析谱图；

图4为不同信噪比下频率估计的均方根误差图；

图5为频率估计的方差统计图；

图6为欠采样速率下的高精度频率测量仪的硬件实施图；

图7为DSP内部程序流图。

具体实施方式

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

对于解决多路欠采样下的余弦信号的频率估计问题的测频难度在于：余弦信号包含两个复指数边带，故对于某一频率的余弦信号，从L路低速采样做FFT谱峰搜索，必然得到2L路余数，怎样从2L个余数中挑选出CRT重构所需的有效的L路余数是个非常棘手的问题。本发明实施例将相位信息作为筛选的依据来解决该问题，并结合闭合形式的中国余数定理，完成欠采样下的余弦信号的频率估计，以填补该空白，参见图1，详见下文描述：

101：对高频模拟余弦信号进行过零点检测，取任一过零点作为高频模拟余弦信号的起始采样位置，并将选取的过零点作为过零检测点；

该步骤具体为：对输入的高频模拟余弦信号x(t)进行过零点检测，对于余弦信号而言，过零点存在如图2所示的两种情况：图2(a)对应为从负波形到正波形过零，这时过零点的瞬间相位为-π/2；图2(b)对应为从正波形到负波形过零，这时过零点的瞬间相位为π/2。过零点处-π/2或π/2的瞬间相位的符号，对于后面步骤的每路信号的FFT峰值位置中选取一个作为余数，起到决定作用。模拟信号经过简单的触发电路可以很容易地确定过零点时刻。

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

即，以过零检测点为起始点采集M点，一共采集L路信号。对于高频模拟信号x(t)＝a·cos(2πf₀t-π/2)，所测频率为f₀，采样频率分别为f_s1～f_sL，则各路采样信号为：

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

其中n＝0,…,M-1，采样速率要求f_s1～f_sL满足具有公约数M，且除以公约数M后是两两互素的。

103：对低速率采样得到的L路信号，分别进行快速傅里叶变换(FFT)，得到幅度谱记为{F_i(k),i＝0,…,L-1；k＝0,…,M-1,}以及相位谱(其中i表示采样路序号，k表示谱线标号)；

该步骤具体为：对各路采样信号进行M点快速傅里叶变换操作，得到对各路低速率采样信号的谱分析结果。

举例如下：直接对x(t)＝a·cos(2πf₀t-π/2)做各路低速欠采样信号，进行FFT操作而取峰值处的频率值，所得到的是一组整数。获得的FFT谱如图3所示(L＝3路采样，a＝2，高速采样频率f₀＝998.4000Hz,低速率采样频率f_s1＝128Hz,f_s2＝192Hz,f_s3＝320Hz,易推出M＝gcd{128,192,320}＝64，“gcd”表示最大公约数。)

图3的FFT幅度谱F_i(k)和相位谱具有如下规律

1)每路幅度谱均有两个谱峰，其位置关于频率轴中心对称。

2)每路信号两个谱峰位置对应的相位谱出现跳变，且其相位谱值是大小相等、符号相反的(例如第1路，谱峰位置为k＝13和k＝51两处，对应的相位谱

3)每路信号的峰值位置的相位信息紊乱，并非真实的相位信息，因此需要对相位信息进行校正。

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

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

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

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

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

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

104：找出各路FFT幅度谱中的峰值位置，每路有两个对称的峰值位置，用三谱线内插法对各路FFT谱分别进行谱校正(包括频率校正和相位校正)，得到谱校正后的L对频率估计值，分别为：(i＝1,...,L)，以及峰值位置校正之后的L对相位信息：(i＝1,...,L)。

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

所以本发明提出“三谱线内插法”用以提高CRT余数精度，该方法需对峰值谱和泄漏出的左右旁谱的两根谱线进行校正，从而可将峰值位置的频率估计值(每路左、右半边频带各包含1个频率校正值)精确到小数。并且为了能够从这L对频率估计值中筛选出L个作为余数，需要以准确峰值位置的相位估计作依据，而峰值位置的相位估计是紊乱的，因此每路的峰值位置的相位估计也需要进行校正，得到(i＝1,...,L)。

“三线内插法”的处理步骤如下：

1)利用幅度谱中峰值位置以及与峰值位置相邻的左右旁瓣的两根谱线的幅度信息求取频率偏差值△k，若峰值谱线处于k＝k^*的位置上，k^*表示峰值位置，F(k^*)为谱峰幅值，F(k^*-1)和F(k^*+1)分别表示峰值两侧的谱线的幅值，Real表示取实部，则

$\mathrm{\&Delta;k}=\frac{\tan (\mathrm{\&pi;}/M)}{\mathrm{\&pi;}/M}\mathrm{Real}\left\{\frac{F(k^{*}-1)-F(k^{*}+1)}{2F\left(k^{*}\right)-F(k^{*}-1)-F(k^{*}+1)}\right\}---\left(2\right)$

2)根据频率偏差值△k进行频率校正和相位校正，获得校正后的和值分别为

$\hat{f}=(k^{*}+\mathrm{\&Delta;k})f_s/M---\left(3\right)$

105：结合过零点类型和经过三谱线内插法校正之后的相位信息，从上述得到的L对频率估计值中挑选出L个频率估计值，得到

该步骤具体为：

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

2)从每路提供的两个校正之后的相位值中，筛选出其中与过零点瞬间相位正负符号一致的相位值；

3)将每路筛选出的相位值所对应的半边频带(左半频带或右半频带)的频率校正值作为CRT所需的余数。

例如：本例中，过零点是从负值到正值的过零点，其瞬间相位为-90°，经过三谱线内插法校正之后的峰值位置的频率和相位估计如表1所示。

表1 校正之后的频率和相位估计(过零点瞬间相位：-90°)

从表1可以很明显地看出与瞬时相位-90°符号相等的相位估计值分别为： 故从对应的半边频带筛选出的频率估计值分别为： ${\hat{f}}_{1,R}=102.3997\mathrm{Hz},$ ${\hat{f}}_{2,L}=38.4005\mathrm{Hz},$ ${\hat{f}}_{3,L}=38.3969\mathrm{Hz},$ 以此作为余数进行下面的操作。

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

上述步骤中，要求f_s1～f_sL为整数，其最大公约数(gcd)为M，且除以公约数M后是两两互素的。上述所得出的各路信号的频率值即为中国余数定理中所需的余数。将各路采样频率f_s1，f_s2，…，f_sL作为CRT的各路模值，即结合最大公约数M＝gcd{f_s1，f_s2，…，f_sL作为CRT的各路模值，从而CRT所需L个的互素整数为Γ_i＝f_si/M,i＝1,...,L。则按照如下闭合解析形式的中国余数定理的算法步骤，估计出高频信号频率值，具体为：

1)从所给的余数计算其中：

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

2)计算模除Γ_i的余数：

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

其中，为Γ₁关于Γ_i的模逆，可以提前算出，Γ₁和Γ_i是前文提到的满足互素的整数。

3)计算

${\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的模逆，且γ₁由γ_i＝Γ₁LΓ_i-1Γ_i+1＝Γ/Γ_i定义。

4)计算(2≤i≤L)

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

5)计算

由上述3)和4)得到的(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)$

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

例如：将表1余数筛选出的3个频率估计值分别为： 结合最大公约数M＝64，及L个模值f_s1＝128Hz,f_s2＝192Hz,f_s3＝320Hz，代入上述步骤可获得高频频率估计值为998.4000000000001Hz,与真实值998.4Hz的误差几乎可以忽略。

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

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

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

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

实验1

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

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

表2 不进行谱校正情况下的频率测量结果统计

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

(2)无噪情况下引入三线谱校正的情况

实验2

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

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

从表3中可以看出，经过三线内插法谱校正后的各路频率值可精确到小数位，从而提高了CRT所需的余数精度，基于此可获得更高精度的频率估计值从表3可看出，的估计精度至少可高达10^-5Hz级，对应的相对误差β₂处于10^-10数量级以上。

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

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

实验3

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

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

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

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

(a)频率测量的均方根误差估计

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

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

实验4

仍采用实验1的测量信号和参数，本情况区别仅在于测量噪声条件，实验1为无噪情况，本实验选取高斯白噪声，其信噪比环境变化范围设为SNR＝1～100dB，每次加噪计算次数为20次，实验中我们取f₀＝1.0×10⁵Hz。按照图1所示的结合三线内插谱校正与CRT的频率测量流程，所得到的均方根误差随信噪比(SNR)的变化如图4所示：

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

(b)频率测量的方差估计

为衡量有噪情况下的测频精度，引入测量方差来度量，对于每路信号的频率信息，其理论的方差表达式为(i＝1,…，L)

$\mathrm{var}\left({\hat{f}}_i\right)=\frac{[{\tan }^2\left(\frac{\mathrm{\&pi;}}{M}\right)+3{\tan }^2\left(\frac{\mathrm{\&pi;}{\mathrm{\&delta;}}_i}{M}\right)]}{2{\mathrm{\&pi;}}^2\mathrm{gMg\&rho;}}---\left(13\right)$

其中，δ_i为每路采样信号的频偏，M为采样点数，为第i路采样信号的校正后的频率估计值，ρ为信噪比表达形式，如下式所示

ρ＝10^SNR/10                       (14)

其中SNR是以分贝(dB)为单位的信噪比。对于余弦信号的频率估计的理论方差，可推导出其表达式为

$\mathrm{var}\left(\hat{f}\right)=\frac{1}{L^2}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}\mathrm{var}\left({\hat{f}}_i\right)g{f_{\mathrm{si}}}^2=\frac{1}{L^2}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}\frac{[{\tan }^2\left(\frac{\mathrm{\&pi;}}{M}\right)+3{\tan }^2\left(\frac{\mathrm{\&pi;}{\mathrm{\&delta;}}_i}{M}\right)]{f_{\mathrm{si}}}^2}{2{\mathrm{\&pi;}}^2\mathrm{gMg\&rho;}}---\left(15\right)$

其中L为采样路数，f_si(i＝1,…，L)为各路的采样频率。

实验5

仍采用实验1的测量信号和参数，本情况区别仅在于测量噪声条件，实验1为无噪情况，实验中取f₀＝9.0317×10⁴Hz，本实验选取高斯白噪声，其信噪比环境变化范围设为SNR＝1～50dB，每次加噪计算次数为200次。按照图1所示的结合三线内插谱校正与CRT的频率测量流程，所得到的频率估计的方差的变化如图5所示：

从图5的统计结果可以看出，频率估计的实测统计方差曲线与式(15)得到的的理论方差曲线几乎是吻合的，验证了式(15)的理论方差表达式的正确性。

一种欠采样速率下的高精度频率测量仪，参见图6，包括：触发电路，模数转化器，DSP器件以及输出驱动及其显示电路，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

参考文献

[1]杨小牛，楼才义，徐建良.软件无线电原理与应用[M].北京：电子工业出版社,2001.

[2]J.H.McClellan and C.M.Rader,Number Theory in Digital Signal Processing.EnglewoodCliffs,NJ:Prentice-Hall,1979.

[3]C.Ding,D.Pei,and A.Salomaa,Chinese Remainder Theorem:Ap-plications in 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 Chinese remainder codes,”inProc.41st IEEE Symp.Foundations Computer Science,Redondo Beach,CA,2000,pp.159–168.

[6]X.-G.Xia and K.Liu,“A generalized Chinese remainder theorem for residue sets with errorsand its application in frequency determination from multiple sensors with low samplingrates,”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-orem with itsapplications in frequency estimation from undersam-pled waveforms,”IEEE Trans.SignalProcess.,vol.57,no.11,pp.4314–4322,Nov.2009.

[8]C.Wang,Q.Y.Yin,and W.J.Wang,“An efficient ranging mehtod for wireless sensornetworks,”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 Chinese remainder theorem,”IEEESignal Process.Lett.,vol.14,no.4,pp.247–250,Apr.2007.

[10]X.W.Li and X.-G.Xia,“A fast robust Chinese remainder theorem based phase unwrappingalgorithm,”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 of SARinterferogram 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-pulse repetition frequencyscheme for mitigating velocity ambiguities of the NOAA P-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,and imaging offast moving targets using multifrequency antenna array SAR,”IEEE Trans.Aerosp.Electron.Syst.,vol.40,no.1,pp.345–355,Jan.2004.

[14]M.Ruegg,E.Meier,and D.Nuesch,“Capabilities of dual-frequency millimeter wave SARwith monopulse processing for ground moving target indication,”IEEE Trans.Geosci.RemoteSens.,vol.45,no.3,pp.539–553,Mar.2007.

[15]Y.M.Zhang and M.Amin,“MIMO radar exploiting narrowband fre-quency-hoppingwaveforms,”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 phaseunwrapping from noisy data:Adaptive local max-imum likelihood approach,”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 of distributedspace-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 estimation of movingtargets 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 moving targets usingnon-uniform linear antenna array,”IEEE Trans.Aerosp.Electron.Syst.,vol.43,no.3,pp.1214–1220,Jul.2007.

[20]X.W.Li and X.-G.Xia,“Multiple-frequency interferomentric velocity SAR location andimaging 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 remainder theorem and itsperformance analysis."Signal Processing,IEEE Transactions on58(11):5655-5666.

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

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

Claims (4)

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

(1)对低速率采样得到的L路信号，分别进行快速傅里叶变换得到幅度谱以及相位谱；

(2)找出各路幅度谱中的峰值位置，用三谱线内插法进行谱校正，得到谱校正后的L对频率估计值，以及峰值位置校正之后的L对相位信息；

(3)结合过零点类型和校正之后的L对相位信息，从L对频率估计值中挑选出L个频率估计值；

(4)利用各路谱校正和L个频率估计值，作为余数，再按照闭合解析形式的中国余数定理对这些余数处理，以重构出原始高频信号的频率f₀。

2.根据权利要求1所述的一种欠采样速率下的高精度频率测量方法，其特征在于，所述用三谱线内插法进行谱校正，得到谱校正后的L对频率估计值，以及峰值位置校正之后的L对相位信息的步骤具体为：

1)利用幅度谱中峰值位置以及与峰值位置相邻的左右旁瓣的两根谱线的幅度信息求取频率偏差值△k，若峰值谱线处于k＝k^*的位置上，k^*表示峰值位置，F(k^*)为谱峰幅值，F(k^*-1)和F(k^*+1)分别表示峰值两侧的谱线的幅值，Real表示取实部，则

$\mathrm{\&Delta;k}=\frac{\tan (\mathrm{\&pi;}/M)}{\mathrm{\&pi;}/M}\mathrm{Real}\left\{\frac{F(k^{*}-1)-F(k^{*}+1)}{2F\left(k^{*}\right)-F(k^{*}-1)-F(k^{*}+1)}\right\}$

2)根据频率偏差值△k进行频率校正和相位校正，获得校正后的频率估计值和相位信息值分别为

$\hat{f}=(k^{*}+\mathrm{\&Delta;k})f_s/M$

3.根据权利要求1所述的一种欠采样速率下的高精度频率测量方法，其特征在于，所述结合过零点类型和校正之后的L对相位信息，从L对频率估计值中挑选出L个频率估计值的步骤具体为：

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

2)从每路提供的两个校正之后的相位值中，筛选出其中与过零点瞬间相位正负符号一致的相位值；

3)将每路筛选出的相位值所对应的半边频带的频率估计值作为CRT所需的余数。

4.一种欠采样速率下的高精度频率测量装置，所述高精度频率测量装置包括：触发电路，模数转化器，DSP器件以及输出驱动及其显示电路，其特征在于，

待测信号首先经过所述触发电路的过零点检测处理来决定其初始相位，然后经过多路采样频率分别为f_s1，f_s2，…，f_sL的所述模数转化器采样得到样本序列，分别以并行数字输入的形式进入所述DSP器件，经过所述DSP器件的内部处理，得到高频信号的频率估计；

最后借助所述输出驱动及其显示电路显示出频率值。