CN104991119A - 一种消除伪峰、谱泄漏效应的互素谱分析方法及其装置 - Google Patents

Abstract

本发明公开了一种消除伪峰、谱泄漏效应的互素谱分析方法及其装置，互素谱分析方法包括以下步骤：对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置。装置包括：数字信号处理器，及输出驱动及显示电路。本发明实现了完全消除伪峰效应和谱泄漏效应，准确地识别出给定的频率，提高了频率分辨率。

Description

一种消除伪峰、谱泄漏效应的互素谱分析方法及其装置

技术领域

本发明涉及数字信号处理领域，尤其涉及一种消除伪峰、谱泄漏效应的互素谱分析方法及其装置，具体涉及信号欠采样、多相滤波过程、信号相关性及功率谱密度估计。

背景技术

谱估计是与许多实际物理量测量(如大气动力学研究中的雷达风速测量^[1]、振动分析中的转速测量^[2]等)紧密相关的信号处理方法^[3,4]，在电子学、电气学和自动化等各个学科都有广泛的应用，是数字信号处理的基础问题。传统谱估计只有满足香农采样定理(即要求一个信号周期内至少采集2个样点)，信号才不会发生失真，才有可能对信号进行有效的恢复。然而，在许多工程应用(如雷达^[5]、阵列波达方向(DOA)估计^[6]等领域)中，受到现场环境、采样设备等条件限制(如最高模数转换速率限制)，难以满足香农采样，并且随着信号频率的提高采样难度和设备成本也会大幅度提高，因而稀疏样本的频谱估计(采样速率fs远小于2倍信号频率f0，即欠采样)成为学术界和工程界迫切需要解决的问题。

针对此问题，国内外学者提出了很多稀疏谱估计方法，如文献[7-9]提出将中国余数定理(Chinese Remainder Theorem,CRT)与离散傅里叶变换(Discrete Fourier Transform,DFT)相结合的频谱测量方法，该方法通过对高频信号进行不同点数的DFT，求取多路欠采样信号的谱余数，将高频信号的问题转换成余数低频信号的问题，利用得到的谱余数来重构未知的高频频率，降低了信号采样难度和对采样设备的要求，然而该方法目前仅仅能实现单频测量，无法胜任多频测量情况；文献[10]提出基于样本嵌入(Nested sample)的测量算法，该方法要求对一个信号进行两种不同速率的采样，将同一信号的高速率样本嵌入到低速率样本中，利用样本位置之间的差分阵列(difference-co-array)增加自由度(freedoms)的性质，通过导出两路样本的互相关特性来获得信号的功率谱估计，该方法在整体上降低了信号采样率，并且对于较密集的频率恢复也有较好的效果，然而该方法在局部无法摆脱高速采样器，仍需进行一定的高速采样，因而它并不是真正意义上的稀疏谱估计方法。因此，找到真正意义上的稀疏谱估计方法，彻底摆脱高速采样器的约束，并且可以适用于更一般的多频信号是一个亟待解决的难题。

为解决稀疏谱估计问题，近年来，一种新型的谱估计方法—互素感知(co-prime sensing)理论^[11-14]受到越来越多的关注，该方法首先对单个模拟输入信号作两路并行的稀疏采样(要求其两个欠采样速率的数值满足互素关系)，然后对得到的两路稀疏样本分别作多相滤波，将单通道信号转换成多通路子信号，再同时对两路并行的多相滤波输出信号做IDFT，最后对IDFT输出的多通道信号做互相关扫描以估计信号频谱位置，从而获得真正意义上的高分辨率谱。该方法结构简单，算法清晰，用到的都是数字信号的基础原理，在一定程度上解决了稀疏谱估计问题。

但文献[11]提出的方法也有自身难以克服的缺陷，为了完成谱估计，需要设计两路低通原型滤波器，通带截止频率分别为π/M和π/N(M、N是互素的整数)，但由于理想滤波器(过渡带为0，截止频率准确落在期望点上)不可能实现(实际的滤波器总存在过渡带、截止频率位置不精确)，因而对于处在滤波器通带中间位置的信号成分(即信号频率接近iΔω，i∈Z⁺，Δω＝2π/(MN))，文献[11]提出的互素谱分析方法可以很好地将之识别出来；但对于接近通带边缘的信号成分(即信号频率接近(i+0.5)Δω，i∈Z⁺)，由于非理想子滤波器相邻边带之间的相互影响，会使得谱分析结果在相邻位置出现泄漏，在其他位置出现伪峰，从而影响最终的信号频率识别，特别是伪峰现象的存在，严重影响了互素谱的可读性，在军事应用中会造成虚警。由于实际信号的复杂性，这种情况很可能大量出现，严重影响频谱估计的性能。因而针对文献[11]提出的互素谱分析器因理想滤波器的不可实现性而造成的伪峰和谱泄漏问题，急需对其互素谱分析结构做必要的改进，对其谱分析算法做优化。

发明内容

