CN109308453A - 基于模式聚类与谱校正的欠采样信号频率估计方法及装置 - Google Patents

Abstract

本发明公开了一种基于模式聚类与谱校正的欠采样信号频率估计方法及装置，方法包括：对L路欠采样信号样本序列做N点加汉宁窗的DFT得到DFT频谱，对DFT频谱进行基于比值法的频谱和相位校正，得到校正后的由频率、相位、幅值组成的参数组，并得到由校正后参数组合而成的向量；对于第m个频率成分的第l路独立向量，根据向量集合选出期望的峰值指标；根据谐波参数构成的峰值指标进行模式聚类，得到聚类之后的余数；利用余数构造频率余数组，并将其带入CRT模型进行重构，得出频率估计值。装置包括：模数转化器、DSP芯片、输出驱动及其显示模块。本发明引入频谱校正算法使频率重构精度得到提高，使用模式聚类提高估计器对噪声的鲁棒性。

Description

基于模式聚类与谱校正的欠采样信号频率估计方法及装置

技术领域

本发明涉及信号分析处理技术领域，尤其涉及对于观测接收到的多频信源信号，对其进行欠采样处理并进行频谱校正、模式聚类与中国余数定理重构，从而实现对多频信源的频率估计。

背景技术

含有多成分信号的频率估计是信号处理领域的研究热点，在雷达、无线通信、电子战等领域有着广泛应用。然而欠采样(即采样速率f_s远小于2倍信号频率f₀)下的高精度频率估计一直是个难题，因为奈奎斯特定理要求一个信号周期内采到2个以上样点，随着信号频率升高，必然要求其采样速率相应升高，这就会对模数转换器(Analog to DigitalConverter，简记为ADC)的转换速率、功耗以及硬件成本提出更高要求，在某些特定场合(如f_s＞10⁹samples/s的采样)，甚至是不可实现的。从而只有引入信号处理新方法，才能根本解决该问题。实际上，通过利用一系列附加的模拟电路，如本地振荡发生器、混频器和中频带通器，这些频率可以从低速率ADCs^[1]中提取。然而这些模拟电路会使整个频率估计器的硬件复杂度升高。因此迫切需要开发一种可以直接处理由低速率ADC收集的样本，无需任何附加模拟电路的新型频率估计器。

中国余数定理(Chinese Remainder Categorization,CRT)和DFT(DiscreteFourier Transform)频谱分析的组合是实现欠采样波形频率估计的可行方法^[2-7]。这种组合的基础原理在于：假设N表示要测量的高频率，M₁～M_L为参考的欠采样率(也称为CRT模值)^[2,4,6]。为重构频率N，仅利用一组相互之间满足互素关系的整数模值：M₁～M_L，及相应的余数r₁～r_L(即满足r_i＝N mod M_i,i＝1,...,L,‘mod’表示模除操作)来估计N。因此，一旦r₁～r_L通过DFT确定，所需频率就可以从L个欠采样波形中恢复。通常，这种组合演变为两类：互素谱分析仪和改进的基于CRT的频率估计器。

Vaidyanathan提出了一种使用DFT滤波器的互素谱分析仪^[5,8]。该分析仪使用L＝2个欠采样路径，其中任一个提供滤波后M点DFT(或N点DFT，M和N是互质整数)输出。通过对所有可能的MN个DFT输出进行交叉相关操作，那些具有较大的互相关值的索引(即频率余数r₁,r₂)将被记录下来，并用传统的CRT方法来恢复频率。然而，这种分析器有两个缺点：

(1)统计互相关运算将导致大的延迟和高的复杂度^[5,9]；

(2)恢复结果有限，由于传统的CRT只能提供MN个可能的重建值，虽然文献[10]提出了0.5-shift操作使互素谱分析仪的分辨率加倍，但是其重建精度仍然较低。

