CN101231315A - 频率估计的多段采样信号融合处理方法 - Google Patents

Abstract

多段采样信号融合处理是一种适用于信号频率估计的频域分析方法，通过为各段采样信号设定不同的频域分析参数和相参补偿因子，本方法能补偿由于各段被测信号之间的频率不等、相位不连续对信号处理结果的影响，实现相关信息融合和有用信息提取，以获得更高的频率估计精度。本方法的应用范围广泛，适用于各段采样信号的初相未知、相邻采样段之间的时间间隔未知、在相邻采样段之间的时间间隔内被采样信号发生未知变化、各次采样过程的采样频率不同等情况，能够实现对多段不同频信号的高精度处理，满足多段同频信号、多段降频信号、多段分频信号等情况下的频率估计需要。本方法给出了理论推导过程和验证实验结果与对比，进行了实时计算量分析。

Description

频率估计的多段采样信号融合处理方法

技术领域

本发明涉及信号频率估计，具体的说，涉及多段采样信号的处理方法。

背景技术

多段采样信号指的是，对同一被测对象，在不同时间段或不同采样参数下，通过数次采样过程，得到了数段采样信号。各次被采样对象在该次采样过程中的频率不变，各段采样所得信号(以下简称为“采样信号”)的频率之间的相对关系已知，相邻段采样之间的时间间隔可以未知，在相邻次采样之间的时间间隔内被采样信号可以发生未知变化，各次采样过程的采样频率可以不同。

多段采样信号尤其是多段同频信号处理问题十分常见，例如，在一些测量设备的校准过程中，需要对已知对象进行多次测量，然后综合评估设备的工作状况；还例如，在LFMCW雷达测距过程中，多个回波信号经过延迟-混频、低通滤波后，表现为同一单频信号的多次采样；在电子对抗过程中，对相雷达进行数字侦收时，也将获得多段(多个脉冲)采样信号。

众所周知，相比于单段信号，多段同频信号应含有更多的有用信息。但一直以来对多段同频信号进行深度信息融合式的综合处理技术不受重视。其原因在于：一方面，该融合处理过程会产生较大的运算量，当时运算处理器的性能相对不足；另一方面，以往的研究中，多将待处理信号近似为某一段时间内的平稳信号且背景噪声的幅度相对较低，从而得到采样持续时间较长的平稳信号，即单段采样数据中包含有足量的、易于提取的信息，已能满足当时的测量要求。在此情况下，虽然多段采样信号十分常见，但其融合处理技术没有得到应有的发展。

近年来，非平稳信号、低信噪比信号的研究日益受到重视，由于被测频率保持相对稳定的时间非常短暂，相当于用一个持续时间很短的窗函数去乘该信号，且采样结果中包含大量噪声，这给信号处理带来了很大的困难。例如，为达到测量精度要求，需要一段信噪比为5dB、持续时间超过0.1毫秒的平稳信号，而实际采样过程只能得到4段信噪比为-5dB、持续时间仅为0.05毫秒的平稳信号。为此，除了进一步改进现有的单段采样信号处理方法和提高硬件性能外，另一种重要的途径就是尽可能利用多段采样信号中包含的信息。近几年来，多段采样信号处理技术逐渐得到关注，已有研究人员开展了初步研究，目前已有的处理方法可分为三种。

(1)直接累加法：对各段采样信号分别进行如DTFT等频域分析，然后将分析结果累加。该方法计算量小，能够应用于各种情况，但是精度较低。

(2)相位关联法：一种对多个相位相参的分段信号进行频谱估计的方法，利用各段信号之间的相位关系，通过相位联系对由单段信号得到的频谱估计做进一步的修正。该方法对多段采样信号提出了较多的限制条件，且推导过程中采用了近似处理较粗略，难以实际应用。

(3)相位积累法：包括直接相位积累法和旋转相位积累法，将分段采样信号进行相位关联，以实现相参积累和提高频率估计精度。直接相位积累法运算量大，并要求各段采样之间的空闲间隔时间已知且在该空闲间隔内被测信号保持稳定不变，应用范围较小。旋转相位积累法是目前较先进的一种方法，但是受噪声影响大，且要求各段采样信号的初相已知，严重制约了该方法的应用与推广。

综上所述，对多段采样信号处理方法的研究具有重要意义，本专利要解决的技术问题是提出一种高精度方法，具有普遍适用性、能够处理各种类型的多段采样信号，以较小的实时计算量获得较高的精度。

发明的内容

由于本专利方法是线性的，具有可加性，因此分析的结果同样适用于多段多频信号；由于本专利方法能够应用于多段不同频信号，因此分析的结果同样适用于多段不同采样频率情况；为便于叙述，本文均以被测信号为单频信号、采样频率相同为例。

1.三类多段采样信号

本专利方法将现有的多段采样信号问题划分为三类：

