CN102175916B - 短样本密集频率信号的参数测量方法 - Google Patents

Abstract

一种短样本密集频率信号的参数测量方法：以采样频率f_s对输入的模拟信号x(t)进行等间隔采样获得2N-1个离散样本x(n)；对x(n)进行apFFT变换后得到复序列X_a(k)，取其模值得振幅谱后搜索出谱峰位置k＝k₀；为缩短频域搜索区间，取以k₀Δω为中心的频段ω∈[(k₀-1)Δω，(k₀+2)Δω]，Δω＝2π/N，将x(n)在该频段上做apDTFT变换，取此变换的相位谱；观察相位谱曲线特征，若相位谱曲线呈现平直分布特征，则为单频情况；若相位谱曲线呈现等幅振荡，则判断为双频分布情况，则进一步估测出两个密集成分的频率、幅值和相位值；除此之外则为成分数大于或等于3个成分的多频情况。本发明具有精度可控的特点。

Description

短样本密集频率信号的参数测量方法

技术领域

本发明涉及一种频率信号的参数测量。特别是涉及一种在采集到短样本时，对包含密集邻近成分的正弦波信号进行频率、初相位和幅值的数字测量的短样本密集频率信号的参数测量方法。

背景技术

多频正弦波的频率、幅值和初相的参数测量是工程界常常遇到的问题。目前，正弦信号的数字测量法以其精度高、方案灵活且无需引入复杂电路等优点成为广泛使用的测量手段。数字测量方法需要以一定采样速率对信号进行采样，再借助数字信号处理算法对采样数据进行分析计算，而得到频率、幅值和初相的估计值。因此，数字信号处理方法的优劣直接决定了最终参数估计精度和测量范围。

然而，在各领域的工程实践中，对于多频正弦波的参数测量，存在一个很棘手的难题。那就是信号的多种频率成分中，经常会遇到相邻的两个频率成分挨得非常近的情况(如在军事无线电环境中，接收载波经常会受到近频位置的敌方强干扰)，而且受采样时间所限，我们常常采集不到足够长的样本信号。如何分辨并测量短样本情况下的近频成分及其频率、幅值和初相参数一直是工程实践中的一个挑战性的难题，迫切需要解决。

近些年来，人们也一直致力于对信号参数测量方法的研究，研究成果层出不穷，比如：

(1)计数相加法(大石和明，石田秀树.频率测量电路[P].日本：008144012000.08.18.)，提出了频率测量单元在相互移位的周期内对参考时钟进行计数。设置加法器来把频率测量单元的计数相加。通过把在移位的周期内测得的计数值这样相加，得到频率测量结果。

(2)自适应频率测量器法(胡前南，曹全喜，龚洪金，刘亚威.自适应频率测量器[P]中国：2008202215832008.12.04.)。此装置包括采样转换电路，采样转换电路连接滤波电路，滤波电路连接数模转换电路，数模转换电路连接数字量处理单片机，数字量处理单片机连接人机接口。通过硬件电路完成比较宽范围的频率测量。

(3)采样统计法(严可国，金存山.一种频率测量方法和装置[P].中国：2008100571922008.01.30.)。此测量装置包括：接收被测信号，所述被测信号的边沿触发中断；当采样完成时，根据当前中断次数获取要求采样次数；根据定时中断次数和定时器时间，获取采样时间；根据所述要求采样次数和采样时间，获取所述被测信号的频率值。

以上种种方法均通过硬件电路实现，过分依赖专用的电路，因而会耗费大量硬件资源而且测量范围有限，精度难把握，并且当遇到两个频率成分挨得很近的情况上述方法显得无能为力。因此，近期人们把信号参数测量方法的由硬件逐渐转向软件，用数学方法通过软件编程实现参数的测量与分析。

DFT(Discrete Fourier Transform，离散傅里叶变换)和DTFT(Discrete Time FourierTransform，离散时间傅里叶变换)是最常用的正弦信号的参数估计法，DFT是DTFT在频域内的离散采样结果。两者都可以直接对采样序列进行变换而得到谱图，在谱图上直接搜索出峰值谱位置，再通过辨识谱峰位置而区分临近成分。

从谱分析的角度看来，一方面，样本长度直接决定了谱分析的阶数，短样本则意味着谱分析阶数N无法取得过高，因而很可能使得信号包含的两近邻频率成分的间距小于DFT的频率分辨率Δω＝2π/N，即属于短样本密集谱情况，显然由短样本而形成的密集谱会增大频率成分的识别难度；另一方面，DTFT与DFT的一个共同缺陷是，两种谱分析都存在很严重的谱泄漏效应。如Sanjit K.Mitra.DIGITAL SIGNAL PROCESSING：A Computer Based Approach，3rdEdition[M].NewYork：McGraw-Hill，2005、高西全，丁玉美.数字信号处理-第三版[M].西安：西安电子科技大学出版社，2008、以及D.Agrez.Improving phase estimation with leakageminimization[J]IEEE Trans.Instrum.Meas.，2005，54(4)：1347-1353所述。

而DFT相比于DTFT，虽然离散实现相对容易，但是存在栅栏效应而降低了谱分析性能。对于短样本密集谱情况，各频率成分的DFT谱和DTFT谱(包括振幅谱和相位谱)都挨得很近而引起非常大的谱间干扰(谢明，丁康.两个密集频率成分重叠频谱的校正方法[J].振动工程学报，1999，12(3)：109-114.方体莲，洪一.利用FFT校正两个密集信号的频率和相位[J].雷达科学与技术，2005，3(6)：378-382.)，这直接导致DFT谱和DTFT谱估计在识别密集频率成分时几乎失效。