近年来，人们提出了许多改进的算法^{[3,4,6,7,11-13]}，其目的是降低CRT的复杂性或提高CRT的重建精度和鲁棒性。这些改进措施大大提高了CRT在频率估计中的适用性。具体而言，降低复杂性是由于CRT重建可以以闭合的方式工作，因此绕过繁琐的搜索操作^[11,12]；提高了分析仪的重建精度和鲁棒性是由于放宽了经典CRT要求L个模值M₁～M_L互素的限制，允许其最大公约数(greatest common divider，简记为gcd)M＝gcd{M₁,...,M_L}为不等于1的正整数^[3,4,6,7,13]。然而，所有这些改进的基于CRT的估计器^{[3,4,6,7,11-13]}只能处理单频欠采样波形的频率估计，不能应用于多频情况。多频检测(假设含有M个频率)的障碍在于将LM个频谱点通过L点DFT频谱归类到M个余数集。

参考文献

[1]R.A.Poisel,Electronic Warfare Receivers and Receiving Systems,Artech House Publishers,2014.

[2]X.-G.Xia,On estimation of multiple frequencies in undersampledcomplex valued waveforms,IEEE Transactions on Signal Processing 47(12)(1999)3417–3419.

[3]W.Wang,X.Li,W.Wang,X.-G.Xia,Maximum likelihood estimation basedrobust Chinese Remainder Theorem for real numbers and its fast algorithm,IEEETransactions on Signal Processing 63(13)(2015)3317–3331.

[4]X.Li,H.Liang,X.-G.Xia,A robust Chinese Remainder Theorem with itsapplications in frequency estimation from undersampled waveforms,IEEETransactions on Signal Processing 57(11)(2009)4314–4322.

[5]P.P.Vaidyanathan,P.Pal,Sparse sensing with co-prime samplers andarrays,IEEE Transactions on Signal Processing 59(2)(2011)573–586.

[6]W.Wang,X.-G.Xia,A closed-form robust Chinese Remainder Theorem andits performance analysis,IEEE Transactions on Signal Processing 58(11)(2010)5655–5666.

[7]L.Xiao,X.G.Xia,A new robust chinese remainder theorem withimproved performance in frequency estimation from undersampled waveforms,Signal Processing 117(2015)242–246.

[8]P.Vaidyanathan,Analog filter banks for sampling:Discretization,polyphase form,and role in compressive sensing,in:Digital Signal Processingand Signal Processing Education Meeting (DSP/SPE),IEEE,2013,pp.24–29.

[9]X.Huang,Z.Yan,S.Jing,H.Fang,L.Xiao,Co-prime sensing-basedfrequency estimation using reduced single-tone snapshots,Circuits,Systems andSignal Processing 35(9)(2016)3355–3366.

[10]X.Huang,Y.Han,Z.Yan,H.Xian,W.Lu,Resolution doubled co-primespectral analyzers for removing spurious peaks,IEEE Transactions on SignalProcessing 64(10)(2016)2489–2498.

[11]X.-G.Xia,G.Wang,Phase unwrapping and a robust Chinese RemainderTheorem,IEEE Signal Processing Letters 14(4)(2007)247.

[12]X.Li,X.-G.Xia,A fast robust Chinese Remainder Theorem based phaseunwrapping algorithm,IEEE Signal Processing Letters(15)(2008)665–668.

[13]G.Li,J.Xu,Y.-N.Peng,X.-G.Xia,An efficient implementation of arobust phase-unwrapping algorithm,IEEE Signal Processing Letters 14(6)(2007)393–396.

[14]F.Zhang,Z.Geng,W.Yuan,The algorithm of interpolating windowed FFTfor harmonic analysis of electric power system,Power Delivery IEEETransactions on 16(2)(2001)160–164.

发明内容

本发明提供了一种基于模式聚类与谱校正的欠采样信号频率估计方法及装置，本发明引入了频谱校正算法使频率重构精度进一步提高，使用模式聚类进一步提高估计器对噪声的鲁棒性，借助中国余数定理实现欠采样情况下的多成分信号的频率估计，详见下文描述：

基于模式聚类与谱校正的欠采样信号频率估计方法，所述方法包括以下步骤：