本发明提供了一种消除伪峰、谱泄漏效应的互素谱分析方法及其装置，本发明实现了完全消除伪峰效应和谱泄漏效应，准确地识别出给定的频率，提高了频率分辨率，详见下文描述：

一种消除伪峰、谱泄漏效应的互素谱分析方法，所述互素谱分析方法包括以下步骤：

对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；

平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；

将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置。

所述对离散傅里叶变换后的输出序列进行互相关扫描求值的步骤之前，所述方法还包括：

对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

所述对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列的步骤具体为：

1)对输入信号进行两路下采样，得到两路稀疏信号；

2)对两路稀疏信号分别进行直接多相滤波，以及分别加平移因子后的多相滤波，输出两对多通道稀疏信号；

3)在每个时刻，对多通道输出序列分别进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

一种消除伪峰、谱泄漏效应的互素谱分析装置，所述装置包括：数字信号处理器，及输出驱动及显示电路，

所述数字信号处理器用于对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置；

所述输出驱动及显示电路用于输出归一化频率在频谱图中的位置。

所述数字信号处理器包括：第一获取模块，

所述第一获取模块，用于对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

所述数字信号处理器还包括：第二获取模块、输出模块和第三获取模块，

所述第二获取模块，用于对输入信号进行两路下采样，得到两路稀疏信号；

所述输出模块，用于对两路稀疏信号分别进行直接多相滤波，以及分别加平移因子后的多相滤波，输出两对多通道稀疏信号；

所述第三获取模块，用于在每个时刻，对多通道输出序列分别进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

本发明提出的可消除伪峰效应和谱泄漏效应的互素谱分析方法，若用于频率估计及实际工程领域，可产生如下有益效果：

第一、消除了原互素谱存在的伪峰效应和谱泄漏效应。

对于传统互素感知频谱估计方法，为了完成谱估计，需要设计两路低通原型滤波器，但由于实际滤波器的非理想性(不可避免的过渡带带宽、截止频率位置不精确)，对于滤波器通带中间位置的信号成分，即信号频率接近iΔω，Δω＝2π/(MN)，能很好的识别出来，对于接近通带边缘的信号成分，即信号频率接近(i+0.5)Δω，受到非理想子滤波器相邻边带的相互干扰及其互素谱输出互相关扫描的影响，如前所述，原互素谱不但会在期望主谱线的左边或右边泄漏出旁谱线，而且会在离期望主谱线较远的位置产生伪峰，从而严重影响频谱估计的性能。本发明提出具有互补性质的0.5平移的互素谱估计器，使得对于每一个频率，至少有一种情况不会出现频谱泄漏，将它与原互素谱估计器结合起来，加入频率指示器和逻辑判决器，可完全消除伪峰效应和谱泄漏效应。

例如实验1和实验2中的频率估计过程，无论是单频还是多频情况，最终都能完全消除伪峰效应和谱泄漏效应，准确地识别出给定的频率。

第二、相比于原互素谱，其频率分辨率提高了1倍。

传统的互素感知频谱估计方法中，对于频率最终重构的频率只能四舍五入到整数，人为地忽略了信号频率的小数偏移，尤其是对于小数偏移接近0.5的频率，会引起较大误差，而本发明提出的谱估计方法由于引入了频率指示器，对于小数偏移接近0.5的频率，可以将真实谱线附近的相邻两根谱线提取出来，再将其定位到两根谱线的中间位置，即可将0.5偏移识别出来，这样使得频率分别率提高了1倍。

由实验1和实验2可以看出，可以将频率为10.51Hz和7.51Hz的频率识别为10.5Hz和7.5Hz而不是四舍五入到11Hz和8Hz，估计精度更高。

第三、适于在低速率应用场合进行频率场合。

这是由采用互素谱感知的结构决定的，由于图1的互素谱感知采用了两级下采样结构(互素下采样和多相滤波下采样)，且IDFT、频率指示、逻辑判决均是在两级下采样后工作，故非常适合于低数据率应用场合，非常有利于降低成本。

附图说明

图1为谱泄漏抑制的高分辨率稀疏谱估计器设计流程图；

图2为多相滤波结构图(u路未加平移)；

图3为单频信号互素感知谱估计示意图；

图4为加入平移因子的多相滤波结构图；

图5为原通道与平移通道互素感知频谱对比图；

(a)δ＝0.03时原通道的频谱图；(b)δ＝0.47时原通道的频谱图；(c)δ＝0.03时平移通道的频谱图；(d)δ＝0.47时平移通道的频谱图。

图6为逻辑判决器算法流程图；

图7为单频频谱估计示意图；

(a)原互素感知频谱估计；(b)平移通路互素感知频谱估计；(c)平移相乘的频率指示；(d)直接相乘频率指示；(e)消除伪峰和谱泄漏的频谱分布。

图8为多频频谱估计示意图；

(a)原互素感知频谱估计；(b)平移通路互素感知频谱估计；(c)平移相乘的频率指示；(d)直接相乘频率指示；(e)消除伪峰和谱泄漏的频谱分布。

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

图10为DSP内部程序流图。