(1)在噪声和干扰背景下，对被测信号进行重复采样，被测信号的频率和各段采样信号的初相未知。

(2)在噪声和干扰背景下，对被测信号进行多次采样得到多段采样信号，各次采样所得信号的频率、初相未知且频率并不相等，但各次采样所得信号频率之间的频率差已知。

(3)在噪声和干扰背景下，对被测信号进行多次采样得到多段采样信号，各次采样所得信号的频率、初相未知且频率并不相等，但各次采样所得信号频率之间的倍数(可为任意小数)已知。

将以上三类情况依次简称为多段同频信号、多段降频信号和多段分频信号，其中后两类可统称为多段变频信号处理。可以看出，本专利方法将多段采样信号处理研究从单纯的重复测量即多段同频信号处理拓展到非同频信号处理。这是因为：在多段变频信号中，虽然各段信号的频率各不相同，但是由于它们的频率之间的相对关系已知，所以本质上仍然可以将它们视作同频信号的重复测量。这种新型的重复测量关系是本专利方法的一个特色，将本专利方法的应用领域从平稳信号扩展到非平稳信号，扩展了本专利方法的应用范围。

2.通用处理流程

由于本专利方法对以上三类多段采样信号进行频率估计的构思基本相同，为便于说明本专利方法的核心内容，特归纳出本专利方法的通用处理流程如附图所示，并给出具体说明如下(为便于叙述，以下数学表达式，若未经特别说明，则均参照MATLAB 7函数规定解释)：

设有M段采样信号，首段信号频率为f₁，各段信号频率间的相对关系已知，采样频率为f_s，各段初相为θ(m)，各段长度为D(m)，前(m-1)段的长度和为B(m)。

(1)根据首段采样信号频率的估计值f₂和最大相对估计误差C确定首段待比较频率序列BF0(k)，k∈[1，nf₂]

(2)根据各段采样信号之间的相对关系，得出各段采样信号的待比较频率序列BF(m，k)。

设Df(m)为各段信号降频幅度，p(m)为各段信号分频比，则：

BF(m，k)＝[BF0(k)-Df(m)]·p(m)。

(3)将各段信号进行BF(m，k)频率处的频域分析，得到Z₀(m，k)。

(4)根据Z₀(m，k)利用直接累加法得到首段采样信号频率的直接累加法估计值f₁₀。

(5)设 $\begin{array}{c}\mathrm{FS}(f_1\left(t\right),h,k)=\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{\mathrm{pD}}\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\\ \·\frac{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{\mathrm{pD}}\·M\}}{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·\stackrel{\‾}{\mathrm{pD}}/f_s\}}\end{array}$ ，其中：

f₁(t)∈[BF0(1)，BF0(nf₂)]，h，k∈[1，nf₂]，pD＝mean(pD)，pD＝p·D。

调出预存的的极大值表S(n，f₁(t)，h)，该表为函数FS(f₁(t)，h，k)的n₁个极大值序列，n∈[1，n₁]。设频率估计区间GF＝[gf₁，gf₂]＝[0.997·f₁₀，1.003·f₁₀]，取该函数表中f₁(t)∈GF部分，得到S₂(n，f₁(t)，h)。

(6)计算 $S_3(n,h)=\underset{i}{\overset{\mathrm{fn}}{\mathrm{\Σ}}}\left|S_2\right[n,f_1(t+1),h]-S_2[n,f_1\left(t\right),h\left]\right|.$ 取S₃(n，h)中较大值对应的h序列h₀(l)，取S₂(n，f₁(t)，h)中h∈h₀(l)部分，得到S₄(n，f₁(t)，h₀(l))。因为一般情况下，n₁＝1即可满足精度要求，为便于说明，取n₁＝1，即S₄(n，f₁(t)，h₀(l))＝S₄(f₁(t)，h₀(l))。

(7)设h₀的长度为nf₁，确定搜索频率序列SF0(l)＝BF0(h₀(l))，l∈[1，nf₁]

(8)部分情况下可以直接根据经验给出θ(m)的估计值θ₁(m)。如果发现该经验值的误差过大，例如θ₁(m)-θ(m)＞pi，其中θ₁(m)，θ(m)∈[0，2*pi]，则可以通过近似计算得到对应于不同信号段的初相估计值θ₂(m)，方法如下：

找到BF(m，k)中最接近f₁₀的频率点BF(m，k₀(l))。

得到θ₂(m)＝angle[Z₀(m，k₀(l))]，其中angle(t)表示复数t的角度。

在接下来的说明中以较复杂的θ₂(m)为例，若可以采用θ₁(m)计算，则只需作相应简化即可，不再赘述。

(9)根据θ₂(m)、BF0(k)和SF0(l)确定相参补偿因子矩阵φ(l，m，k)。

其中，针对多段同频信号和多段降频信号，有：

φ(l，m，k)＝θ₂(m)-θ₂(1)-2·pi·B(m)·[SF0(l)-BF0(k)]/f_s；