对L路欠采样信号样本序列做N点加汉宁窗的DFT得到DFT频谱，对DFT频谱进行基于比值法的频谱和相位校正，得到校正后的由频率、相位、幅值组成的参数组，并且得到由校正后参数组合而成的向量；

对于第m个频率成分的第l路独立向量，根据向量集合选出期望的峰值指标；根据谐波参数构成的峰值指标进行模式聚类，从而得到一个聚类之后的余数；

利用余数构造频率余数组，并将其带入CRT模型进行重构，得出频率估计值。

其中，所述根据谐波参数构成的峰值指标进行模式聚类，从而得到一个聚类之后的余数的步骤包括：

1)对每个序列x_l(n)采用加窗的DFT操作，对每类DFT频谱簇进行频谱校正，以便从序列x_l(n)的加汉宁窗DFT结果X_l(k)中精确地提取谐波参数对{a_m,θ_m}；

2)利用向量作为模式聚类指标，得到LM个校正后的余数将其分类到基于CRT的频率重建的M个余数集，从而得到余数集并代入CRT重构模型进行频率估计；

3)频率恢复范围。

进一步地，所述对每类DFT频谱簇进行频谱校正的步骤具体为：

计算出DFT谱峰在X_l(k_m,l)处的幅值和相邻谱线中幅值最大的一根谱线的幅值的比值v：

估计频率偏差

得到校正后的频率估计的幅度估计和相位估计

所述向量具体为：

另一实施例，基于模式聚类与谱校正的欠采样信号频率估计装置，所述装置包括：

将采集到的观测信号x(t)经模数转化器欠采样得到样本序列x(n)，以并行数字输入的形式进入DSP芯片；

经过DSP芯片的内部算法处理，完成信号的加窗DFT变换，频谱校正，模式聚类，利用中国余数定理恢复出频率估计结果；

输出驱动及其显示模块显示频率的估计值。

本发明提出的基于模式聚类与谱校正的多信源欠采样频率估计方法，若将其应用于实际工程领域将具有以下功能优势:

1)降低硬件成本

本发明所提出的联合估计算法，不需要进行统计互相关计算，减少了计算复杂度，同时，本发明打破了传统奈奎斯特采样速率限制，不需要进行下面变频即可实现对高频段信号的采样，进一步降低了硬件要求。

2)提高数据利用率

本发明只需要两个欠采样信号、按照所提算法对其进行处理，即可得出频率估计结果，避免了多次采样与长时间采样。

3)测量精度高

由于本发明采用了频谱校正与模式聚类并结合中国余数定理，具有很高的估计精度。

附图说明

图1为基于模式聚类与谱校正的欠采样信号频率估计方法的流程图；

图2为两个欠采样信号的DFT幅值谱的示意图；

图3为频谱校正前的矢量图；

图4为频谱校正和模式聚类后的矢量图；

图5为本发明的硬件实施图；

图6为DSP的内部流程图。

具体实施方式

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

为了实现欠采样多频波形的频率重构，本发明实施例提出了一种基于谐波参数聚类和中国余数定理的估计器。

首先，通过修正两个加汉宁窗的DFT欠采样波形，可以构造两个精确的谐波参数集。其次，通过对这些谐波参数产生的矢量模式进行聚类，信号中的频率成分得到匹配。最后，利用中国余数定理可以得到最终的频率估计。由于硬件成本低并且适用于多频检测，所提出的估计器在雷达、仪器测量等方面具有巨大的潜力。

实施例1

本发明实施例提供了一种基于模式聚类与谱校正的欠采样信号频率估计方法，参见图1，该方法包括以下步骤：

101：对L路欠采样信号样本序列做N点加汉宁窗的DFT得到DFT频谱；

其中，L路欠采样信号样本序列用x₁(n),...,x_L(n)表示，DFT频谱用X₁(k),...,X_L(k)(k＝0,...,N-1)表示。对任意第l路频谱(l＝0,...,L)，记DFT谱峰点为k_m,l,m＝1,...,M。M的取值为正整数。