具体实施方式

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

针对以上不足，本发明提出一种互补型的互素谱分析器，它是在原分析器的低速率端加上一个等效的0.5Δω平移得到的，这样，对于任意信号成分，至少有一种谱分析器(原始谱分析器或互补谱分析器)不会产生伪峰和谱泄漏，然后对两个输出结果做逻辑与运算，就可以将泄漏出来的谱线屏蔽掉，再经过一些必要的判定处理，最终将所需的频谱提取出来。这种方法通过两种谱估计器对不同信号成分分析时的互补特性，绕过了滤波器的非理想性对谱估计的影响，从而从根本上消除了伪峰问题。

本发明可以解决现有互素感知理论的谱伪峰问题，在稀疏样本的雷达、波达方向估计和波束成型等领域有广阔的应用前景。

图1中，本发明的处理分三个部分，包括信号的互素感知处理、频率指示器以及逻辑判决器，详见下文描述：

101：对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列，对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；

1)对输入信号x(n)进行两路下采样，下采样因子分别为M和N(M和N为任意的互素的整数，具体的取值根据实际应用中的需要进行设定)，得到两路稀疏信号x_u(n)和x_v(n)；

2)对两路稀疏信号x_u(n)和x_v(n)分别进行直接多相滤波，以及分别加平移因子λ_u和λ_v后的多相滤波，输出两对多通道稀疏信号；

其中，一对多通道稀疏信号为x_u,p(n)与另一对多通道稀疏信号为x_v,q(n)与0≤p≤M-1，0≤q≤N-1。

3)在每个时刻n，对多通道输出序列{x_u,p(n),p＝0,...,M-1}和{p＝0,...,M-1}分别进行M点IDFT(离散傅里叶变换)，得到离散傅里叶变换后的输出序列{u_k(n),k＝0,...,M-1}和{k＝0,...,M-1}。

即，多通道输出序列由多通道稀疏信号组成。

类似地，对多通道输出序列{x_v,q(n),q＝0,...,N-1}和{q＝0,...,N-1}分别进行N点IDFT，得到离散傅里叶变换后的输出序列{v_l(n),l＝0,...,N-1}和{l＝0,...,N-1}。

4)对离散傅里叶变换后的输出序列和{k＝0,...,M-1}，{v_l(n),l＝0,...,N-1}和{l＝0,...,N-1}分别进行互相关扫描求值，得到原频谱序列和平移后的频谱序列

即，原频谱序列是由{u_k(n),k＝0,...,M-1}和{v_l(n),l＝0,...,N-1}进行互相关扫描求值得到的；平移后的频谱序列是由{k＝0,...,M-1}和{l＝0,...,N-1}进行互相关扫描求值得到的；ω_i表示在归一化频率位置i处的角频率，即ω_i＝iΔω，Δω为设定的频率分辨率。

102：平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；

1)对平移后(对时间进行操作)的频谱序列向左频移(对频率进行操作)一个单位，得到频移后的频谱序列

2)将原频谱序列与频移后的频谱序列对应元素相乘，得到具有频率指示作用的第一频率指示序列ρ₁(i)；

将原频谱序列和平移后的频谱序列对应元素相乘，得到具有频率指示作用的第二频率指示序列ρ₂(i)，其中i为归一化频率在频谱图中的位置，0≤i≤MN-1。

103：将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置。

通过上述的步骤101-103，本发明实现了完全消除伪峰效应和谱泄漏效应，准确地识别出给定的频率，提高了频率分辨率。

实施例2

本实施例2结合说明书附图、具体的计算公式对实施例1中的步骤进行详细的介绍，详见下文描述：

201：信号互素感知处理；

由于多频信号的互素谱感知最终可分解为多个单频信号互素谱感知分别进行处理。因而，先以单频信号为例，来分析信号经过互素谱分析中的各个节点的过程，然后再将之推广到多频情况。

(1)直接通路的互素谱感知信号流程分析；

单频信号中，a₀是信号幅度，f₀是归一化频率，f₀＝(m+δ)Δf，Δf＝1/(MN)，m∈[0，MN-1)，δ是小数频率偏移(-0.5≤δ＜0.5)，θ₀是初始相位。

对于u路处理，经过N倍下采样后的信号为x_u(n)，有如下形式

$x_u\left(n\right)=x\left(Nn\right)=a_0\·e^{j(2{\mathrm{\πf}}_{0}Nn+{\mathrm{\θ}}_0)}---\left(1\right)$

同理，v路M倍下采样后的信号x_v(n)有如下形式

$x_v\left(n\right)=x\left(Mn\right)=a_0\·e^{j(2{\mathrm{\πf}}_{0}Mn+{\mathrm{\θ}}_0)}---\left(2\right)$

然后进行多相滤波过程，其中u路未加平移因子的多相滤波的结构图如图2所示。

其中E_p(z)(0≤p≤M-1)为截止频率为π/M的低通滤波器H(z)的线性多相子滤波器，满足