对于密集谱参数识别问题，现代谱估计(Petre Stoica，Randolph Moses.Sepctral Analysis ofSignals[M].Upper Saddle River，Prentice Hall，2005.)方法(如AR模型方法(Kay S M.ModernSpcetral Estimation：Theroy and Application[M].Englewood Cliffs.NJ：Prentice-Hall，1988.)，最大熵谱估计方法(Burg J P.Maximum entropy spectral analysis[D].Ph.D dissertation，Dept.ofGeophysics，Stanford Univ.，CA，1975.))同样无法胜任。原因是现代谱估计是构筑在统计意义上的基于模型参数的谱分析方法，必须对大量样本进行统计分析(如计算自相关函数)才能得出精确的模型参数；另外，现代谱估计方法都是指功率谱估计(Petre Stoica，Randolph Moses.Sepctral Analysis of Signals[M].Upper Saddle River，Prentice Hall，2005.)，完全丢弃了信号的相位信息。

故以上各种方法过度强调了振幅谱(或功率谱)，单一地依靠搜索振幅谱峰的办法来获知频率位置。而对于密集成分样本不足的情况，由于各成分相互存在严重的谱间干扰，很可能会导致振幅谱出现单峰形状而无法辨别各子谱。因而密集谱估计一直是信号处理的难题。

为从谱图上来解释短样本和长样本下密集谱的区别，假定要估计x(t)＝cos(2πf₁t+θ₁)+2cos(2πf₂t+θ₂)的各频率成分，其中f₁＝3.3Hz，f₂＝3.5Hz，其x(t)在10s范围内的波形如图1(a)所示，以f_s＝32Hz的采样速率分别对x(t)采样，采集N＝32个样本得到的波形如图1(b)所示，采集N＝256个样本得到的波形如图1(c)所示。对这两段样本分别进行DFT和DTFT，得到的谱图分别如图1(d)和图1(e)所示。

观察图1(d)中余弦信号x(t)采样点数N＝32时的DFT和DTFT谱图：因为样本短，两频率分布在一个分辨率以内，频谱泄漏严重，只有一个谱峰，显然我们不能从这样的图谱中读出想要的频率成分，相位信息则被完全掩盖。图1(e)为余弦信号x(t)采样点数N＝256时的DFT和DTFT密集谱图，图中大致能区分出两密集成分所对应的两个谱峰，代价是耗费的采样点数是前者的8倍。

从图1可以看出，我们可以延长采样时间，或者搜集更多的样本来提高分辨率，从而减弱频率泄漏提高参数估计的准确率。但是，在工程实践中常常是得不到足够多的样本：比如，在军事上，时间就是生命，数据处理的过程中对速度和准确率又有极高的要求，很难保证充足的时间采集足够的样本做参数估计；又如在当今临床医学上，越来越要考虑到病人的忍受时限，激励信号通常不能持续过长，因而通常采集不到足够长的有用的响应数据，所以临床诊断分析也必须从少量仅有的样本入手。因此，对于短样本密集谱我们不得不选择其他途径，开辟新的观测角度来估测密集谱参数值。

事实上，信号的特征既体现在振幅谱中也包含在相位谱中。对于密集谱分布情况，在振幅谱分析失效的情况下，有必要考虑从相位谱中提取信息。全相位FFT(all-phase DFT，apFFT^{[12，13，14]})是近年来提出的一种新型的谱分析方法，文献[10]指出，全相位FFT具有比传统DFT更好的抑制谱泄漏特性，而且可直接提取出相位信息。

发明内容

本发明所要解决的技术问题是，提供一种计算量较小，估计效率高，资源耗费少的短样本密集频率信号的参数测量方法。

本发明所采用的技术方案是：一种短样本密集频率信号的参数测量方法，包括如下步骤：

1)首先以采样频率f_s对输入的模拟信号x(t)进行等间隔采样获得2N-1个离散样本x(n)；

1)对x(n)进行apFFT变换后得到复序列X_a(k)，取其模值得振幅谱后搜索出谱峰位置k＝k₀；

3)为缩短频域搜索区间，取以k₀Δω为中心的频段ω∈[(k₀-1)Δω，(k₀+2)Δω]，Δω＝2π/N，将x(n)在该频段上做apDTFT变换，取此变换的相位谱；

4)观察相位谱曲线特征，若相位谱曲线呈现平直分布特征，则判断为单频情况，这时不存在密集谱识别问题；若相位谱曲线呈现等幅振荡，且具有宽过渡带特征，则判断为双频分布情况；除此之外则为成分数大于或等于3个成分的多频情况；

5)若判断为双频情况，则进一步估测出两个密集成分的频率、幅值和相位值。

所述的apFFT变换，是用长为(2N-1)的卷积窗w_c对输入数据加权，然后将间隔为N的数据两两叠加，中间元素除外，而形成数据y(0)，y(1)，...，y(N-1)，再对这N个新数据进行FFT即得到apFFT结果X_a(0)，X_a(1)，...，X_a(N-1)，其中w_c由两个为N的矩形窗卷积而成。

所述的apDTFT变换，是把2N-1个输入数据分为N个子分段x₀、x₁、x₂、..、x_N-1。其中子分段x_m＝[x(-m)，x(-m+1)，...，x(-m+N-1)]^T，m＝0，1，...，N-1；进而分别求取各子分段的DTFT谱