针对多段分频信号，有：

$\left\{\begin{array}{c}\mathrm{\φ}(l,m,k)={\mathrm{\θ}}_2\left(m\right)-{\mathrm{\θ}}_2\left(1\right)-2\·\mathrm{pi}\·\mathrm{pB}\left(m\right)\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\\ \mathrm{pB}\left(m\right)=\mathrm{sum}\left(\mathrm{pD}(1:m-1)\right)\end{array}\right.$

(10)计算如下：

$\left\{\begin{array}{c}{Z}_1(l,m,k)=\exp [-j\·\mathrm{\φ}(l,m,k)]\·Z_0(m,k)\\ Z_2(l,k)=\left|\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{Z}_1(l,m,k)\right|\end{array}\right.$

(11)取Z₂(l，k)中的谱峰频率位置Z₃(l)＝BF(MAX[Z₂(l，k)])。

(12)根据Z₃和S₄，计算 $\left\{\begin{array}{c}{Z}_4(f_1\left(t\right),l)=|Z_3\left(l\right)-S_4(f_1\left(t\right),h_0\left(l\right))|\\ Z_5\left(f_1\left(t\right)\right)=\underset{l}{\mathrm{\Σ}}{Z}_4(f_1\left(t\right),l)\end{array}\right.,$ 取Z₅的最小值对应f₁(t₀)，得到f₃＝f₁(t₀)；

(13)若f₃位于频率区间GF的左端处即f₃-gf₁≤Δ(Δ根据实际情况取值)，则左移频率区间GF和相应修改S₄，然后重复上述流程中的第(12)步。若f₃位于频率区间GF的右端处gf₂-f₃≤Δ(Δ根据实际情况取值)，则右移频率区间GF和相应修改S₄，然后重复上述流程中的第(12)步。

(14)重复上述流程中的第(13)步，直至f₃不再位于频率区间GF的左右端，该频率f₃即为本专利方法的频率估计结果。

3.理论证明

为更好说明本专利方法的优越性，特给出理论证明过程如下：

(1)针对多段同频信号的理论证明

凡是对同一平稳信号进行多次测量的场合都会产生多段同频信号，其一般形式如下：设有M段同频信号，信号频率为f₁，采样频率为f_s，各段长度为D(m)，各段初相为θ(m)，即：

$\left\{\begin{array}{c}B\left(1\right)=0,a\∈[1,D\left(m\right)],m\∈[1,M]\\ B\left(s\right)=\underset{m=1}{\overset{s-1}{\mathrm{\Σ}}}D\left(m\right),s\∈[2,M+1]\\ x_0[B\left(m\right)+a]=\cos [\mathrm{\θ}\left(m\right)+2\·\mathrm{pi}\·f_1\·(a-1)/f_s]\end{array}\right.$

根据上述通用处理流程，可知：

$\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\exp [-j\·\mathrm{\φ}(l,m,k)]\·\underset{a=1}{\overset{D\left(m\right)}{\mathrm{\Σ}}}\exp \{j\·[\mathrm{\θ}\left(m\right)+2\·\mathrm{pi}\·f_1\·a/f_s\left]\right\}\·\exp \{-j\·2\·\mathrm{pi}\·\mathrm{BF}0\left(k\right)\·a/f_s\}$

$\≈\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·D\left(m\right)\}}{2*\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}$

$\·\exp [j\·\{\mathrm{\θ}\left(m\right)-\mathrm{\φ}(l,m,k)+\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]\·[D\left(m\right)+1]|/f_s\}]$

$\≈\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·D\left(m\right)\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\·\exp [j\·\{\mathrm{\θ}\left(1\right)+$

$\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·[2\·B\left(m\right)+D\left(m\right)+1]/f_s+\mathrm{pi}\·[f_1-\mathrm{SF}0\left(l\right)]\·[D\left(m\right)+1]/f_s\left\}\right]$

$\≈\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\·\frac{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\·M\}}{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·\stackrel{\‾}{D}/f_s\}}$