$H\left(z\right)=\underset{p=0}{\overset{M-1}{\Σ}}{z}^{-p}{E}_p\left(z^M\right)---\left(3\right)$

由于各路多相子滤波器E_p(z)的长度相等，故各相子滤波会产生一近似相等的群延时τ，该群延时会引起时延相位但由于各子通道的时延相位近似相等，使得各相子滤波输出信号之间的相对相位，相比于输入前，基本维持不变。从而各相子滤波输出信号x_u,p(n)为

$\begin{array}{c}{x}_{u,p}\left(n\right)\≈a_0\·e^{j\[2{\mathrm{\πf}}_{0}N(n-p)+{\mathrm{\θ}}_0+{\stackrel{\‾}{\mathrm{\θ}}}_0\]}{|}_{n=Mn,f_0=(m+\mathrm{\δ})/\left(MN\right)}\\ =a_0\·e^{j\[2\mathrm{\π}\mathrm{\δ}n-(m+\mathrm{\δ})2\mathrm{\π}p/M+{\mathrm{\θ}}_0+{\stackrel{\‾}{\mathrm{\θ}}}_0\]},0\≤p\≤M-1\end{array}---\left(4\right)$

然后对u路输出的M路子信号按列进行IDFT(即对变量p做IDFT)，得到如下形式的信号：

$\begin{array}{c}{u}_k\left(n\right)=\frac{1}{M}\underset{p=0}{\overset{M-1}{\Σ}}{x}_{u,p}\left(n\right)\·e^{j\frac{2\mathrm{\π}}{M}pk}\\ =e^{j2\mathrm{\π}\mathrm{\δ}n}U\left(k\right)\end{array}---\left(5\right)$

其中，

$U\left(k\right)=e^{j\[({\mathrm{\θ}}_0+{\stackrel{\‾}{\mathrm{\θ}}}_0)-\frac{M-1}{M}(m+\mathrm{\δ}-k)\mathrm{\π}\]}\·\frac{a_{0}sin\[(m+\mathrm{\δ}-k)\mathrm{\π}\]}{M\sin \[(m+\mathrm{\δ}-k)\mathrm{\π}/M\]}---\left(6\right)$

同理，v路经IDFT变换后的输出信号为

$\begin{array}{c}{v}_l\left(n\right)=\frac{1}{N}\underset{q=0}{\overset{N-1}{\Σ}}{x}_{v,q}\left(n\right)\·e^{j\frac{2\mathrm{\π}}{N}qk}\\ =e^{j2\mathrm{\π}\mathrm{\δ}n}V\left(l\right)\end{array}---\left(7\right)$

其中，

$V\left(l\right)=e^{j\[({\mathrm{\θ}}_0+{\stackrel{\‾}{\mathrm{\θ}}}_0)-\frac{N-1}{N}(m+\mathrm{\δ}-l)\mathrm{\π}\]}.\frac{a_{0}sin\[(m+\mathrm{\δ}-l)\mathrm{\π}\]}{Nsin\[(m+\mathrm{\δ}-l)\mathrm{\π}/N\]}---\left(8\right)$

由式(5)～(8)，可得经互相关平均后的频谱

$\begin{array}{c}{S}_{xx}\left(e^{{\mathrm{j\ω}}_i}\right)\≈E\[e^{j2\mathrm{\π}\mathrm{\δ}n}U\left(k\right)\·e^{-j2\mathrm{\π}\mathrm{\δ}n}{V}^{*}\left(l\right)\]\\ =U\left(k\right)\·V^{*}\left(l\right)E\[e^{j2\mathrm{\π}\mathrm{\δ}n}\·e^{-j2\mathrm{\π}\mathrm{\δ}n}\]\\ =U\left(k\right)\·V^{*}\left(l\right)\end{array}---\left(9\right)$

对于u路而言，由式(5)可以看出，由于IDFT是针对x_u,p(n)的变量p而言的，且式(4)中x_u,p(n)的数字角频率(m+δ)2π/M很可能超出2π，故结合DFT的循环周期性，只有当式(5)的谱序号k取值为(m+δ)模除M的四舍五入取整结果时，即满足

k＝[(m+δ)mod M]＝m mod M       (10)

IDFT谱线才出现峰值。

同理，对于v路而言，只有当式(7)的谱序号l取值为(m+δ)模除N的四舍五入取整结果时，即满足

l＝[(m+δ)mod N]＝m mod N        (11)

IDFT谱线才出现峰值。

由式(9)可以看出，最终频谱只与u路谱余数位置k和v路谱余数位置l有关，扫描i(0≤i≤MN-1)，可将信号频谱提取出来。

然而原互素感知过程有一个问题，当信号频域小数偏移δ分别接近0和接近0.5时，会出现截然不同的分析效果。

具体说来，当信号归一化频率f₀＝(m+δ)Δf中的频偏值δ接近于0时，可以得到期望的谱估计输出，即在i＝m处，获得大的幅值，当谱序号i≠m时，的值趋于0。