102：对DFT频谱进行基于比值法的频谱和相位校正，得到校正后的由频率、相位、幅值组成的参数组，并且得到由校正后参数组合而成的向量；

其中，参数组用表示；校正后参数用表示；向量用表示。

103：对于第m个频率成分的第1路独立向量根据集合选出期望的峰值指标c_m,l；

根据谐波参数构成的峰值指标可以进行模式聚类，从而可以得到一个聚类之后的余数

其中，F₁和F_l为第1路和第l路的欠采样频率，为校正后的频率。

104：利用步骤103中确定的匹配参数(即余数)构造频率余数组(m＝1,...,M)，并将其带入CRT模型进行重构，得出频率估计值

其中，CRT重构模型主要用于单成分信号，本发明实施例在单成分情况的基础上进行扩展，将重构模型用于多成分信号的频率估计。

综上所述，本发明实施例通过上述步骤101-步骤104使用模式聚类进一步提高估计器对噪声的鲁棒性，借助中国余数定理实现欠采样情况下的多成分信号的频率估计。

实施例2

为使本发明的目的、技术方案和优点更加清楚，下面结合具体的数学公式、图1对实施例1中的方案作进一步地详细介绍，详见下文描述：

一、单频估计的CRT重构模型

对于一个单频信号x(t)，

其中，a、θ₀、f₀代表幅度、原始相位和待估频率。

在低硬件成本下估计频率f₀，单频信号x(t)需要用L个欠采样率F₁,...,F_L＜＜f₀进行离散化，并且需要统一欠采样率的最大公约数N＝gcd{f_s1,...,f_sL}，从而得到互素整数集合Γ₁,...,Γ_L：

Γ_l＝F_l/N,l＝1,...,L (4)

从而可以获得L个长度为N的序列x_l(n),l＝1,...L，

其中，F_l是第l路的欠采样率，gcd为取最大公约数操作。

由于F₁,...,F_L＜＜f₀，可以建立以下方程：

其中，n₁,...,n_L是未知折叠整数，ε₁,...,ε_L是L个欠采样序列的归一化频率。理论上，ε_l可以表示为：

ε_l＝(k_l+δ_l)/N,l＝1,...L (7)

其中，k_l∈{0,1,...,N-1}，δ_l是区间[-0.5,0.5)内的小数。因此，CRT所需的理想的L个余数可以表示为：

r_l＝ε_lF_l＝(k_l+δ_l)F_l/N,l＝1,...,L (8)

理想的情况下，折叠整数n₁,...,n_L和待估频率f₀可以通过替换F₁,...,F_L精确恢复。余数组合r₁,...,r_L带入CRT重建程序的详细过程参考文献[6]。

实际上，对于单个序列x_l(n)，整数k_l可以精确地由DFT频谱点索引确定，但是δ_l(也指频率偏移)很可能由于噪声干扰而产生错误。因此，余数r_l的估计为：

需要强调的是，在文献[4,6,7,11,12]中所提出的估计器中忽略了的影响，从而不可避免地产生一些重建误差。

二、多频估计的CRT重构模型

结合单频信号的CRT模型，本发明实施例扩展了多频信号的重构模型，并且进行了精确的频率估计。对于一个多频信号x(t)，

其中，a_m为幅度，f_m为待估频率，θ_m为相位。

对x(t)进行相同欠采样率的采样后得到L个离散序列x_l(n),l＝1,...L：

因此，对于信源中含有M个频率成分，可以建立M个方程：

其中，对于第l个欠采样序列的第m个归一化频率ε_m,l表示为：

ε_m,l＝(k_m,l+δ_m,l)/N,m＝1,...,M,l＝1,...L (13)

其中，δ_m,l为区间[-0.5,0.5)内的小数。

需要强调的是，相比于单频估计，多频估计中余数获取更加困难。对于单频估计，DFT频谱在任何重构路径下都只包含一个谱峰，因此可以很容易地收集L个余数通过CRT重建来恢复f₀。然而，对于多频估计，在任何重构路径上的DFT频谱肯定包含M个谱峰。此外，由于是欠采样，对于第M路独立频率信号，它的谱峰索引k_m,1,...,k_m,L很可能发生变化。因此将LM个谱峰索引通过L点DFT频谱分类到M个余数组合是急需解决的问题。