$\·\exp [j\·\{\mathrm{\θ}\left(1\right)+[B(M+1)+1]\·\mathrm{pi}/[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s+(\stackrel{\‾}{D}+1)\·\mathrm{pi}\·[f_1-\mathrm{SF}0\left(l\right)]/f_s\left\}\right]$

证毕。

(2)针对多段降频信号的理论证明

产生多段降频信号的情况也比较多，例如LFMCW雷达测距过程中，对同一目标使用不同延时长度的延迟线组合则可以得到多段降频信号，其一般形式如下：设有M段降频信号，首段信号频率为f₁，各段信号降频幅度为Df(m)，采样频率为f_s，各段长度为D(m)，各段初相为θ(m)，即：

$\left\{\begin{array}{c}B\left(1\right)=0,\mathrm{Df}\left(1\right)=0,a\∈[1,D\left(m\right)],m\∈[1,M]\\ B\left(s\right)=\underset{m=1}{\overset{s-1}{\mathrm{\Σ}}}D\left(m\right),s\∈[2,M+1]\\ x_0[B\left(m\right)+a]=\cos \{\mathrm{\θ}\left(m\right)+2\·\mathrm{pi}\·{[f}_1-\mathrm{Df}\left(m\right)]\·(a-1)/f_s\}\end{array}\right.$

根据上述通用处理流程，可知：

$\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\exp [-j\·\mathrm{\φ}(l,m,k)]\·\underset{a=1}{\overset{D\left(m\right)}{\mathrm{\Σ}}}\begin{array}{c}\exp [j\·\{\mathrm{\θ}\left(m\right)+2\·\mathrm{pi}\·[f_1-\mathrm{Df}\left(m\right)]\·a/f_s\left\}\right]\\ \·\exp \{-j\·2\·\mathrm{pi}\·[\mathrm{BF}0\left(k\right)-\mathrm{Df}\left(m\right)]\·a/f_s\}\end{array}$

$\≈\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·D\left(m\right)\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}$

$\·\exp [j\·\{\mathrm{\θ}\left(m\right)-\mathrm{\φ}(l,m,k)+\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]\·[D\left(m\right)+1]|/f_s\}]$

$\≈\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·D\left(m\right)\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\·\exp [j\·\{\mathrm{\θ}\left(1\right)+\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]$

$\·[2\·B\left(m\right)+D\left(m\right)+1]/f_s+\mathrm{pi}\·[f_1-\mathrm{SF}0\left(l\right)]\·[D\left(m\right)+1]/f_s\left\}\right]$

$\≈\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\·\frac{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\·M\}}{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·\stackrel{\‾}{D}/f_s\}}$

$\·\exp [j\·\{\mathrm{\θ}\left(1\right)+[B(M+1)+1]\·\mathrm{pi}/[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s+(\stackrel{\‾}{D}+1)\·\mathrm{pi}\·[f_1-\mathrm{SF}0\left(l\right)]/f_s\left\}\right]$

证毕。

(3)针对多段分频信号的理论证明

产生多段分频信号的情况也比较多，例如LFMCW雷达测距过程中，对同一目标使用不同电压调频斜率的压控振荡器则可以得到多段分频信号，其一般形式如下：设有M段分频信号，各段信号的频率为f(m)＝p(m)·f₁，p(1)＝1，p(m+1)≥p(m)，采样频率为f_s，各段长度为D(m)，各段初相为θ(m)，即：

$\left\{\begin{array}{c}B\left(1\right)=0,a\∈[1,D\left(m\right)],m\∈[1,M]\\ B\left(s\right)=\underset{m=1}{\overset{s-1}{\mathrm{\Σ}}}D\left(m\right),s\∈[2,M+1]\\ x_0[B\left(m\right)+a]=\cos [\mathrm{\θ}\left(m\right)+2\·\mathrm{pi}\·{f\left(m\right)\·f}_1\·(a-1)/f_s]\end{array}\right.$

根据上述通用处理流程，可知：

$\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\exp [-j\·\mathrm{\φ}(l,m,k)]\·\underset{a=1}{\overset{D\left(m\right)}{\mathrm{\Σ}}}\begin{array}{c}\exp [j\·\{\mathrm{\θ}\left(m\right)+2\·\mathrm{pi}\·f_1\·p\left(m\right)\·a/f_s\left]\right\}\\ \·\exp [-j\·2\·\mathrm{pi}\·\mathrm{BF}0\left(k\right)\·p\left(m\right)\·a/f_s]\end{array}$

$\≈\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·p\left(m\right)\·D\left(m\right)\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}$

$\·\exp [j\·\{\mathrm{\θ}\left(m\right)-\mathrm{\φ}(l,m,k)+\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]\·p\left(m\right)\·[D\left(m\right)+1]|/f_s\}]$

$\≈\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\mathrm{pD}\left(m\right)\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\·\exp [j\·\{\mathrm{\θ}\left(1\right)+\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]$

$\·[2\·\mathrm{pB}\left(m\right)+\mathrm{pD}\left(m\right)+1]/f_s+\mathrm{pi}\·[f_1-\mathrm{SF}0\left(l\right)]\·[\mathrm{pD}\left(m\right)+1]/f_s\left\}\right]$

$\≈\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{\mathrm{pD}}\}}{2\·\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\·\frac{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{\mathrm{pD}}\·M\}}{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·\stackrel{\‾}{\mathrm{pD}}/f_s\}}$