当频偏值δ接近于0.5时，因滤波器无法实现期望的理想特性，则会在相邻谱位置(如主谱峰i₁＝m处和旁谱i₂＝m+1处)时都会出现大的幅值，即出现频谱泄漏效应，更严重的是，由于互素谱是对u、v两路IDFT各输出做互相关扫描得到，因而这会导致主谱峰i₁＝m处和旁谱i₂＝m+1处所对应的IDFT各输出支路出现交叉项干扰。

具体而言，与主峰i₁＝m对应的u路和v路IDFT输出序号分别为k₁＝m modM，l₁＝m mod N，与旁谱i₂＝m+1对应的u路和v路IDFT输出序号分别为k₂＝(m+1)mod M，l₂＝(m+1)mod N，则当扫描至将u路的第k₁个IDFT输出和v路的第l₂个IDFT输出做互相关时，必然会在另一i₃处产生大的幅值(满足k₁＝i₃mod M，l₂＝i₃mod N)；同理，当扫描至将u路的第k₂个IDFT输出和v路的第l₁个IDFT输出做互相关时，必然会在另一i₄处产生大的幅值(满足k₂＝i₄mod M，l₁＝i₄mod N)，很显然，谱序号i₃、i₄的位置必然远离主峰i₁＝m的位置，故称之为“伪谱峰”，这种因IDFT支路输出交叉产生的大幅值谱必然会大大降低谱分析的可读性，使得谱分析性能变得非常差。

为了说明这个问题，以两个单频信号x₁(n)＝exp[j(2π35.03Δfn)]和x₂(n)＝exp[j(2π35.47Δfn)]为例进行说明，其中，δ₁＝0.03，δ₂＝0.47，取M＝11，N＝7，使用长度为221的Remez滤波器作为原型滤波器。其互素谱分析结果分别如图3(a)、(b)所示。

图3(a)是δ₁＝0.03时的谱估计，可以看出，互素感知较准确的估计出信号的频率35；图3(b)是δ₂＝0.47时的谱估计，由于δ接近0.5，故产生了谱泄漏和伪谱峰，即除了在期望的i₁＝m＝35处，获得大的幅值外，在相邻谱位置(即i₂＝m+1＝36处)，也泄漏出现大幅值旁谱；此外，由于与主峰i₁＝35对应的u路和v路IDFT输出序号分别为k₁＝35mod11＝2，l₁＝35mod7＝0，与旁谱i₂＝36对应的u路和v路IDFT输出序号分别为k₂＝36mod 11＝，3l₂＝36mod 7＝，1故在i₃＝14(因k₂＝14mod11＝3，l₁＝14mod7＝0)和i₄＝57(因k₁＝57mod11＝2，l₂＝57mod7＝1)处出现了伪谱峰。

202：平移通路的互素谱感知信号流程分析；

针对这种情况，本发明对原信号进行0.5倍分辨率平移后再进行一次互素感知处理，使得无论信号频率的小数偏移是否接近0.5，对于各个频率都可以保证至少有一种情况没有谱泄漏。设频移后的原信号为

$\tilde{x}\left(n\right)=x\left(n\right)\·e^{j2\mathrm{\π}0.5\mathrm{\Δ}fn}=a_0\·e^{j\[2\mathrm{\π}(f_0+0.5\mathrm{\Δ}f)n+{\mathrm{\θ}}_0\]}---\left(12\right)$

以u路为例，经过互素感知处理后的多相支路信号为

${\tilde{x}}_{u,p}\left(n\right)={(-1)}^n\·a_0\·e^{j\[2\mathrm{\π}\mathrm{\δ}n-(m+\mathrm{\δ})2\mathrm{\π}p/M+{\mathrm{\θ}}_0+{\mathrm{\θ}}_1\]}\·e^{j\frac{\mathrm{\π}}{M}p}---\left(13\right)$

根据式(13)可知，为了减小计算复杂度，本发明经过等效转换，可将频偏加在多相滤波过程中的低数据率端。

加入平移因子的多相滤波流程如下所示，图4中，相应的平移因子λ_u由下式得出

λ_u＝(-1)ⁿ·e^jpπ/M                          (14)

与之类似，经过互素感知，v路的平移前多相支路信号x_v,q(n)和平移后的多相支路信号为

$x_{v,q}\left(n\right)=a_0\·e^{j\[2\mathrm{\π}\mathrm{\δ}n-(m+\mathrm{\δ})2\mathrm{\π}q/N+{\mathrm{\θ}}_0+{\mathrm{\θ}}_1\]},0\≤q\≤N-1---\left(15\right)$

$\begin{array}{c}{\tilde{x}}_{v,q}\left(n\right)=a_0\·e^{j\[2\mathrm{\π}\mathrm{\δ}n-(m+\mathrm{\δ})2\mathrm{\π}q/N+{\mathrm{\θ}}_0+{\mathrm{\θ}}_1\]}\·{\left(-1\right)}^n\·e^{jq\mathrm{\π}/N}\\ =x_{v,q}\left(n\right)\·{\mathrm{\λ}}_v\end{array}---\left(16\right)$

其中，