三、基于谐波参数聚类的余数分类

1)余数采集

对于任意l，将式(12)代入式(11)，得到：

从式(14)中可以得到，对于第m个频率成分，它在所有L个重建路径上的频谱分量共享一个共同的谐波参数对{a_m,θ_m}，这为余数分类提供了基础。然而，受其他频率干扰和频谱泄漏的影响，第m个频率成分的谐波参数对{a_m,θ_m}很可能不能准确地从其谱峰中获取。

下面两个措施被用于提高余数采集精度：

(1)对每个序列x_l(n)采用加窗的DFT操作，而不是直接进行DFT。

其中，采用加窗的DFT操作为本领域技术人员所公知，本发明实施例对此不做赘述。

(2)对每类DFT频谱簇进行频谱校正。

从文献[14]中，可以导出一个基于比值的频谱校正器，以便从序列x_l(n)的加汉宁窗DFT结果X_l(k)中精确地提取谐波参数对{a_m,θ_m}。

其中，上述步骤(2)的过程具体为：

步骤1：计算出DFT谱峰在X_l(k_m,l)处的幅值和相邻谱线(即，用X_l(k_m,l-1)或X_l(k_m,l+1)表示)中幅值最大的一根谱线的幅值的比值v：

步骤2：估计频率偏差

步骤3：得到校正后的频率估计的幅度估计和相位估计

其中，ang[·]指的是取角操作，可以得到当前谱线的相位。

2)基于余数分类的多成分频率估计

一方面，对于第l个重构路径中的第M个频率成分，三重谐波参数实际上是结合在一起的。

另一方面，如上所述，所有L个重构路径的第M个频率成分的频谱分量共享一个公共的谐波参数对{a_m,θ_m}。因此，在任何DFT频谱簇中，校正参数应结合起来去构造一个向量

利用向量作为模式聚类指标，可以得到LM个校正后的余数(m＝1,...,M,l＝1,...L)，将其分类到基于CRT的频率重建的M个余数集，从而得到式(2)中的余数集并代入CRT重构模型进行频率估计。

3)频率恢复范围

正如文献[6]指出的，CRT重建的上段等于所有模的最小公倍数(Least CommonMultiple,LCM)。因此，本发明实施例所提出的估计器的频率恢复范围在0和f_max之间：

其中，F₁和F_L是第1路和第L路的欠采样频率，Γ_l是由式(4)得到的互素整数。

综上所述，本发明实施例引入了频谱校正算法使频率重构精度进一步提高，使用模式聚类进一步提高估计器对噪声的鲁棒性，借助中国余数定理实现欠采样情况下的多成分信号的频率估计。

实施例3

下面结合图2-图4，以及表1对实施例1和2中的方案进行可行性验证，详见下文描述：

在本节中，使用本发明实施例提出的估计器处理欠采样多频复值信号，如式(11)所示，实验中定义M＝4,L＝2和N＝512,Γ₁＝3301,Γ₂＝3307(因此F₁＝1.690112Msamples/s、F₂＝1.693184M samples/s)。

对于多频情况，DFT频谱中含有M个谱峰，而且由于欠采样，对于第M个频率成分，它的谱峰在DFT频谱X₁(k)和X₂(k)中的位置可能会有偏差(如图2所示)。对于相同的信号，幅度和相位参数是配对的，并且上述特征可以用来识别在两个ADC采样序列上相同源的相应峰值位置。同时，为了克服由非零频偏引起的DFT频谱泄漏，采用加窗DFT变换。图2展示出归一化的加窗DFT频谱。

一、频谱校正和聚类算法的性能分析

每个传感器(加性高斯白噪声)的信噪比(SNR)被设置为13dB。图3展示出了没有频谱校正的聚类结果，可以看出结果并不理想。图4是频谱校正和聚类后的矢量图，从图4中可以看出，这两个信号向量匹配成对并聚集成4个方向，频谱校正后聚类方向更加清晰，聚类性能显著提高。