$\·\exp [j\·\{\mathrm{\θ}\left(1\right)+[\mathrm{pB}(M+1)+1]\·\mathrm{pi}/[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s+(\stackrel{\‾}{\mathrm{pD}}+1)\·\mathrm{pi}\·[f_1-\mathrm{SF}0\left(l\right)]/f_s\left\}\right]$

证毕。

4.实验证明

本专利方法给出了理论证明，而且通过大量实验对本专利方法的精度和适应性进行了验证，这一点在后文“具体实施方式”中的两个应用实例的实验结果分析中得到了体现。

附图说明

图本专利方法的具体处理流程

具体实施方式

为进一步说明本专利方法的处理流程，特给出如下应用实例。

1.多段同频信号处理的一个实例及实验结果分析

设有四段同频信号，已知各段采样信号的频率相同，首段采样信号的频率的初始估计值f₂＝3.15MHz，其最大相对估计误差为0.1，即首段采样信号频率f₁∈[2.863636MHz，3.5MHz]。已知采样信号的频率为8MHz，各段采样信号的长度依次为[50 44 48 44]。为精确测量首段采样信号频率，根据本专利方法，可采用如下处理过程：

(1)确定首段待比较频率序列BF0＝linspace(2.863636MHz，3.5MHz，207)。设定该序列长度为207以减小后续的Chirp_Z变换的运算量。

(2)确定各段采样信号的待比较频率序列BF(m，k)＝BF0(k)。

(3)将各段同频信号进行BF0(k)频率处的Chirp-Z变换，得到Z₀(m，k)。

(4)根据Z₀(m，k)利用求和取平均的方法得到首段采样信号频率的直接累加法估计值f₁₀。计算过程如下：

for k＝1:207

z1(k)＝abs(sum(Z₀(1:M，k)))；

end

[z21，z22t]＝max(z1)；f₁₀＝BF0(z22t)；

(5)设 $\begin{array}{c}\mathrm{FS}(f_1\left(t\right),h,k)=\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\}}{2*\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\\ \·\frac{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\·M\}}{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·\stackrel{\‾}{D}/f_s\}}\end{array},$ 其是D=47，

f₁(t)∈[2.863636MHz，3.5MHz]，h，k∈[1，207]。调出预存的的极大值表S(f₁(t)，h)，该表为函数FS(f₁(t)，h，k)的最大值序列。取该函数表中f₁(t)∈(0.997·f₁₀，1.003·f₁₀)部分，得到S₂(f₁(t)，h)。

(6)计算 $S_3\left(h\right)=\underset{t}{\mathrm{\Σ}}\left|S_2\right[f_1(t+1),h]-S_2[f_1\left(t\right),h\left]\right|.$ 取S₃(h)中最大的10个值所对应的h序列h₀(l)，取S₂(f₁(t)，h)中h∈h₀(l)部分，得到S₄(f₁(t)，h₀(l))。

(7)确定搜索频率序列。SF0(l)＝BF0(h₀(l))。