m＝0，1，...，N-1；再对各子谱X_m(jω)进行累加求和，即得到apDTFT相位谱X_a(jω)。

所述的apFFT变换，是将输入信号x(n)经过三角窗加权后分别与e^-jωn，n∈[-N+1，N-1]相乘，再累积求和即得到等效的apDTFT相位谱X_a(jω)。

对双频情况，则进一步估测出两个密集成分的频率、幅值和相位值的过程如下：将频率观测区间进一步缩小为ω∈[(k₀+1)Δω，(k₀+2)Δω]，在该区间内继续做极值点搜索，搜到的第1个极值点ω₁则对应成分1的频率估计

该频率的极值相位则为成分2的相位估计

搜到的第2个极值点ω₂则对应成分2的频率估计

该频率的极值相位则为成分1的相位估计

然后将估计得到的两个频率成分的数字角频率值ω₁、ω₂及整数N值代入幅值估计公式

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\ω}}_2+\mathrm{\Δ\ω})\left]\right|}{\frac{{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\ω}}_1+\mathrm{\Δ\ω})\right]}{\frac{{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}}$

即分别得到这两个频率成分的幅值估计。

本发明的短样本密集频率信号的参数测量方法的有益效果如下：

1、填补了专门针对短样本密集谱情况，对信号初相位、频率等参数做有效精确估计方法的空白。

原有的数字信号处理的频谱估计大多是在有充足样本情况下，通过观测幅度谱或功率谱峰值的方式做频谱估计。它们没有把握相位谱中包含的信号参数信息，短样本密集谱情况下的参数估计问题，恰恰可以由相位谱解决，经查阅，现有文献没有涉及，本发明填补了这一空白。

2、本发明方法具有精度可控的特点。

实际应用中，可以根据实际工程中的精度需要，通过适当改变频率扫描步长调整分辨率来控制精度。

3、计算量较小，估计效率高，资源耗费少，大大节省了硬件成本。

在参数估计过程中做峰点搜索，将计算范围缩小到峰点k₀附近的一小段，对其余部分不做处理，减少了很多不必要的计算，这样大大加快了计算速度。

4、对于密集双频谱参数测量情况，有较高的微弱信号检测能力。

这是因为前面已经分析出，在获知密集谱属于双频谱的情况下，为精确得出两个密集成分的频率值，仅需搜索一个Δω区间内的apDTFT谱的相位谱的两个极值点位置即可实现，与两成分的幅值没有关系。这样在高信噪比情况下，对于一个强成分周围存在一个弱成分的情况，非常有利于对弱成分的检测。

5、有一定的抗噪声能力。

从表4可以表明，即使是对于图11所示的受到噪声严重干扰的信号，本发明方法仍可以完成精确的频率估计和较精确的相位估计。只要在全相位DTFT相位谱图中准确定位中心频率位置和极值点位置，即可快速准确的做出参数估计值。

附图说明

图1是某双频余弦信号不同样本长度的时域波形图和谱图；

图2是本发明的技术方案流程图；

图3是全相位FFT谱分析简化(N＝4)示意图；

图4是双频复指信号x(n)的apFFT幅度谱；

图5是从传统DTFT到apDTFT的衍生(N＝4)示意图；

图6是全相位DTFT谱分析简化(N＝4)示意图；

图7是直接DTFT谱和全相位DTFT谱(N＝8)波形图；

其中：

(a)由图3直接得到的DTFT振幅谱|Y(jω)|

(b)由图3直接得到的DTFT相位谱

(c)由图6简化的apDTFT流程得到的振幅谱|X_a(jω)|

(d)由图6简化的apDTFT流程得到的相位谱

图8是单频、双频、多频信号的全相位DTFT谱比较(N＝16)波形图；

其中：

(a)单频信号的apDTFT振幅谱    (b)单频信号的apDTFT相位谱

(c)双频信号的apDTFT振幅谱    (d)双频信号的apDTFT相位谱

(e)多频信号的apDTFT振幅谱    (f)多频信号的apDTFT相位谱

图9是双频信号的apDTFT谱，N阶DTFT谱和2N阶DTFT谱(N＝8)波形图；

其中：

(a)apDTFT振幅谱     (b)apDTFT相位谱

(c)N阶DTFT振幅谱    (d)N阶DTFT相位谱

(e)2N阶DTFT振幅谱   (f)2N阶DTFT振幅谱

图10是apDTFT相位谱的极值点和过渡带波形图；

其中：

(a)宽过渡带附近的apDTFT相位谱

(b)相位谱过渡带右侧第1极值点附近的局部放大图

(c)相位谱过渡带右侧第2极值点附近的局部放大图

图11是噪声干扰下相位谱图(信噪比SNR＝24dB)；

其中：

(a)宽过渡带附近的apDTFT相位谱

(b)宽过渡带附近的放大后的apDTFT相位谱

(c)宽过渡带附近的进一步放大后的apDTFT相位谱

图12是本发明的硬件实施例原理图；

图13是DSP内部程序流程图。

具体实施方式

下面结合实施例和附图对本发明的短样本密集频率信号的参数测量方法做出详细说明。

如图2所示，本发明的短样本密集频率信号的参数测量方法，首先以采样频率f_s对输入的模拟信号x(t)进行等间隔采样获得2N-1个离散样本x(n)；对x(n)进行apFFT变换后得到复序列X_a(k)，取其模值得振幅谱后搜索出谱峰位置k＝k₀；为缩短频域搜索区间，取以k₀Δω为中心的频段ω∈[(k₀-1)Δω，(k₀+2)Δω]，Δω2π/N，将x(n)在该频段上做apDTFT变换，取此变换的相位谱，观察其相位谱曲线特征：若相位谱曲线呈现平直分布特征，则判断为单频情况，这时不存在密集谱识别问题；若相位谱曲线呈现等幅振荡，且具有宽过渡带特征，则判断为双频分布情况；除此之外则为成分数大于或等于3个成分的多频情况。

若判断为双频情况，还可进一步估测出两个密集成分的频率、幅值和相位值。估测过程如下：将频率观测区间进一步缩小为ω∈[(k₀+1)Δω，(k₀+2)Δω]，在此在该区间内继续做极值点搜索，搜到的第1个极值点ω₁则对应成分1的频率估计

其极值相位则为成分2的相位估计

搜到的第2个极值点ω₂则对应成分2的频率估计

其极值相位则为成分1的相位估计

然后将估计得到的两个频率成分的数字角频率值ω₁、ω₂及整数N值代入幅值估计公式

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\ω}}_2+\mathrm{\Δ\ω})\left]\right|}{\frac{{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\ω}}_1+\mathrm{\Δ\ω})\right]}{\frac{{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}}$

即分别得到这两个频率成分的幅值估计。

下对图2所示的步骤的内部处理过程以及所涉及的理论原理逐一做说明，并对实验结果进行分析。

(1)全相位FFT

其过程如图3所示：只需用长为(2N-1)的卷积窗w_c对输入数据加权，然后将间隔为N的数据两两叠加(中间元素除外)而形成数据y(0)，y(1)，...，y(N-1)，再对这N个新数据进行FFT即得到apFFT结果X_a(0)，X_a(1)，...，C_a(N-1)。其中w_c由两个为N的矩形窗卷积而成。

本发明需借助全相位FFT(apFFT)以确定谱峰位置，进而缩小全相位DTFT的频率扫描范围，降低计算量。

严可国、金存山在“一种频率测量方法和装置”[P].中国：200810057192 2008.01.30中指出，全相位FFT具有很优良的抑制谱泄漏特性，故非常方便用于谱峰位置搜索，根据它的apFFT幅度谱来确定频率的大致范围是可信的。本发明的技术方案中就利用apFFT谱泄漏小的优良特性，对x(n)做变换后得其频谱序列，然后经过峰点搜索，确定整数频点位置k₀，则x(n)包含k₀Δω附近的密集频率，设置频段ω∈[(k₀-1)Δω，(k₀+2)Δω]为下面apDTFT变换的计算区间。我们用峰点搜索把频段锁定在一个很小的范围内，这样可以大大减少计算量。例如：以N＝8为例，对ω₁＝3.2Δω，ω₂3.5Δω，θ₁＝20°，θ₂＝70°作峰点搜索。由于|ω₁-ω₂|＜Δω，故属于密集谱分布情况。对x(n)做apFFT得到的归一化振幅谱|X_a(k)|如图4所示，对图4的谱序列做峰点搜索，可知谱峰位置k₀＝3，对应的频率为ω₀＝3Δω，因而x(n)的频率成分在ω₀＝3Δω附近。由此，设置频段区间ω∈[(k₀-1)Δω，(k₀+2)Δω]＝[2Δω，5Δω]，为下一步apDTFT相位谱分析的分析区间。

(2)全相位DTFT谱分析过程

如图5所示，全相位离散时间傅里叶变换(all-phase Discrete Time Fourier Transform，apDTFT)的过程是：把2N-1个输入数据分为N个子分段x₀、x₁、x₂、...、x_N-1。其中子分段x_m＝[x(-m)，x(-m+1)，...，x(-m+N-1)]^T，m＝0，1，...，N-1；进而分别求取各子分段的DTFT谱m＝0，1，...，N-1；再对各子谱X_m(jω)进行累加求和，即得到apDTFT谱X_a(jω)。

图5过程可简化为如图6流程所示：将输入信号x(n)经过三角窗加权后分别与e^-jωn，n∈[-N+1，N-1]相乘，再累积求和即得到等效的apDTFT谱X_a(jω)。

图5相比于图6其计算复杂度较高，而且存取这N个子DTFT谱需要较大的内存空间。而经过图6过程简化后，将N次DTFT运算降为1次，故仅需少量乘法和加法，从而大大节省了计算量和存储空间。需指出的是，由于生成的是连续谱，故图5和图6的中的频率变量ω都是需要以一定的频率步长进行扫描而得到的。

以上过程引入apDTFT的原因是：全相位FFT同传统FFT一样，是离散谱，存在栅栏效应，无法观测到单个频率分辨率内的谱分布情况(包括振幅谱和相位谱)。然而对于密集谱分布情况，关心的恰恰是一个频率分辨率内的谱分布信息，故本发明提出的apDTFT谱分析方法，这种新型谱分析方法可以在全频段上获得泄漏程度很小的振幅谱和直观的相位谱，非常易于提取密集谱信息。

注意apDTFT谱分析不能从图3直接引出，图3得到的输出X_a(k)是离散的，即只能在ω＝kΔω＝k2π/N(k＝0，1，...，N-1)上观测到谱值。为得到连续谱，一个容易直接想到的办法就是对图3中经全相位预处理后长度为N的数据y(0)～y(N-1)直接进行传统DTFT，即：

然而这样的直接做法是错误的。以下举例说明其错误性。

以信号

n∈[-N+1，N-1]，ω₀＝3.3Δω，Δω＝2π/N，θ₀＝20°，A＝1为例。由图3生成的N个数据y(0)～y(N-1)直接计算得到的Y(jω)的振幅谱|Y(jω)|如图7(a)所示，其相位谱

如图7(b)所示。

从图7(a)(b)可看出，虽然在ω＝kΔω上分别对|Y(jω)|和采样后即可得到apDFT振幅谱|X_a(k)|和相位谱

但是在ω≠kΔω的其他频点上，|Y(jω)|谱泄漏很大，

也不能直接反映真实初相值20°。因而图7(a)(b)的DTFT谱不是apDTFT谱。

经图6流程处理后，即可得到图7(c)所示的归一化后的apDTFT振幅谱|X_a(jω)|和图7(d)所示的相位谱

与图7(a)(b)一样，分别对|X_a(jω)|和相位谱

进行等间隔采样后，即得apDFT振幅谱|X_a(k)|和相位谱图7(c)与图7(a)不同的是，由于各个子谱在叠加过程中正负相互抵消，在ω≠kΔω的其他所有频点上，|X_a(jω)|的谱泄漏都很小，另外在全频段上，图7(d)的

都等于真实初相值20°。

(3)基于全相位DTFT相位谱的密集谱类型识别

该密集谱识别过程如下：观察由apDTFT得到的相位谱曲线特征：若相位谱曲线呈现平直分布特征，则判断为单频情况，这时不存在密集谱识别问题；若相位谱曲线呈现等幅振荡，且具有宽过渡带特征，则判断为双频分布情况；除此之外则为成分数大于或等于3个成分的多频情况。

举例来比较信号包含单频、双频、多频三种情况下的全相位DTFT相位谱的差别。假定单频、双频、多频三个信号分别是x₁(n)、x₂(n)、x₃(n)：

$x_1\left(n\right)=A_1{e}^{j({\mathrm{\ω}}_{1}n+{\mathrm{\θ}}_1)},$ ω₁＝(k₀+Δk₁)Δω；

$x_2\left(n\right)=A_1{e}^{j({\mathrm{\ω}}_{1}n+{\mathrm{\θ}}_1)}+A_2{e}^{j({\mathrm{\ω}}_{2}n+{\mathrm{\θ}}_2)},$ ω₁＝(k₀+Δk₁)Δω，ω₂＝(k₀+Δk₂)Δω；

$x_3\left(n\right)=A_1{e}^{j({\mathrm{\ω}}_{1}n+{\mathrm{\θ}}_1)}+A_2{e}^{j({\mathrm{\ω}}_{2}n+{\mathrm{\θ}}_2)}+,A_3{e}^{j({\mathrm{\ω}}_{3}n+{\mathrm{\θ}}_3),}$ ω₁＝(k₀+Δk₁)Δω，ω₂＝(k₀+Δk₂)Δω，

ω₃＝(k₀+Δk₃)Δω；

其中N＝16，Δω＝2π/16；A₁＝1，A₂＝2，A₃＝2；ω₁＝3.2Δω，ω₂＝3.5Δω，ω₃＝3.8Δω。

按照图6的简化分析流程分别对这三个信号作apDTFT和apFFT，得到的振幅谱和相位谱如图8所示。

这里信号分单频成分、双频成分和多频成分三种情况来讨论。

由于各密集成分频率间隙过小，全部挤在大小为Δω＝2π/N的小区间内，而且各成分的强度各有差异，故apDTFT谱跟传统DTFT谱一样，仅存在一个谱峰(如图8(a)(c)(d)的各振幅谱所示)。单纯依靠振幅谱去区分这三种情况则完全失效。因而有必要考虑相位谱，下面对三种频率分布的相位谱特征进行描述。

(i)单频情况

apDTFT的相位谱特征很明显，如图8(b)所示：相位谱

在全频率轴上呈现一根平直的线。在图8(b)上可直接读出单频信号的相位值θ₁＝20°，由单频信号振幅谱谱峰可读得频率成分ω₁＝3.2Δω。最后经过简单计算就得到信号的幅值参数A＝1。

(ii)双频密集谱分布情况

从图8(d)的相位谱图可看出，双频密集频率信号的相位谱有它自身的特点：宽过渡带、极值点周期分布，只有两个不同的极值(一个最大值，一个最小值)，如图中相位谱所示。我们也用这一特点判断一个信号是否双频。

由于相位谱图呈现周期分布和局部宽过渡带特征，因而在密集谱类型识别时，无需在整个频率区间内频率扫描，仅需在ω∈[(k₀-1)Δω，(k₀+2)Δω]内进行频率扫描即可。从图8(d)可看出，在该区间ω∈[(k₀-1)Δω，(k₀+2)Δω]＝[2Δω，5Δω]内，相位谱同时具备如下特征：(1)在ω＝k₀Δω＝3Δω附近呈现宽过渡带特性；(2)截断的周期性，即在ω∈[2Δω，5Δω]的宽度为3Δω内本应出现3个极值，由于中间宽过渡带的影响，减为2个极大值和2个极小值，两个极大或极小值的间隔为2Δω。若符合以上特征，即可断定为双频密集谱。

(iii)多频情况

若得到的apDTFT相位谱不符合以上描述的单频和双频成分的相位谱特征，如图8(f)所示，就可以断定信号包含3个或3个以上的密集频率成分。

(4)双频谱情况的两个成分频率、幅值和相位的参数估测

其过程如下：判断密集谱分布为双频谱情况后，将频率观测区间进一步缩小为ω∈[(k₀+1)Δω，(k₀+2)Δω]，在此在该区间内继续做极值点搜索，搜到的第1个极值点ω₁则对应成分1的频率估计其极值相位则为成分2的相位估计

搜到的第2个极值点ω₁则对应成分2的频率估计

其极值相位则为成分1的相位估计

然后将估计得到的两个频率成分的数字角频率值ω₁、ω₂及整数N值代入幅值估计公式

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\ω}}_2+\mathrm{\Δ\ω})\left]\right|}{\frac{{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\ω}}_1+\mathrm{\Δ\ω})\right]}{\frac{{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}}$

即分别得到这两个频率成分的幅值估计。

假定两个位于ω＝k₀Δω附近的密集频率分别是ω₁＝(k₀+Δk₁)Δω，ω₂＝(k₀+Δk₂)Δω，k∈z，k≠0不妨设ω₁＜ω₂，则有0＜Δk₁＜Δk₂＜1。可以证明，在ω＝(k+k₀+Δk₁)Δω，k≠0处，apDTFT的相位谱值即为第二个频率成分的初相值θ₂；同样地，在ω＝(k+k₀+Δk₂)Δω，k≠0处，apDTFT的相位谱值即为第一个频率成分的初相值θ₁；

为测量方便，在确知振幅峰值位置k₀后，令k＝1，则在ω＝(1+k₀+Δk₂)Δω处读出的相位谱值即可作为ω₁成分的初相值θ₁估计，在ω＝(1+k₀+Δk₁)Δω处读出的相位谱值即可作为ω₁成分的初相值θ₂估计。

因而在已经确知双频密集谱分布情况下，且已判断出峰值谱线位置为k＝k₀的情况下，为进一步准确地估计出双频谱成分的两个频率ω₁和ω₂值，可将频率扫描区间进一步缩小，仅只需在ω∈[(k₀+1)Δω，(k₀+2)Δω]扫描内即可。

在识别出谱分布为双频谱分布后，为估测密集谱参数，需做进一步的分析。为说明其原理，我们将全相位DTFT与传统DTFT做谱分析对比，对双频信号x₂(n)进行8阶apDTFT得到的振幅谱如图9(a)所示，相位谱如图9(b)所示；对双频信号x₂(n)进行8阶DTFT得到的振幅谱如图9(c)所示，相位谱如图9(d)所示；对双频信号x₂(n)进行8阶DTFT得到的振幅谱如图9(e)所示，相位谱如图9(f)所示；

从图9的振幅谱图可看出，由于两个成分的频率间隔只有0.3Δω，故图9(a)的apDTFT和图9(c)的8阶传统DTFT振幅谱都只有1个谱峰，即使把传统DTFT谱阶数增大1倍(相应耗费的样点数要多于apDTFT)，图9(e)的16阶DTFT振幅谱峰仍只有1个，并且谱泄漏没有得到明显改善。故依靠振幅谱分析结果，是无法区分出信号到底包含一个、两个还是多个频率成分。

从而不得不借助相位谱来识别密集谱成分。传统DTFT的相位谱图9(d)和图9(f)比较紊乱，没有明显特征。然而我们从图9(b)的apDTFT相位谱图中可看出，它具有周期性，在中心频率处有明显的过渡带，并且只有两个极值(一个最大值，一个最小值)。相位谱的这些特征表示它是双频谱，因此适合使用本发明中的参数估计方法进行估计。图9(b)的apDTFT相位谱图显示：在ω＝3Δω附近

存在一个较宽的过渡带，故可判断密集谱的整数频点位置为k₀＝3。这一点可以从图10更清楚地观察到。

图10(a)是双频信号apDTFT相位谱图的过渡带，对过渡带右侧的两个极值点局部放大，得图10(b)(c)。图10(b)观测到过渡带右侧的第一个极值点在ω＝4.2Δω处，而前面已分析该位置的理论值是ω＝(1+k₀+Δk₁)Δω，因而可确定Δk₁＝0.2，从图10(b)还可读出该频点

进一步可推知θ₂＝70°；此外图10(c)还可观测到过渡带右侧的第二个极值点位置在ω＝4.5Δω处，而前面已分析该位置的理论值是ω＝(1+k₀+Δk₂)Δω，因而可确定Δk₂＝0.5；从图10(c)还可读出该频点故同样可推知θ₁＝20°。最后，确定出θ₁＝20°，ω₁＝(k₀+Δk₁)Δω＝3.2Δω，θ₂＝70°，ω₂＝(k₀+Δk₂)Δω＝3.5Δω。从以上分析过程中可发现，在确定各成分的频率值的过程中，只需寻找极值点位置即可实现，也就是说两个密集成分的频率值的确定完全不受两成分的幅值影响。

根据以上过程测出两个频率值ω₁和ω₂后，将阶数N和Δω＝2π/16的值代入下式便可测算出两个频率成分的幅值，即

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\ω}}_2+\mathrm{\Δ\ω})\left]\right|}{\frac{{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_2+\mathrm{\Δ\ω}-{\mathrm{\ω}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\ω}}_1+\mathrm{\Δ\ω})\right]}{\frac{{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}{N^2{\sin }^2[({\mathrm{\ω}}_1+\mathrm{\Δ\ω}-{\mathrm{\ω}}_2)N/2]}}$

从图9(a)振幅谱图中可读出|X_a[j(ω₂+Δω)]|＝0.0262，|X_a[j(ω₁+Δω)]|＝0.1131，从而根据上式即可计算出A₁＝1，A₂＝2。

至此，完成对密集双频信号x₂(n)的频率、初相位和振幅等参数的估计。从整个过程可以看出本发明所述参数估计方法对短样本密集谱有很强的针对性，能快速准确地对单频双频信号做出精确的参数估计，并判断出其频率成分是否多于2个。同时，在不同场合，根据对精确度的不同要求可以通过综合调节频段ω长度和扫描步长控制测量精度。

因此本发明是个节省计算量、精度可控的优良参数估计方法，具有很强的工程应用性。

(5)实验结果

工程应用中，信号总是受到噪声干扰的，也可能是两相邻密集成分的能量差别很大。另外，实际工程中遇到的信号通常不会是理想的复指数信号，而是实余弦信号(包含有两个共轭的复指数信号成分)。因而有必要通过实验研究来分别比较无噪声和有噪声条件下单频复指数信号，双频密集复指数信号和双频密集余弦信号的频率和初相估计精度的差异。

实验1：

对信号

$x\left(n\right)=e^{j({\mathrm{\ω}}_{0}n+{\mathrm{\θ}}_0)},$ ω₀＝3.3Δω，Δω＝2π/N，θ₀＝20°； $y\left(n\right)=e^{j({\mathrm{\ω}}_{1}n+{\mathrm{\θ}}_1)}+2e^{j({\mathrm{\ω}}_{2}n+{\mathrm{\θ}}_2)};$

z(n)＝cos(ω₁n+θ₁)+2cos(ω₂n+θ₂)，ω₁＝3.2Δω，ω₂3.5Δω，θ₁＝20°，θ₂＝70°。

在无噪情况下按本发明方法进行频率和相位的估测，表1给出了不同谱分解阶数N情况对应的参数估计结果。

表1无噪情况下本方法的参数测量结果(扫描步长：Δω/10³，Δω＝2π/N)

从表1中可看出，用本发明的方法进行的参数测量，对于单频复指数信号x(n)，在谱阶数N＝8，即仅仅用了2N-1＝15个样本时，就能得到精确的频率估计结果3.3Δω和相位估计结果20°；对于双频复指数信号y(n)，在谱阶数N＝8，即仅仅用了2N-1＝15个样本时，也能得到精确的频率估计结果3.2Δω、3.5Δω和初相位估计结果20°、70°；对于双频余弦信号z(n)，由于余弦信号存在两个边带的相互干扰，因而其精度相比于复指数情况稍降低些(频率估测误差仅为0.03Δω，相位估测误差仅为2°左右)，但其频率成分的测量误差处于10^-2倍分辨率Δω，在工程实践中也常常是符合要求的。

实验2：

将扫描步长改为：Δω/1007，Δω＝2π/N，重复实验1，其得到的实验数据如表2所示。

表2无噪情况下本方法的参数测量结果(扫描步长：Δω/(1007)，Δω＝2π/N)

表2得到的测量结果相对于表1有一定的偏差，这是由于扫描步长限制，扫描频率ω没有遍历到理想的频率位置的缘故。比较表1与表2中的实验结果可以得出结论，在确定Δω后，频率扫描步长最好选择为Δω/10ⁿ，n∈z，以利于精度分析。这里注意Δω与模拟频率分辨率的对应关系，即Δω对应的模拟频率分辨率为Δf＝f_s/N，即由采样频率与谱分析阶数共同决定。

实验3

上述实验中的双频信号参数测量均是在两个种频率成分的幅值差别不大的情况下进行的，当两个频率成分的幅值相差很大时，传统振幅谱会出现强成分的泄漏掩盖弱信号成分现象，导致幅值大的强信号很容易掩盖幅值小的弱信号的参数信息。

为此，估测如下信号的频率和幅值参数值

$y\left(n\right)=A_1{e}^{j({\mathrm{\ω}}_{1}n+{\mathrm{\θ}}_1)}+A_2{e}^{j({\mathrm{\ω}}_{2}n+{\mathrm{\θ}}_2)}{e}^{,}$ z(n)＝A₁cos(ω₁n+θ₁)+A₂cos(ω₂n+θ₂)

表3记录了用本发明的方法分析无噪情况下，A₁、A₂具有不同幅值差异的双频谱信号的频率和幅值检测结果。

表3无噪情况下本方法的参数测量结果(扫描步长：Δω/10³，Δω＝2π/N，N＝16)

从表3可看出，对于双频复指数信号，无论A₂为A₁的2倍、100倍甚至是1000倍，本方法都能很精确地估计出频率和幅值参数值；对于双频余弦信号，由于余弦信号的两个边带间的相互影响，相比于复指数情况，其估计精度有所降低，其中弱信号成分参数相比于强信号成分的参数估计精度要更低些：在A₂为A₁的1～15倍的范围内，其弱成分频率估测精度可控制在0.1Δω以内，幅值估计误差能控制在10％左右，基本能符合工程实践要求，而其强信号成分的频率估计基本无偏差，幅值估计误差也很小。

本发明之所以能在两个密集成分能量差异很大的情况下仍能完成参数估测，是由于本发明依据相位谱进行估测的缘故，避开了传统方法仅依据振幅谱在弱信号检测时的缺陷。

实验4

以上讨论仅限于无噪情况。对于存在噪声干扰情况，以双频复指数信号和双频余弦信号为例，仍用本发明方法对实验1的信号进行参数检测与误差分析。其中双频密集谱信号y(n)的apDTFT相位谱图如图11所示。

实验所加噪声为高斯白噪声，信噪比24dB。从图11中可以看出，该噪声对相位谱图产生了影响，但是过渡带依然很明显，从图11可以读出相关频偏值和初相位。实验还进行了有噪条件下的单频复指数信号、双频余弦信号的参数估计，对三种不同类型信号的估计结果统计后填入表4。

表4噪声情况下本方法的参数估计结果(扫描步长为π/2000，信噪比SNR＝24dB)

从表4可看出，本测量方法具有一定的抗噪能力，表现在对频率分量的估测上。对于双频复指数信号，当谱阶数N＝8时(样本数目为15)，其弱成分频率估计误差为0.07Δω左右，强成分的频率估计误差由于自身能量强，抵御噪声干扰能力大，其误差仅为0.003Δω。当谱阶数N＝16时(样本数目为31)，强弱信号的频率估计精度进一步提高。对于双频余弦信号，由于存在两个边带的相互干扰，因而在同样的阶数下，频率估计精度比双频复指数信号更低些。

另外，从表4还可看出，所有情况的幅值估计精度比频率估计精度都要低，这很容易从式(4)的可以看出，幅值是由频率估计算出来的，因而频率估计的误差会代入到幅值估计中去。

下面对实施本发明的硬件予以简单说明。如图12所示，待测信号经过A/D(模数转化器)采样得到样本序列x(n)，以并行数字输入的形式进入DSP器件，经过DSP器件的内部算法处理，得到信号的参数估计；最后借助输出驱动显示及其显示模块显示出调频率和中心频率的估计值。

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

(1)调用核心算法，完成接收信号的参数估计处理；

(2)根据实际需要调整采样率f_s(改变图12的CP₂时钟频率)，使之在该采样率条件下的频率分辨率Δω，尽量地满足实际需要；

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

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

需指出的是，本发明所涉及的是短样本参数测量，因而并不需要很大的存储空间存储采集到的样本数据(仅需DSP的数据内部RAM的空间)。图12之所以接外部RAM是因为在样本量很少的情况下，要想得到足够精确的谱图，必须将频率扫描点取得很密集，这样就需消耗一定的存储空间用于存放apDTFT处理后的数据。

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

本发明将所提出的“新型短样本密集谱参数估计方法”这一核心估计算法植入DSP器件内，基于此完成高精度、低复杂度、高效的频率、幅值和相位参数估计。

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

(1)首先需根据具体应用要求(如医学和军事等的具体测量要求)，粗略估计信号频率分辨率Δω并根据具体需要设定相位精度要求。该步骤是从工程方面提出具体需求，以使得后续流程有针对性地进行处理。

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

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

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

(5)判断本发明方法是否满足工程需求，若不满足，程序返回，重新根据要求设定频率分辨率。

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

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

Claims (4)

1.一种短样本密集频率信号的参数测量方法，其特征在于，包括如下步骤： 

1）首先以采样频率f_s对输入的模拟信号x(t)进行等间隔采样获得2N-1个离散样本x(n)； 

2）对x(n)进行apFFT变换后得到复序列X_a(k)，取其模值得振幅谱后搜索出谱峰位置k=k₀； 

3）为缩短频域搜索区间，取以k₀Δω为中心的频段ω∈[(k₀-1)Δω,(k₀+2)Δω]，Δω=2π/N，将x(n)在该频段上做apDTFT变换，取此变换的相位谱； 

4）观察相位谱曲线特征，若相位谱曲线呈现平直分布特征，则判断为单频情况，这时不存在密集谱识别问题；若相位谱曲线呈现等幅振荡，且具有宽过渡带特征，则判断为双频分布情况；除此之外则为成分数大于或等于3个的多频情况； 

5）若判断为双频情况，则进一步估测出两个密集成分的频率、幅值和相位值。 

2.根据权利要求1所述的短样本密集频率信号的参数测量方法，其特征在于，所述的apFFT变换，是用长为(2N-1)的卷积窗w_c对输入数据加权，然后将间隔为N的数据两两叠加，中间元素除外，而形成数据y(0),y(1),…,y(N-1)，再对这N个新数据进行FFT即得到apFFT结果X_a(0),X_a(1),…,X_a(N-1)，其中w_c由两个为N的矩形窗卷积而成。 

3.根据权利要求1所述的短样本密集频率信号的参数测量方法，其特征在于，所述的apDTFT变换，是把2N-1个输入数据分为N个子分段x₀、x₁、x₂、…、x_N-1，其中子分段x_m=[x(-m),x(-m+1),…,x(-m+N-1)]^T，m=0,1,…,N-1；进而分别求取各子分段的DTFT谱 

，m=0,1,…,N-1；再对各子谱X_m(jω)进行累加求和，即得到apDTFT相位谱X_a(jω)。 

4.根据权利要求1所述的短样本密集频率信号的参数测量方法，其特征在于，所述的apFFT变换，是将输入信号x(n)经过三角窗加权后分别与e^-jωn,n∈[-N+1,N-1]相乘，再累积求和即得到等效的apDTFT相位谱X_a(jω)。 

Images