二、频率估计能力

m＝4个源的信号参数和本发明所提出的频率估计器的最终结果显示在表1中。实验结果表明，该算法能够准确估计源信号的频率，且具有较高的估计精度。

表1频率估计的最终结果

实施例4

本发明实施例提供了一种基于模式聚类与谱校正的欠采样信号频率估计装置，参见图5和图6，该装置是与实施例1和2中的方法部分相对应，详见下文描述：

硬件实施图如图5所示，将采集到的观测信号x(t)经过A/D(模数转化器)欠采样得到样本序列x(n)，以并行数字输入的形式进入DSP芯片，经过DSP芯片的内部算法处理，得到频率的估计；最后借助输出驱动及其显示模块显示频率的估计值。

其中，图5的DSP(Digital Signal Processor，数字信号处理器)为核心器件，在信号参数估计过程中，完成如下主要功能：

(1)调用核心算法，完成信号的加窗DFT变换，频谱校正，模式聚类；

(2)利用中国余数定理恢复出频率估计结果，将结果输出至驱动和显示模块。

DSP器件的内部程序流程如图5所示。本发明实施例将所提出的“一种鲁棒的基于模式聚类与谱校正的多信源欠采样频率估计”这一核心估计算法植入DSP器件内，基于此完成高精度、低复杂度、高效的频率估计。

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

1)首先需根据具体应用要求，设置参数初始值；

2)其次，CPU主控器从I/O端口读采样数据，进入内部RAM；

3)最终，按照本发明实施例1和2中的处理过程进行频率估计，并将频率估计结果通过外部显示装置进行显示。

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

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

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

Claims (5)

1.基于模式聚类与谱校正的欠采样信号频率估计方法，其特征在于，所述方法包括以下步骤：

对L路欠采样信号样本序列做N点加汉宁窗的DFT得到DFT频谱，对DFT频谱进行基于比值法的频谱和相位校正，得到校正后的由频率、相位、幅值组成的参数组，并且得到由校正后参数组合而成的向量；

对于第m个频率成分的第l路独立向量，根据向量集合选出期望的峰值指标；根据谐波参数构成的峰值指标进行模式聚类，从而得到一个聚类之后的余数；

利用余数构造频率余数组，并将其带入CRT模型进行重构，得出频率估计值。

2.根据权利要求1所述的基于模式聚类与谱校正的欠采样信号频率估计方法，其特征在于，所述根据谐波参数构成的峰值指标进行模式聚类，从而得到一个聚类之后的余数的步骤包括：

1)对每个序列x_l(n)采用加窗的DFT操作，对每类DFT频谱簇进行频谱校正，以便从序列x_l(n)的加汉宁窗DFT结果X_l(k)中精确地提取谐波参数对{a_m,θ_m}；

2)利用向量作为模式聚类指标，得到LM个校正后的余数将其分类到基于CRT的频率重建的M个余数集，从而得到余数集并代入CRT重构模型进行频率估计；

3)频率恢复范围。

3.根据权利要求2所述的基于模式聚类与谱校正的欠采样信号频率估计方法，其特征在于，所述对每类DFT频谱簇进行频谱校正的步骤具体为：

计算出DFT谱峰在X_l(k_m,l)处的幅值和相邻谱线中幅值最大的一根谱线的幅值的比值v：

估计频率偏差

得到校正后的频率估计的幅度估计和相位估计

4.根据权利要求3所述的基于模式聚类与谱校正的欠采样信号频率估计方法，其特征在于，所述向量具体为：

5.基于模式聚类与谱校正的欠采样信号频率估计装置，其特征在于，所述装置包括：

将采集到的观测信号x(t)经模数转化器欠采样得到样本序列x(n)，以并行数字输入的形式进入DSP芯片；

经过DSP芯片的内部算法处理，完成信号的加窗DFT变换，频谱校正，模式聚类，利用中国余数定理恢复出频率估计结果；

输出驱动及其显示模块显示频率的估计值。