λ_v＝(-1)ⁿ·e^jqπ/N                            (17)

经IDFT后求互相关，可得加平移因子后的频谱假设输入信号以及其他参数与图3仿真相同，分别用原互素感知和加平移因子的互素感知对频率进行估计，如下图所示：

图5(a)为δ＝0.03时原通道的频谱图，可以看出没有发生谱泄漏，(c)为对应的平移通道频谱图，由于加了平移，出现了谱泄漏，(b)为δ＝0.47时原通道的频谱图，出现了谱泄漏，(d)为对应的平移通道频谱图，加平移后消除了谱泄漏。由图4可以得出，无论信号频率的小数偏移δ是多少，总是至少能找到一种情况(原直通路或平移通路)没有谱泄漏，这样，根据后续步骤的方法，无论信号频率小数偏移如何，都可以将信号频率准确地识别出来。

203：频率指示器；

通过以上描述以及图5可以看出，对于任意信号成分，由于本发明提出了平移通路的概念，无论它的频率小数偏移δ是多少，总能得到没有谱泄漏的频谱图。基于此，本发明分三种情况进行讨论：

当|δ|→0.5时，平移通道频谱没有谱泄漏，原通道在其他多个位置有谱泄漏，如图5(b)和(d)所示。为了得到较准确的频率估计，需要将真实频率以及与它相邻的泄漏出来的谱线提取出来，作为后续逻辑判决的依据；

当|δ|→0时，原通道没有谱泄漏，平移通道频谱在其他多个位置有谱泄漏，如图5(a)和(c)所示。此时只需要将真实谱线提取出来，作为后续逻辑判决的依据；

当|δ|位于(0,0.5)的中间位置时，原通道和平移通道频谱都没有谱泄漏，需要将真实谱线提取出来，作为后续逻辑判决的依据。

为此，首先将平移通道频谱向左平移一个单位得再构造以下两个乘式ρ₁(i)和ρ₂(i)：