(8)计算θ₂(m)＝angle[Z₀(m，z22t]。

(9)确定 $\left\{\begin{array}{c}\mathrm{\φ}(l,m,k)={\mathrm{\θ}}_2\left(m\right)-{\mathrm{\θ}}_2\left(1\right)-2\·\mathrm{pi}\·B\left(m\right)\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\\ B\left(m\right)=\mathrm{sum}\left(D(1:m-1)\right)\end{array}\right.$

(10)计算如下：

$\left\{\begin{array}{c}{Z}_1(l,m,k)=\exp [-j\·\mathrm{\φ}(l,m,k)]\·Z_0(m,k)\\ Z_2(l,k)=\left|\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{Z}_1(l,m,k)\right|\end{array}\right.$

(11)取Z₂(l，k)中的谱峰频率位置Z₃(l)＝BF0(MAX[Z₂(l，k)])。

(12)根据Z₃和S₄计算得到频率f₃。 $\left\{\begin{array}{c}{Z}_4(f_1\left(t\right),l)=|Z_3\left(l\right)-S_4(f_1\left(t\right),h_0\left(l\right))|\\ Z_5\left(f_1\left(t\right)\right)=\underset{l}{\overset{{\mathrm{nf}}_1}{\mathrm{\Σ}}}{Z}_4(f_1\left(t\right),l)\end{array}\right.,$ 取Z₅的最小值对应f₁(t₀)，得到f₃＝f₁(t₀)；

(13)若f₃位于频率区间GF的左端处即f₃-gf₁≤0.01MHz，则左移频率区间GF和相应修改S₄，然后重复上述流程中的第(12)步。若f₃位于频率区间GF的右端处gf₂-f₃≤0.01MHz，则右移频率区间GF和相应修改S₄，然后重复上述流程中的第(12)步。

(14)重复上述流程中的第(13)步，直至f₃不再位于频率区间GF的左右端，该频率f₃即为本专利方法的频率估计结果。

接下来，结合该例，对本专利方法和目前适应性最好的直接累加法、处理精度最高的旋转相位积累法(为满足旋转相位积累法的特殊应用条件，专为该方法提供其所需特别信息，其他两种方法不使用这些特别信息)进行了全面对比，给出了定量结果。仿真实验中使用的参数设置如下：信噪比SNR＝-5dB；噪声类型为高斯白噪声；各段采样信号的初相均满足幅度为2·pi的高斯分布。

(1)抗噪性仿真

该仿真过程包括500次实验，每次实验的参数除随机函数和信噪比外完全相同，每进行50次实验后信噪比改变一次，得到实验结果如表1所示。表1中的信噪比SNR单位为dB，标准差单位为KHz。可以看出，在多段同频信号、相同信噪比情况下，本专利方法的精度大约提高到其他方法的4～5倍。

表1多段同频信号情况下三种方法的抗噪性能对比

(2)初始估计值f₂的相对误差对处理结果的影响仿真

该仿真过程包括500次实验，每次实验的参数除随机函数和初始估计值外完全相同，每进行50次实验后初始估计值改变一次，得到实验结果如表2所示。表2中的初始估计值f₂的相对误差C＝(f₂-f₁)/f₁，其单位为千分之一，标准差单位为KHz。可以看出，在多段同频信号、相同初始估计值情况下，本专利方法的精度大约提高到其他方法的3～4倍。

表2多段同频信号情况下初始估计误差对三种方法的精度影响对比

(3)对不同首段采样信号频率f₁的适应性仿真

该仿真过程包括500次实验，每次实验的参数除随机函数和首段采样信号频率外完全相同，每进行50次实验后首段采样信号频率改变一次，得到实验结果如表3所示。表3中的首段采样信号频率f₁的单位为MHz，标准差单位为KHz。可以看出，在多段同频信号、相同首段采样信号频率情况下，本专利方法的精度大约提高到其他方法的4～5倍。

表3多段同频信号情况下首段采样信号频率对三种方法的精度影响对比

(4)对不同采样频率f_s的适应性

该仿真过程包括500次实验，每次实验的参数除随机函数和采样频率外完全相同，每进行50次实验后采样频率改变一次，得到实验结果如表4所示。表4中的采样频率f_s的单位为MHz，标准差单位为KHz。可以看出，在多段同频信号、相同采样频率情况下，本专利方法的精度大约提高到其他方法的4～5倍。

表4多段同频信号情况下采样频率对三种方法的精度影响对比

(5)对不同单段信号长度平均值D的适应性

该仿真过程包括500次实验，每次实验的参数除随机函数和单段信号长度平均值外完全相同，每进行50次实验后单段信号长度平均值改变一次，得到实验结果如表5所示。表5中的单段信号长度平均值D的单位为点，标准差单位为KHz。可以看出，在多段同频信号、相同单段信号长度平均值情况下，本专利方法的精度大约提高到其他方法的4～5倍。

表5多段同频信号情况下单段信号长度平均值对三种方法的精度影响对比

(6)对不同信号总段数的适应性

该仿真过程包括500次实验，每次实验的参数除随机函数和信号总段数外完全相同，每进行50次实验后信号总段数改变一次，得到实验结果如表5所示。表5中的信号总段数M的单位为段，标准差单位为KHz。可以看出，在多段同频信号、相同信号总段数情况下，本专利方法的精度大约提高到其他方法的4～5倍。

表6多段同频信号情况下信号总段数对三种方法的精度影响对比

(7)实时计算量对比

直接累加法实时计算量：4·M·(L·log₂(L)+L+1)，其中L≥nf₂+max(D(m))。

旋转相位积累法实时计算量：4·M·D·nf₂

本专利方法实时计算量：

4·M·(L·log₂(L)+L+1)+4·(M+nf₁·M·nf₂)+n·M·nf₁·0.02/0.0005，其中，nf₁一般取30即可；n表示重复上述流程中第(13)步骤的次数，一般取2即可；0.02表示上述流程中第(13)步骤中的频率搜索区域的长度；0.0005表示频率搜索中的的频率步进幅度。

在本例多段同频信号的一般处理情况下，三种方法的实时计算量依次为36880、165600和184256，本专利方法的实时计算量增幅很小。

2.多段分频信号处理的一个实例及实验结果分析

设有四段分频信号，已知首段采样信号的频率的初始估计值f₂＝3.15MHz，其最大估计误差为0.1，即首段采样信号频率f₁∈[2.863636MHz，3.5MHz]。已知各段采样信号频率与首段采样信号频率成正比，其分频比例关系依次为[1 0.9 0.94 0.88]。已知采样信号的频率为8MHz，各段采样信号的长度依次为[50 51 49 50]。为精确测量首段采样信号频率，根据本专利方法，可采用如下处理过程：

(1)确定首段待比较频率序列BF0＝linspace(2.863636MHz，3.5MHz，207)。设定该序列长度为207以减小后续的Chirp-Z变换的运算量。

(2)确定各段采样信号的待比较频率序列BF(m，k)。

其中，BF(1，k)＝BF0(k)；BF(2，k)＝0.9·BF0(k)；

BF(3，k)＝0.94·BF0(k)；BF(4，k)＝0.88·BF0(k)。

(3)将各段同频信号进行BF(m，k)频率处的Chirp-Z变换，得到Z₀(m，k)。

(4)根据Z₀(m，k)利用求和取平均的方法得到首段采样信号频率的直接累加法估计值f₁₀。计算过程如下：

for k＝1:207

z1(k)＝abs(sum(Z₀(1:M，k)))；

end

[z2l，z22t]＝max(z1)；f₁₀＝BF0(z22t)；

(5)设 $\begin{array}{c}\mathrm{FS}(f_1\left(t\right),h,k)=\frac{\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\}}{2*\sin \{\mathrm{pi}\·[f_1-\mathrm{BF}0\left(k\right)]/f_s\}}\\ \·\frac{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\·\stackrel{\‾}{D}\·M\}}{\sin \{\mathrm{pi}\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]\·\stackrel{\‾}{D}/f_s\}}\end{array}$ ，其中pD=47，

f₁(t)∈[2.863636MHz，3.5MHz]，h，k∈[1，207]。调出预存的的极大值表S(f₁(t)，h)，该表为函数FS(f₁(t)，h，k)的最大值序列。取该函数表中f₁(t)∈(0.997*f₁₀，1.003*f₁₀)部分，得到S₂(f₁(f)，h)。

(6)计算 $S_3\left(h\right)=\underset{t}{\mathrm{\Σ}}\left|S_2\right[f_1(t+1),h]-S_2[f_1\left(t\right),h\left]\right|.$ 取S₃(h)中最大的10个值所对应的h

序列h₀(l)，取S₂(f₁(t)，h)中h∈h₀(l)部分，得到S₄(f₁(t)，h₀(l))。

(7)确定搜索频率序列SF0(l)＝BF0(h₀(l))。

(8)计算θ₂(m)＝angle[Z₀(m，z22t]。

(9)确定 $\left\{\begin{array}{c}\mathrm{\φ}(l,m,k)={\mathrm{\θ}}_2\left(m\right)-{\mathrm{\θ}}_2\left(1\right)-2\·\mathrm{pi}\·\mathrm{pB}\left(m\right)\·[\mathrm{SF}0\left(l\right)-\mathrm{BF}0\left(k\right)]/f_s\\ \mathrm{pB}\left(m\right)=\mathrm{sum}\left(\mathrm{pD}(1:m-1)\right)\end{array}\right.$

(10)计算如下：

$\left\{\begin{array}{c}{Z}_1(l,m,k)=\exp [-j\·\mathrm{\φ}(l,m,k)]\·Z_0(m,k)\\ Z_2(l,k)=\left|\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{Z}_1(l,m,k)\right|\end{array}\right.$

(11)取Z₂(l，k)中的谱峰频率位置Z₃(l)＝BF0(MAX[Z₂(l，k)])。

(12)根据Z₃和S₄计算得到频率f₃。 $\left\{\begin{array}{c}{Z}_4(f_1\left(t\right),l)=|Z_3\left(l\right)-S_4(f_1\left(t\right),h_0\left(l\right))|\\ Z_5\left(f_1\left(t\right)\right)=\underset{l}{\overset{{\mathrm{nf}}_1}{\mathrm{\Σ}}}{Z}_4(f_1\left(t\right),l)\end{array}\right.$ 取Z₅的最小值对应f₁(t₀₎，得到f₃＝f₁(t₀)；

(13)若f₃位于频率区间GF的左端处即f₃-gf₁≤0.01MHz，则左移频率区间GF和相应修改S₄，然后重复上述流程中的第(12)步。若f₃位于频率区间GF的右端处即gf₂-f₃≤0.01MHz，则右移频率区间GF和相应修改S₄，然后重复上述流程中的第(12)步。

(14)重复上述流程中的第(13)步，直至f₃不再位于频率区间GF的左右端。

接下来，结合该例，对本专利方法和目前适应性最好的直接累加法进行了全面的精度对比，给出了定量结果。仿真实验中使用的参数设置如下：噪声类型为高斯白噪声，信噪比SNR＝-5dB；各段采样信号的初相均满足幅度为2·pi的高斯分布。

(1)抗噪性仿真

该仿真过程包括500次实验，每次实验的参数除随机函数和信噪比外完全相同，每进行50次实验后信噪比改变一次，得到实验结果如表7所示。表7中的信噪比SNR单位为dB，标准差单位为KHz。可以看出，在多段分频信号、相同信噪比情况下，本专利方法大约能将精度提高到现有方法的6倍。

表7多段分频信号情况下两种方法的抗噪性能对比

(2)初始估计值f₂的相对误差对处理结果的影响仿真

该仿真过程包括500次实验，每次实验的参数除随机函数和初始估计值外完全相同，每进行50次实验后初始估计值改变一次，得到实验结果如表8所示。表8中的初始估计值f₂的相对误差C＝(f₂-f₁)/f₁，其单位为千分之一，标准差单位为KHz。可以看出，在多段分频信号、相同初始估计值情况下，本专利方法大约能将精度提高到现有方法的6倍。

表8多段分频信号情况下初始估计误差对两种方法的精度影响对比

(3)对不同首段采样信号频率f₁的适应性仿真

该仿真过程包括500次实验，每次实验的参数除随机函数和首段采样信号频率外完全相同，每进行50次实验后首段采样信号频率改变一次，得到实验结果如表9所示。表9中的首段采样信号频率f₁的单位为MHz，标准差单位为KHz。可以看出，在多段分频信号、相同首段采样信号频率情况下，本专利方法大约能将精度提高到现有方法的6倍。

表9多段分频信号情况下首段采样信号频率对两种方法的精度影响对比

(4)对不同采样频率f_s的适应性

该仿真过程包括500次实验，每次实验的参数除随机函数和采样频率外完全相同，每进行50次实验后采样频率改变一次，得到实验结果如表10所示。表10中的采样频率f_s的单位为MHz，标准差单位为KHz。可以看出，在多段分频信号、相同采样频率情况下，本专利方法大约能将精度提高到现有方法的6倍。

表10多段同频信号情况下采样频率对两种方法的精度影响对比

(5)对不同单段信号长度平均值D的适应性

该仿真过程包括500次实验，每次实验的参数除随机函数和单段信号长度平均值外完全相同，每进行50次实验后单段信号长度平均值改变一次，得到实验结果如表11所示。表11中的单段信号长度平均值D的单位为点，标准差单位为KHz。可以看出，在多段分频信号、相同单段信号长度平均值情况下，本专利方法大约能将精度提高到现有方法的6～7倍。

表11多段分频信号情况下单段信号长度平均值对两种方法的精度影响对比

(6)对不同信号总段数的适应性

该仿真过程包括500次实验，每次实验的参数除随机函数和信号总段数外完全相同，每进行50次实验后信号总段数改变一次，得到实验结果如表12所示。表12中的信号总段数M的单位为段，标准差单位为KHz。可以看出，在多段分频信号、相同信号总段数情况下，本专利方法大约能将精度提高到现有方法的6倍。

表12多段分频信号情况下信号总段数对两种方法的精度影响对比

(7)实时计算量对比

直接累加法实时计算量：4·M·(L·log₂(L)+L+1)，其中L≥nf₂+max(D(m))。

旋转相位积累法实时计算量：4·M·D·nf₂

本专利方法实时计算量：

4·M·(L·log₂(L)+L+1)+4·(M+nf₁·M·nf₂)+n·M·nf₁·0.02/0.0005，其中，nf₁一般取30即可；n表示重复上述流程中第(13)步骤的次数，一般取2即可；0.02表示上述流程中第(13)步骤中的频率搜索区域的长度；0.0005表示频率搜索中的频率步进幅度。

在本例多段分频信号的一般处理情况下，三种方法的实时计算量依次为36880、165600和184256，本专利方法的实时计算量增幅很小。

Claims (7)

1.一种信号频率估计方法，其特征在于，给出频域分析参数矩阵、相参补偿因子序列和搜索频率序列的计算公式，实现对多段同频信号、多段降频信号和多段分频信号的频率估计。

2.如权利要求1所述的多段同频信号，其特征在于，每个采样段内被测信号的频率不变，不同采样段的被测信号的频率相等。

3.如权利要求1所述的多段降频信号，其特征在于，每个采样段内被测信号的频率不变，不同采样段的被测信号的频率之间的差已知。

4.如权利要求1所述的多段分频信号，其特征在于，每个采样段内被测信号的频率不变，不同采样段的被测信号的频率之间的倍数已知，该频率之间的倍数可为任意小数。

5.如权利要求1所述的频域分析参数矩阵，其特征在于，为了处理多段不同频信号，根据首段被测信号的频率与其他段被测信号频率的对应关系，生成一个频域分析参数矩阵，根据该矩阵对其他段被测信号的频域分析过程所需的频域分析参数做相应调整，以实现同频化效果。

6.如权利要求1所述的相参补偿因子序列，其特征在于，以首段被测信号为基准生成一个相参补偿因子序列，再对其他段被测信号的频域分析结果乘以一个对应的相参补偿因子，能够补偿由于各段被测信号之间的频率不等、相位不连续对信号处理结果的影响。

7.如权利要求1所述的搜索频率序列，其特征在于，为获得首段被测信号频率计算值而生成一个搜索频率序列，由此可得到对应的首段被测信号频率计算值序列，根据该计算值序列可以求得首段采样信号频率以及其他段采样信号的频率。

Images