$\{\begin{array}{c}{\mathrm{\ρ}}_1\left(i\right)=S_{xx}\left(e^{{\mathrm{j\ω}}_i}\right)\·{\tilde{S}}_{xx}\left(e^{{\mathrm{j\ω}}_{i+1}}\right)\\ {\mathrm{\ρ}}_2\left(i\right)=S_{xx}\left(e^{{\mathrm{j\ω}}_i}\right)\·{\tilde{S}}_{xx}\left(e^{{\mathrm{j\ω}}_i}\right)\end{array}---\left(18\right)$

表1为δ→0.5时原互素谱和平移互素谱的输出位置大小以及相对应的两个频率指示器的大小。

表1 4个位置的峰值大小(δ→0.5)

204：逻辑判决器。

得到两个频率指示器ρ₁(i)和ρ₂(i)后，分三步进行判定。在此之前，如表1所示，当频率小数偏移δ→0.5时，由于频谱泄漏使得原互素谱在四个位置出现较大谱峰，而又由于频谱泄漏导致能量不均等地散布在这四个位置，使得无法判断具体的数值大小而只能粗略地判断出泄漏出来的谱与真实谱大小相当，而在整个频谱图中还存在由于噪声产生的与泄漏谱相比小得多的干扰谱，因此，判定之前需要设定一个阈值ξ，将噪声干扰谱去除而只留下感兴趣的真实谱和泄漏谱。

首先筛选真实存在的谱峰(包括真实谱峰与干扰出来的伪峰)。设定一个阈值ξ，i＝0扫描ρ₁(i)和ρ₂(i)，选择ρ₁(i)＞ξ或ρ₂(i)＞ξ的谱峰进行下一步判定；

判定频率小数偏移δ是否接近0.5。由表1所示，需要ρ₁(i)＞ξ，ρ₂(i+1)＞ξ，并且ρ₂(i)＜ξ，当δ接近-0.5时也有相同的判断。

判定成功后，将的值作为真实的输出，又因为δ→0.5，因此，真实的频谱偏移应该位于δ＝0.5附近，因此，可将谱峰定位在归一化频率(i+0.5)Δf位置(由于谱序号要求为整数，故其谱序号为j＝(2i+1))，提高了频率分辨率。

对于频率小数偏移δ远离0.5的情况(即ρ₁(i)＞ξ，ρ₂(i+1)＞ξ，ρ₂(i)＜ξ，三个条件中只要有1个不满足)，可将谱峰定位在归一化频率iΔf位置(由于谱序号要求为整数，故其谱序号为j＝2i)，提高了频率分辨率。

若不满足ρ₁(i)＞ξ或ρ₂(i)＞ξ的条件(即这时ρ₁(i)≤ξ且ρ₂(i)≤ξ)，说明当前扫描iΔf和(i+0.5)Δf位置的两个频率位置均没有大的谱成分，故将这两个位置的谱输出均赋为0值。

根据以上三个步骤可去除由于|δ|→0.5造成的伪峰干扰，将给定信号的归一化频率在频谱图中的位置识别出来，并且提高了频率分辨率。需要强调的是，在每次判定后，要将该位置上已用过的大幅值指示器ρ₁(i)、ρ₂(i)或ρ₂(i+1)清零，以免影响下一个位置的判定。逻辑判决器算法流程如图6所示。

实施例3

(1)单频谱泄漏抑制实验

实验1

设定两个互素的整数M＝5，N＝4，单频信号x(t)＝1.2exp(j2π·10.51t)，δ＝-0.49，理想采样率F_s＝MN＝20Hz(相应的采样间隔T＝1/F_s＝1/MN)，频率分辨率Δf＝F_s/MN＝1Hz，则根据互素谱分析方法，两路实际欠采样值分别为f_s1＝1/NT＝M＝5Hz，f_s2＝1/MT＝N＝4Hz，均远低于信号频率10.51Hz。用长度为121的Remez方法设计两路原型低通滤波器，图7(a)～7(e)给出了原互素谱平移互素谱指示器ρ₁(i)、指示器ρ₂(i)及其本发明得到的改进的分辨率提高1倍的互素谱Y(j)。

图7(a)为原互素谱估计器的频谱输出，由于δ＝-0.49，在4个位置出现谱线，分别位于i＝6，10，11，15，而事实上只有i＝10或11才是真实的谱线，其它的为伪谱峰；(b)为0.5平移的互素谱估计器的频谱输出，只有在i＝1处出现一根谱线；(c)和(d)为频谱指示器，分别提取出i＝10和i＝11处的真实谱线；(d)为经逻辑判决器后的输出频率jΔω/2＝10.5Δω，可见本发明提出的频谱估计器不仅能将真实谱线提取出来，消除伪峰效应和谱泄漏效应，而且其频率分辨率相对于原始互素谱也提高了1倍。

(2)多频谱泄漏抑制实验

实验2

本实验采用与实验1相同的参数和信号模型，只是输入信号由单频变为多频，三个频率分别为f₁＝7.51Hz，f₂＝15.02Hz，f₃＝13.23Hz，图8为本发明提出的频谱估计器的估计过程。

由图8可以看出，由于三个频率小数偏移δ分别为-0.49，0.02和0.23，致使无论是原互素谱输出还是0.5平移互素谱输出都在很大频率范围内出现了伪峰，但两个频率指示器只在真实的谱峰位置附近出现谱线(ρ₁(7)、ρ₁(13)、ρ₁(15)和ρ₂(8))，根据逻辑判决，可将正确的频谱提取出来，分别为7.5Δω，13Δω和15Δω(对应j＝15，26和30)。

实施例4

一种消除伪峰、谱泄漏效应的互素谱分析装置，参见图9，该装置包括：数字信号处理器，及输出驱动及显示电路，

数字信号处理器，用于对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置；

输出驱动及显示电路用于输出归一化频率在频谱图中的位置。

其中，数字信号处理器包括：第一获取模块，

其中，第一获取模块，用于对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

其中，数字信号处理器还包括：第二获取模块、输出模块和第三获取模块，

其中，第二获取模块，用于对输入信号进行两路下采样，得到两路稀疏信号；

其中，输出模块，用于对两路稀疏信号分别进行直接多相滤波，以及分别加平移因子后的多相滤波，输出两对多通道稀疏信号；

其中，第三获取模块，用于在每个时刻，对多通道输出序列分别进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

实施例5

在图9中，首先将实际采集的信号以及滤波器系数存入外部RAM中，再将它们实时输入到DSP中，经过DSP内部核心算法，对信号进行下采样、多相滤波、IDFT、频率指示以及构建逻辑判决器等，对信号频率进行估计，最后借助输出驱动显示及其显示模块显示出频率值。

其中，图9的DSP(Digital Signal Processor，数字信号处理器)为核心器件，在频率重构的过程中，完成如下主要功能：

调用内部核心算法，完成实际采集信号的下采样、多相滤波、IDFT、频率指示和逻辑判决等过程；

控制下采样值以及信号样本，实时对其进行调整，使其符合实际需要；

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

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

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

本发明将所提出的核心算法植入DSP器件内，基于此完成高精度、低复杂度、高效的频率估计。

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

首先根据实际需要，设置信号的下采样值(M和N，互素的整数)，设计所需的滤波器系数，用于多相滤波过程，并设定逻辑判决所需的阈值；

然后，CPU主控器从I/O端口读取参设定的参数，进入内部RAM；

按图1本发明的处理过程进行频率估计器的设计是DSP算法最核心的部分，运行该算法后，即可得到所采集的信号的频率估计；

判断本发明方法是否满足实际需求，若不满足，程序返回，重新根据要求设定信号参数；

直至设计结果符合实际要求，然后通过DSP的输出总线输出至外部显示驱动设备，将频率估计结果进行数码显示。

需指出，由于采用了DSP实现，使得整个频谱估计器设计变得更为灵活快捷，可根据频谱估计器设计过程中的实际需要，灵活变换所需参数，使之最终符合工程需要。

参考文献

[1]青海银,张援农,周晨,等.基于MST雷达垂直风速的大气温度剖面反演[J].物理学报,2014,63(9):94301-094301.

[2]Ding Z,Yao X S,Liu T,et al.Long-range vibration sensor based on correlation analysis ofoptical frequency-domain reflectometry signals[J].Optics express,2012,20(27):28319-28329.

[3]黄翔东,孟天伟,丁道贤,等.前后向子分段相位差频率估计法[J].物理学报,2014,63(21):214304-214304.

[4]丛超,李秀坤,宋扬.一种基于新型间歇混沌振子的舰船线谱检测方法[J].物理学报,2014,63(6):64301-064301.

[5]Wu Q,Liang Q.Coprime sampling for nonstationary signal in radar signal processing[J].EURASIP Journal on Wireless Communications and Networking,2013,2013(1):1-11.

[6]Zhang Y D,Qin S,Amin M G.DOA estimation exploiting coprime arrays with sparse sensorspacing[J].Proc.IEEE ICASSP,Florence,Italy,2014:2267-2271.

[7]Wang W,Xia X G.A closed-form robust Chinese remainder theorem and its performanceanalysis[J].Signal Processing,IEEE Transactions on,2010,58(11):5655-5666.

[8]Li X,Liang H,Xia X G.A robust Chinese remainder theorem with its applications infrequency estimation from undersampled waveforms[J].Signal Processing,IEEE Transactionson,2009,57(11):4314-4322.

[9]黄翔东,丁道贤,南楠,等.基于中国余数定理的欠采样下余弦信号的频率估计[J].物理v学报,2014,63(19):198403-198403.

[10]Pal P,Vaidyanathan P P.Nested arrays:a novel approach to array processing with enhanceddegrees of freedom[J].Signal Processing,IEEE Transactions on,2010,58(8):4167-4181.

[11]Vaidyanathan P P,Pal P.Sparse sensing with co-pprime samplers and arrays[J].IEEETransactions on Signal Processing,2011,59(2):573-586.

[12]Vaidyanathan P P,Pal P.Theory of sparse coprime sensing in multiple dimensions[J].SignalProcessing,IEEE Transactions on,2011,59(8):3592-3608.

[13]Vaidyanathan P P,Pal P.Sparse co-prime sensing with multidimensional latticearrays[C]//Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop(DSP/SPE),2011IEEE.IEEE,2011:425-430.

[14]Vaidyanathan P P,Pal P.Coprime Sampling and Arrays in One and Multiple Dimensions[M]//Multiscale Signal Analysis and Modeling.Springer New York,2013:105-137.

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

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

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

Claims (6)

1.一种消除伪峰、谱泄漏效应的互素谱分析方法，其特征在于，所述互素谱分析方法包括以下步骤：

对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；

平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；

将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置。

2.根据权利要求1所述的一种消除伪峰、谱泄漏效应的互素谱分析方法，其特征在于，所述对离散傅里叶变换后的输出序列进行互相关扫描求值的步骤之前，所述方法还包括：

对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

3.根据权利要求1所述的一种消除伪峰、谱泄漏效应的互素谱分析方法，其特征在于，所述对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列的步骤具体为：

1）对输入信号进行两路下采样，得到两路稀疏信号；

2）对两路稀疏信号分别进行直接多相滤波，以及分别加平移因子后的多相滤波，输出两对多通道稀疏信号；

3）在每个时刻，对多通道输出序列分别进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

4.一种消除伪峰、谱泄漏效应的互素谱分析装置，所述装置包括：数字信号处理器，及输出驱动及显示电路，其特征在于，

所述数字信号处理器用于对离散傅里叶变换后的输出序列进行互相关扫描求值，得到原频谱序列和平移后的频谱序列；平移后的频谱序列向左频移一个单位，原频谱序列与频移后的序列的对应元素相乘，获取第一频率指示序列；将原频谱序列和平移后的频谱序列的对应元素相乘，获取第二频率指示序列；将第一频率指示序列和第二频率指示序列输入逻辑判决器，根据两个频率指示序列中同一频率的对应关系，获取归一化频率在频谱图中的位置；

所述输出驱动及显示电路用于输出归一化频率在频谱图中的位置。

5.根据权利要求4所述的一种消除伪峰、谱泄漏效应的互素谱分析装置，其特征在于，所述数字信号处理器包括：第一获取模块，

所述第一获取模块，用于对由多通道稀疏信号组成的多通道输出序列进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。

6.根据权利要求4所述的一种消除伪峰、谱泄漏效应的互素谱分析装置，其特征在于，所述数字信号处理器还包括：第二获取模块、输出模块和第三获取模块，

所述第二获取模块，用于对输入信号进行两路下采样，得到两路稀疏信号；

所述输出模块，用于对两路稀疏信号分别进行直接多相滤波，以及分别加平移因子后的多相滤波，输出两对多通道稀疏信号；

所述第三获取模块，用于在每个时刻，对多通道输出序列分别进行离散傅里叶变换，得到离散傅里叶变换后的输出序列。