CN101833035B - 线性调频信号参数估计方法及其实施装置 - Google Patents

Abstract

本发明属于数字信号处理技术领域。为提供一种高精度、低复杂度、快速、高效的线性调频信号参数估计方法及其实施装置，本发明采用的技术方案是，借助基于能量重心插值估计和Radon模糊变换RAT的调频率估计方法，具体为：对于给定输入信号s(n)，首先求取输入信号s(n)的模糊函数AF_s(τ，ξ)，然后不断更换扫描直线ξ＝kτ，借助Radon变换计算得出AF_s(τ，ξ)在扫描直线ξ＝kτ上的投影，并记录所有投影值，构造全景离散谱，接着搜索出该离散谱的谱峰位置，估计实际峰值投影值，最后得出信号调频率的估计值

然后，搜索出分数阶傅立叶谱的峰值谱，计算出谱能量重心位置

再借助

即可得到信号重心频率的估计。本发明主要应用于数字信号处理中的线性调频信号参数估计。

Description

线性调频信号参数估计方法及其实施装置

技术领域

本发明属于数字信号处理技术领域。具体涉及一种高精度地估计线性调频信号两个最基本参数——调频率和中心频率的方法，即线性调频信号参数估计方法及其实施装置。

背景技术

线性调频(linear frequency modulation，LFM)信号是一种具有特殊意义的非平稳信号。作为一种大时间-带宽积的扩频信号，LFM信号广泛地应用于各种信息系统中。在雷达、声纳等探测系统中，目标的多普勒频率与目标的速度近似成正比，当目标作等加速运动时，回波即为线性调频信号。机动目标的回波在一段短的时间内也可用LFM信号作为其一阶近似。在通信系统中，线性调频可作为实现扩频通信的一种手段。此外，在地质勘探、生物医学和机械设计与加工等领域中，许多信号本身就具有线性调频性质。因此，对LFM信号的分析与处理方法的研究具有重要的意义。

在LFM信号的处理中，信号的检测与参数估计是一个重要的研究课题，这一问题在本质上可归结为非平稳信号的瞬时频率估计问题。因此，有必要首先对瞬时频率的定义及其估计方法的发展进行简要的回顾。

解析形式的非平稳信号可以表示为[1]

s(t)＝a(t)e^jφ(t)         (1)

其中a(t)为实的包络函数，φ(t)为瞬时相位函数。信号的瞬时频率定义为瞬时相位函数对时间的导数，即

$f\left(t\right)=\frac{1}{2\mathrm{\π}}\·\frac{\mathrm{d\φ}\left(t\right)}{\mathrm{dt}}---\left(2\right)$

式(1)、式(2)的离散形式为

s(n)＝a(n)e^jφ(n)           (3)

$f\left(n\right)=\frac{1}{2\mathrm{\π}}\·\mathrm{\φ}\left(n\right)*d\left(n\right)---\left(4\right)$

式(4)中的d(n)是一个FIR差分滤波器的冲击响应，*表示卷积运算。

实际系统中，人们得到的观测信号不可避免地受到各种噪声的污染，而所谓瞬时频率的估计即为由有限长的观测信号对瞬时频率的变化规律进行估计。长期以来，这个问题一直是信号处理中备受关注的研究课题之一，有关的研究成果也极为丰富[2]。归纳起来，瞬时频率估计的各种方法可分为两类，即非参数化估计和参数化估计。

非参数化估计直接利用式(1)和式(2)给出的定义来估计信号瞬时频率，常用的非参数化估计方法有相位差分估计、过零检测估计和自适应估计等。由相位差分估计可以得到瞬时频率的无偏估计，其运算量也较小，但差分运算导致了高频噪声的增大，从而影响了估计的精度[3-4]。过零检测估计的计算也非常简单，对于平稳信号具有较好的估计性能，但对于非平稳信号，其估计精度不高[5]。自适应估计包括锁相环法[2]、扩展Kalman滤波法[6]、LMS[7]和RLS[8]自适应滤波法等，其主要特点是对于频率变化较慢的准平稳信号具有较高的估计精度，并且算法简单，便于实现；但当瞬时频率的变化速率较快时，其性能下降。此外，非参数化估计的一个致命的缺陷是不能适用于多分量信号。近年来，随着时频分析方法研究的不断深入，基于时频分析工具的瞬时频率估计方法也不断出现。文献[9-10]提出，利用某些时频分布(如维格纳分布WVD和短时傅立叶变换STFT)的一阶矩可有效地估计非平稳信号的瞬时频率，但由于要计算二维的时频分布及其一阶矩，这种方法的计算量非常大。

参数化瞬时频率估计依赖于被估计信号的某些先验知识，通常是用某种确定的模型来描述信号，通过对模型参数的估计来实现对瞬时频率变化规律的估计。最常用的模型形式是多项式相位(Polynomial Phase)模型，在这种模型中，用时间变量的有限次多项式函数描述信号的瞬时相位，即

$\mathrm{\φ}\left(t\right)=a_0+a_{1}t+a_2{t}^2+\·\·\·+a_P{t}^P=\underset{k=0}{\overset{P}{\mathrm{\Σ}}}{a}_k{t}^k---\left(5\right)$

或

$\mathrm{\φ}\left(t\right)=a_0+a_{1}n+a_2{n}^2+\·\·\·+a_P{n}^P=\underset{k=0}{\overset{P}{\mathrm{\Σ}}}{a}_k{n}^k---\left(6\right)$

通过对参数a_k(k＝0，1，...，p)的估计即可得到φ(t)或φ(n)的估计。本发明中仅讨论具有恒定幅度的LFM信号，即限定P＝2，a(t)或a(n)为常数。

对于LFM信号的参数估计，最大似然估计是解决这一问题的主要途径，其中基于时频平面峰值检测的算法具有较好的性能。常用的时频分析工具包括线性时频表示，如STFT(短时傅立叶变换，Short Time Fourier Transform)和小波变换(WT，wavelet Transform)[11]，以及二次型的非线性时频分布，如WVD(Wigner-Ville Distribution，维格纳分布)[12-13]、广义WVD[14]、广义模糊函数[15]、多项式相位变换(Polynomial-Phase Transform)[16]等。但STFT窄的观察窗和WT宽度变化的时间窗影响了时频域的分辨率。而在非线性的时频分布中，WVD对LFM信号具有最好的时频聚集性，非常适合于对LFM信号的处理。但由于其变换过程的非线性，在利用这类方法处理多分量信号时，必然会受到交叉项的困扰，虽然可通过选择合适的核函数来抑制交叉项，但同时也降低了信号的时频聚集性[17]。文献[18-19]提出了基于WVD-HT的多分量LFM信号的检测和参数估计方法，借助于WVD-HT(WVD-Hilbert变换，维格纳分布-希尔伯特变换)的线积分过程，该方法可有效地抑制交叉项的干扰。这些方法在本质上大都可归结为一个多变量的最优化问题，但常因计算量太大，使得算法的工程实现较为困难，文献[20-21]提出的基于相位解缠绕技术的快速算法可显著地降低运算量，但不适用于多分量信号。

近年来，一种新的时频分析工具——分数阶Fourier变换(fractional Fourier transform，以下简称FRFT)引起了越来越多的关注[22]，由于FRFT是一维线性变换，与二次型时频分布相比，它不受交叉项困扰，而且可以理解为chirp基分解[23]，具有计算量与FFT相当的快速算法[24]，所以非常适合于LFM信号的参数估计，文献[25]也证明了基于FRFT的LFM信号参数估计方法的有效性。该算法利用FRFT对LFM信号良好的能量聚焦性，将信号的参数估计问题转化为分数阶Fourier域上的二维谱峰搜索过程，具有鲁棒性强、物理意义明确等优点。但在实际应用中，只能得到有限长样本的离散FRFT谱，这就不可避免存在由于时域截断产生的能量泄漏，使谱峰值变小、分辨率降低。另外，离散谱分析必然会受到栅栏效应的影响，导致峰值谱线位置偏离真实位置，若直接以离散谱峰扫描的方法估计谱峰位置，最大估计误差可能达到谱线间隔的50％。在对LFM信号参数精度要求比较高的应用场合，只能细化扫描步长、增大信号的采样率，但这样做会使计算量呈级数增加，不利于工程实现。

发明内容

本发明的目的在于克服现有LFM信号参数估计技术的各种缺陷，提供一种高精度、低复杂度、快速、高效的LFM信号参数估计方法及其实施装置。

为达到上述目的，本发明采用的技术方案如下：

线性调频信号参数估计方法，包括下列步骤：

借助基于能量重心插值估计和Radon模糊变换RAT的调频率估计方法，具体为：对于给定输入信号s(n)，首先求取输入信号s(n)的模糊函数AF_s(τ，ξ)，然后从μ＝0开始以Δμ为步长等间隔地在模糊域(τ，ξ)平面内构造过原点的扫描直线(τ，ξ分别为模糊参数)，假定第k次斜率为kΔμ的扫描直线为ξ＝(kΔμ)τ，然后将模糊函数AF_s(τ，ξ)在该直线上进行投影得到投影值η(k)，并记录所有投影值η(k)以构造全景离散谱，接着搜索该离散谱峰位置再借助能量重心插值公式

μ_k＝k*+kΔμ(其中n为以峰值谱为中心的左右投影值数目)即得出信号调频率的估计值

然后，用得到的调频率估计值

作为参数对信号做分数阶傅立叶变换，而得到分数阶傅立叶变换的离散谱S(u_i)，其中u＝iΔu，Δμ为离散谱间隔，进而找到其峰值谱位置u＝mΔu，借助

u_i＝iΔu计算出谱能量重心位置，再借助即可得到信号重心频率的估计。

高效的LFM信号参数估计实施装置，包括：

信号调理电路，用于对输入信号进行模拟预处理，以对信号幅度范围进行必要调整，并去除外干扰噪声；

A/D模数转化器，用于采样得到样本序列s(n)，以并行数字输入的形式送入DSP器件；经过DSP器件的内部的算法处理，而得到信号参数的估计，

数字信号处理器，用于：

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

(2)根据实际需要调整采样率f_s，使得在该采样率下，尽量高精度地估计出信号参数；

(3)内部RAM存储数据不足时，将处理数据与外部RAM进行数据交换，以配合核心算法处理；

(4)将相位估计结果实时输出到驱动和显示模块；

显示模块及其输出驱动，用于显示出调频率和中心频率的估计值。

数字信号处理器进一步包括：①数据预处理模块；②RAT变换处理模块；③RAT谱峰的搜索和校正模块；④分数阶Fourier变换FRFT变换处理模块；⑤分数阶Fourier变换FRFT谱峰的搜索和校正模块。

本发明提出的基于RAT和FRFT谱校正的LFM信号参数估计方法，应用在各工程领域，可以产生如下的有益效果。

第一，具有较高的估计精度，适合于对LFM信号参数估计精度要求比较高的应用场合：

本发明将能量重心插值估计措施引入到RAT域和FRFT域上的两次谱峰搜索过程中，克服了离散谱分析所固有的栅栏效应、截断效应对谱峰位置估计精度的影响，实现了谱峰位置的超分辨率估计，直接提高了LFM信号的调频率和中心频率的估计精度。

第二，估计效率高、速度快，适合于各种实时快速估计的应用场合：

通过前面的理论分析可以看到，本发明提出的估计方法将LFM信号的调频率和中心频率的估计问题转化为信号的RAT域和FRFT域上的两次一维谱峰搜索过程，并将能量重心谱插值估计引入两次搜索过程。算法简单、参数估计的速度较快，尤其可以通过对谱峰位置的实时变化情况连续估计接收信号的参数变化情况，具有比较好的实时估计性能。

第三，适合于多分量LFM信号的参数估计应用场合：

实际应用中，接收信号通常包含多个LFM分量，本发明提出的估计方法利用FRFT对LFM信号的良好能量聚集特性，实现了各个LFM信号分量之间的分离，然后逐一估计出各信号分量的参数。

第四，具有较强的抗噪声能力：

由FRFT的定义可知，一个LFM信号在适当的分数阶域中是一个冲击函数，即FRFT在某个分数阶Fourier域中对给定的LFM信号具有非常强的能量聚集特性，而白噪声的能量则均匀地分布在整个时频平面上。这一特性将使得该方法具有非常优良的抗噪声性能。

第五，可大大节省硬件成本：

通过上述的理论分析可知，本发明提出的估计方法的运算量主要集中在对接收信号进行的RAT以及FRFT两次变换过程中。实际工程中，通常利用模糊函数和分数阶相关的关系来简化RAT的运算过程[26]，从而将RAT表示为由FRFT和Fourier变换组成的等价形式。另外，由于FRFT具有计算量与FFT相当的快速算法[24]，且最后的谱线校正措施几乎不需要消耗计算量，所以本发明提出的估计方法的整体所需计算量很小，其核心的运算模块只需一块中、低端的DSP芯片即可实现，这大大节省了硬件成本。

附图说明

图1离散谱校正的能量重心法示意图。

图2基于能量重心插值与RAT变换的调频率估计流程。

图3基于能量重心插值与分数阶傅立叶变换的中心频率估计流程。

图4本发明硬件框图。

图5DSP器件的内部程序流程。

图6单分量LFM信号的模糊函数。

图7单分量LFM信号的RAT谱。

图8信号能量最聚集域的FRFT谱。

具体实施方式

1、基于Radon变换的LFM信号参数估计原理：

首先给出理想的无限长LFM信号的数学模型

s(t)＝a₀exp[jπ(2f₀t+μ₀t²)]，-∞≤t≤∞           (7)

其中a₀、f₀和μ₀为未知参数，分别代表信号的幅度、中心频率和调频率，其中中心频率和调频率为工程实践中比较关心的参数。根据定义，信号s(t)的模糊函数(ambiguity function，AF)表达式为

${\mathrm{AF}}_s(\mathrm{\τ},\mathrm{\ξ})={\∫}_{-\∞}^{+\∞}s(t+\mathrm{\τ}/2)s*(t-\mathrm{\τ}/2)e^{-j2\mathrm{\π\ξt}}\mathrm{dt}---\left(8\right)$

相应的(τ，ξ)平面通常称为模糊域。将式(7)代入上述定义式并取模

|AF_s(τ，ξ)|＝|a₀ ²exp{j2πf₀τ}δ(ξ-μ₀τ)|＝a₀ ²δ(ξ-μ₀τ)     (9)

由上式可知，理想LFM信号的AF的模在模糊域上是一条过原点的直线型冲击函数，且该直线的斜率就是LFM信号的调频率。对式(9)做过原点斜率为k的直线ξ＝kτ上的Radon变换(Radon Transform)，不断变化直线的斜率k，形成如下的检测统计量

$\mathrm{\η}\left(k\right)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}\left|{\mathrm{AF}}_s(\mathrm{\τ},\mathrm{\ξ})\right|\mathrm{\δ}(\mathrm{\ξ}-\mathrm{k\τ})\mathrm{d\ξd\τ}={\∫}_{-\∞}^{\∞}\left|{\mathrm{AF}}_s(\mathrm{\τ},\mathrm{k\τ})\right|\mathrm{d\τ}---\left(10\right)$

显然当直线斜率为μ₀时，η(k)的值最大。利用这一特性，可以通过对η(k)的谱峰位置搜索实现LFM信号的调频率估计。用公式表示为

${\widehat{\mathrm{\μ}}}_0=\arg \underset{k}{\max }\mathrm{\η}\left(k\right)---\left(11\right)$

以上就是基于Radon-Ambiguity Transform(以下简称RAT)的LFM信号调频率估计原理。实际工程中，通常利用模糊函数和分数阶相关的关系来简化RAT的运算过程[26]。

以上推导是理想情况下的检测结果。实际工程实现时，扫描曲线ξ＝kτ只能是以一定的扫描间隔Δμ通过多次扫描遍历而寻最优，则第k次扫描的曲线为ξ＝kΔμτ，故这很有可能所有整数倍扫描斜率kΔμ都没法落在理想斜率μ₀上，故本专利采用能量中心插值方法来估计μ₀，即找出峰值投影值周围几根投影值的能量中心位置作为最终的调频率估计

由于我们已经得到了观测信号调频率的估计值，根据文献[27]的结论，可以计算出该LFM信号能量最聚集域的FRFT角度为：

$\mathrm{\α}={\mathrm{\α}}_0=-\mathrm{arccot}\left({\widehat{\mathrm{\μ}}}_0\right)---\left(12\right)$

对信号做角度为α₀的FRFT，得到能量最聚集域上的FRFT谱

。通过对

的谱峰位置搜索即可实现LFM信号的中心频率估计，用公式表示为

${\hat{f}}_0={\hat{u}}_0\csc {\mathrm{\α}}_0---\left(13\right)$

${\hat{u}}_0=\arg \underset{u}{\max }\left|S_{{\mathrm{\α}}_0}\left(u\right)\right|---\left(14\right)$

以上就是基于FRFT的LFM信号中心频率估计原理。

综上所述，LFM信号的调频率和中心频率估计问题可以转化为式(11)和式(14)所示的两个一维谱峰搜索问题。实际应用中，我们只能得到信号RAT域和FRFT域上的两个离散谱。由于受到栅栏效应的影响，在真实的谱峰位置未对准谱线间隔的整数倍位置时，若采用传统的扫描方法，以峰值谱线的位置作为真实谱峰位置的估计，则最大可能造成50％谱线间隔的误差。在对LFM信号参数精度要求比较高的应用场合，只能细化扫描步长、增大信号的采样率，但这样做会使计算量呈级数增加，不利于工程实现。所以需要寻找一种更加简单快捷的超分辨率谱峰搜索方法。

2、离散谱校正的能量重心法

离散谱校正理论是为了解决正弦信号的频率和相位估计问题而提出的。由于计算机只能对正弦信号的有限个样本进行运算，FFT和谱分析也只能在有限区间内进行，这就不可避免地存在由于时域截断而产生的能量泄漏，使谱峰值变小、精度降低，即经FFT得到的离散频谱的幅值、相位和频率都可能产生较大的误差。另外，由点数为N的FFT处理得到的离散谱受频率分辨率Δω＝2π/N的限制，很难保证信号的真实频率处在Δω的整数倍位置，所以无法通过峰值搜索获得准确的频率和相位值。

文献[28]通过对常用窗函数的能量谱特性的分析，推导出各种对称窗函数谱的能量重心都位于它们的对称轴上的结论。本发明将该方法进行了一些修正，描述如下。

假设X(u)为u域上的一个连续谱函数，若该谱函数在某一支撑区[a，b]内是关于谱峰点u＝u₀对称的，如图1所示，则显然X(u)在该支撑区内的能量重心也位于该峰值点上。用公式表示为

$u_0=\frac{{\∫}_a^b\left|X\left(u\right)\right|\×\mathrm{udu}}{{\∫}_a^b{\left|X\left(u\right)\right|}^2\mathrm{du}}---\left(15\right)$

离散情况下，对X(u)进行等间隔采样，采样间隔为Δu，得到

X(u_i)|_u＝iΔu＝X(u)|_u＝iΔu，i∈Z                (16)

假设序号某一序号m∈z的对应的谱线X(u_m)＝X(mΔu)处于峰值谱线的位置。则此时能量重心位置可以由所有谱线能量的位置加权和与总的谱线能量之比得到，实际应用中，在谱峰能量比较集中的情况下，可以在峰值谱线两边各取能量较大的n根谱线做近似计算，即

$\widehat{u_0}=\frac{\underset{i=m-n}{\overset{m+n}{\mathrm{\Σ}}}{\left|X\left(u_i\right)\right|}^2\·u_i}{\underset{i=m-n}{\overset{m+n}{\mathrm{\Σ}}}{\left|X\left(u_i\right)\right|}^2},$ u_i＝iΔu               (17)

式(17)说明，若一个连续谱在某一支撑区内关于谱峰对称，则可以通过其离散谱线的能量重心位置来恢复真实的谱峰位置。以上就是离散谱校正的能量重心法原理。

3、基于RAT和FRFT的能量重心插值估计的理论基础

根据前面的结论，我们将证明有限长LFM信号的RAT和FRFT连续谱分别关于其谱峰点对称，从而可以通过计算两个离散谱的能量重心位置来恢复真实的谱峰位置，实现谱峰位置的超分辨率估计。具体的证明过程如下：

首先证明RAT谱η(k)关于峰值点k₀的对称性，给出时间长度为T的单位能量LFM信号表达式

$s\left(t\right)=\frac{1}{\sqrt{T}}\exp \left\{\mathrm{j\π}(2f_{0}t+{\mathrm{\μ}}_0{t}^2)\right\},$ $-\frac{T}{2}\≤t\≤\frac{T}{2}---\left(18\right)$

取该信号的AF并取模

$\left|{\mathrm{AF}}_s(\mathrm{\τ},\mathrm{\ξ})\right|=\left|\frac{\sin \left[\mathrm{\π}(T-|\mathrm{\τ}\left|\right)(\mathrm{\ξ}-{\mathrm{\μ}}_0\mathrm{\τ})\right]}{\mathrm{\πT}(\mathrm{\ξ}-{\mathrm{\μ}}_0\mathrm{\τ})}\right|,$ -T≤τ≤T(19)

对上式做直线ξ＝kτ上的Radon变换

$\mathrm{\η}\left(k\right)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}\left|{\mathrm{AF}}_s(\mathrm{\τ},\mathrm{\ξ})\right|\mathrm{\δ}(\mathrm{\ξ}-\mathrm{kt})\mathrm{d\ξd\τ}$

$={\∫}_{-T}^T\left|\frac{\sin \left[\mathrm{\π}(T-|\mathrm{\τ}\left|\right)(k-{\mathrm{\μ}}_0)\mathrm{\τ}\right]}{\mathrm{\πT}\left({k-\mathrm{\μ}}_0\right)\mathrm{\τ}}\right|\mathrm{d\τ}---\left(20\right)$

考察函数g(ε)＝η(k+ε)

$g\left(\mathrm{\ϵ}\right)={\∫}_{-T}^T\left|\frac{\sin \left[\mathrm{\π}(T-|\mathrm{\τ}\left|\right)(k+\mathrm{\ϵ}-{\mathrm{\μ}}_0)\mathrm{\τ}\right]}{\mathrm{\πT}(k+\mathrm{\ϵ}-{\mathrm{\μ}}_0)\mathrm{\τ}}\right|\mathrm{d\τ}---\left(21\right)$

将k＝k₀＝μ₀代入上式并化简，得

$g\left(\mathrm{\ϵ}\right)={\∫}_{-T}^T\left|\frac{\sin \left[\mathrm{\π}(T-|\mathrm{\τ}\left|\right)\mathrm{\ϵ\τ}\right]}{\mathrm{\πT\ϵ\τ}}\right|\mathrm{d\τ}---\left(22\right)$

显然式(22)是关于ε的偶函数，则可以说明η(k)关于轴k＝μ₀对称。

另外证明FRFT谱|S_α(u)|关于峰值点u₀的对称性，假设由RAT方法得到了如式(18)所示的信号调频率μ₀的无偏估计，对该信号做式(12)所示的能量最聚集域上的FRFT

$S_{{\mathrm{\α}}_0}\left(u\right)={\∫}_{-\∞}^{\∞}s\left(t\right)K_{{\mathrm{\α}}_0}(t,u)\mathrm{dt}$

$={\∫}_{-T/2}^{T/2}\frac{1}{\sqrt{T}}\exp \left\{\mathrm{j\π}(2f_{0}t+{\mathrm{\μ}}_0{t}^2)\right\}\sqrt{1-j{\cot \mathrm{\α}}_0}\exp \left\{\mathrm{j\π}\right[t^2\cot {\mathrm{\α}}_0-2\mathrm{tu}\csc {\mathrm{\α}}_0+u^2\cot {\mathrm{\α}}_0\left]\right\}\mathrm{dt}$

$=\sqrt{\frac{1-j\cot {\mathrm{\α}}_0}{T}}\exp \left\{\mathrm{j\π}{u}^2\cot {\mathrm{\α}}_0\right\}\frac{\sin \left[\mathrm{\π}(f_0-u\csc {\mathrm{\α}}_0)T\right]}{\mathrm{\π}(f_0-u\csc {\mathrm{\α}}_0)}$

                                    (23)

对上式取模

$\left|S_{{\mathrm{\α}}_0}\left(u\right)\right|=\left|\sqrt{\frac{1-j\cot {\mathrm{\α}}_0}{T}}\exp \right\{j{\mathrm{\πu}}^2\cot {\mathrm{\α}}_0\left\}\frac{\sin \left[\mathrm{\π}(f_0-u\csc {\mathrm{\α}}_0)T\right]}{\mathrm{\π}(f_0-u\csc {\mathrm{\α}}_0)}\right|$

$=\sqrt[4]{\frac{1+{\left(\cot {\mathrm{\α}}_0\right)}^2}{T^2}}T\left|\sin c\right[\mathrm{\π}(f_0-u\csc {\mathrm{\α}}_0)T\left]\right|---\left(24\right)$

考察函数 $h\left(\mathrm{\ϵ}\right)=\left|S_{{\mathrm{\α}}_0}(u+\mathrm{\ϵ})\right|$

$h\left(\mathrm{\ϵ}\right)=\sqrt[4]{\frac{1+{\left(\cot {\mathrm{\α}}_0\right)}^2}{T^2}}T\left|\sin c\right\{\mathrm{\π}[f_0-(u+\mathrm{\ϵ})\csc {\mathrm{\α}}_0]T\left\}\right|---\left(25\right)$

将u＝f₀/cscα₀代入上式并化简，得

$h\left(\mathrm{\ϵ}\right)=\sqrt[4]{\frac{1+{\left(\cot {\mathrm{\α}}_0\right)}^2}{T^2}}T\left|\sin c\right\{-\mathrm{\π\ϵ}\csc {\mathrm{\α}}_{0}T\left\}\right|---\left(26\right)$

显然式(26)是关于ε的偶函数，则可以说明

关于轴u＝u₀＝f₀/cscα₀对称，证明结束。

4、具体的调频率和中心频率的估计流程

综上所述，关于chirp信号的调频率参数估计的流程图如图2所示。

给定输入信号s(n)，首先求取其模糊函数AF_s(τ，ξ)，然后不断更换扫描直线ξ＝kτ，借助Radon变换计算得出AF_s(τ，ξ)在扫描直线ξ＝kτ上的投影，并记录所有投影值，构造全景离散谱，接着搜索出该离散谱的谱峰位置，并借助式(17)所示的能量重心插值公式，估计实际峰值投影值，最后根据式(11)得出信号调频率的估计。

得出信号调频率估计值

后，其中心频率

的估计流程如图3所示。

图3中，借助本专利提出的基于能量重心插值估计和RAT的调频率估计方法，用得到的调频率估计值

对信号做分数阶傅立叶变换，进而搜索出分数阶傅立叶谱的峰值谱，借助式(17)计算出谱能量重心位置，再借助

即可得到信号重心频率的估计。

下面首先对实施本发明的硬件予以简单说明。参见图4，为精确估计出LFM信号x(t)的调频率和中心频率参数，需借助信号调理电路对输入信号进行模拟预处理，以对信号幅度范围进行必要调整，并去除外干扰噪声等；再经过A/D(模数转化器)采样得到样本序列x(n)以并行数字输入的形式进入DSP器件，经过DSP器件的内部的算法处理，而得到信号参数的估计，最后借助输出驱动及其显示模块显示出调频率和中心频率的估计值，即图4的整个系统构成了一个“高精度LFM信号参数估计仪”。

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

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

(2)根据实际需要调整采样率f_s，使得在该采样率下，尽量高精度地估计出信号参数；

(3)内部RAM存储数据不足时，将处理数据与外部RAM进行数据交换，以配合核心算法处理；

(4)将相位估计结果实时输出到驱动和显示模块。

本发明将所提出的“基于RAT和FRFT能量重心插值估计”的核心估计算法植入DSP器件内，基于此完成高精度、高速、高效的LFM信号参数估计。需指出，由于采用了数字化的估计方法，因而决定图4系统的复杂度、实时程度和稳定度的主要因素并不是图4中的DSP器件的外围连接，而是DSP内部程序存储器所存储的核心估计算法。DSP器件的内部程序流程如图5所示。

为了从数值上说明本发明所提出的LFM信号参数估计方法的有效性，构造一个单分量LFM信号，参数选取：幅度a₀＝1，中心频率f₀＝10.725Hz，调频率μ₀＝1.858Hz/s，采样点数N＝501，采样频率f_s＝1000Hz。在10dB高斯白噪声环境中，分别考察传统扫描方法和本发明方法的调频率和中心频率的估计精度。图6为该观测信号的模糊函数。相应的信号RAT谱和能量最聚集域的FRFT谱如图7、图8所示。

按照前文所述步骤对该信号的参数进行估计，重复进行100次Monte Carlo仿真，统计结果如表1所示。

表1单分量情况的仿真结果

另外构造一个由两个LFM分量组成的复合信号，参数选取：幅度a₀＝1.5，a₁＝1；中心频率f₀＝6.352Hz，f₁＝12.885Hz；调频率μ₀＝2.145Hz/s，μ₁＝1.293Hz/s；采样点数N＝501，采样频率f_s＝1000Hz。在10dB高斯白噪声环境中，分别考察传统扫描方法和本发明方法的调频率和中心频率的估计精度。按照前文所述步骤对该复合信号的参数进行估计，重复进行100次Monte Carlo仿真，统计结果如表2所示。

表2多分量情况的仿真结果

表1实验数据表明，在单分量情况下，特别是在真实谱峰位置未对准扫描步长的整数倍位置时，本文方法的估计精度远优于传统方法，对信号的RAT和FRFT谱峰位置(μ₀和u₀)的估计误差分别达到了扫描步长的1.55％和4.94％；对于两个分量的情况，由于分量之间的交叉项干扰，使得估计精度略有下降，但本文方法仍然优于传统的扫描方法。

图3流程中各处理阶段所涉及的具体硬件资源耗费情况，分别作如下阐述：

(1)对输入信号x(t)进行采样而得到离散数据x(n)

这里对采样频率f_s值并不需要作严格限制，不要求信号频率为频率分辨率(由采样值f_s除以信号长度N而决定)的整数倍。另外，对采样器的转换精度也没有提出过高要求，选用常用的字长数大于8bit以上的采样器即可。事实上，一般DSP器件所配备的最低的A/D转换精度就很高(如TMS320LF2407与TMS320F2812内部自带的A/D转换器为10bit的转换精度)，因而这足以保证本发明的算法在硬件实现时具有足够高的参数估计精度。

(2)内部数据的存储

可把LFM信号参数估计的程序存放在DSP器件的程序存储器内。由上面的描述可知，本发明提出的估计算法需要N个输入数据进行存储处理。需指出的是，通过能量重心插值措施，在不增加采样频率和样点个数的前提下保证了很高的参数估计精度；反过来考虑，这意味着在同样的估计精度的前提下，本发明提出的方法所需的样点数比传统方法要少得多。通过上面的仿真实验可知，该算法在采样点数N＝501时，已经具有了很高的估计精度，就算以16位的采样精度进行采样，在每个数据占用2个字节(2Bytes)的存储空间，这意味着仅需1002个字节(不到1kBytes)的内部RAM空间就可获得很高的估计精度，这对通用DSP芯片的内部RAM来说(如TMS320c2407的内部RAM为5kBytes，TMS320c2812的内部RAM为20kBytes，TMS320c5407的内部RAM为40kBytes，TMS320c6x的内部RAM为256kBytes)，足以存下这些数据；另外，该算法的主要处理是在谱生成和谱峰搜索环节，而在处理过程中，无需开辟新的数据存储空间进行算法处理。

(3)算法程序调用

而就算法程序本身而言，包括以下部分：①数据预处理过程；②RAT变换处理过程；③RAT谱峰的搜索和校正过程；④FRFT变换处理过程；⑤FRFT谱峰的搜索和校正过程。共5个步骤组成，程序简单。另外，在RAT变换的数据运算时，利用了模糊函数和分数阶相关的关系来对运算过程进行简化，即采用FRFT以及Fourier变换的运算来实现RAT变换的运算，而FRFT和Fourier变换具有快速算法，运算复杂度仅为O(NlgN)。对取得的谱线做峰值搜索时，只需记录其序号即可，然后取峰值谱线两旁的能量较为集中的几根谱线参与校正，这部分的运算复杂度很小，可以忽略。总而言之，整个LFM信号参数估计程序非常简单，无需外部扩展程序存储器进行程序调用。

(4)计算结果输出

在由DSP硬件计算得到LFM信号的调频率和中心频率的估计值后，直接可通过DSP的输出总线输出至外部显示驱动设备进行数码显示。

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

参考文献：

[1]L.科恩.时频分析：理论与应用，白居宪译，西安：西安交通大学出版社。

[2]B.Boashash.Estimating and Interpreting the Instantaneous Frequency of a Signal.Part 2：Algorithms and Application，Proceeding of the IEEE，80(4)，1992：540-568。

[3]B.Boashash，G.Jones，P.O.Shea.Instantaneous Frequency of Signals：Concepts，Estimation Techniques and Applications.Proc.of SPIE，vol.1152，1989：382-400。

[4]B.Boashash，P.O.Shea，M.Arnold.Algorithms for Instantaneous Frequency Estimation：AComparative Study.Proc.of SPIE，vol.1348，1990。

[5]B.Ferguson.A ground-based narrow-band passive acoustic technique for estimating thealtitude and speed of a propeller-driven aircraft.Journal of the Acoustical Society of America，92(3)，1992：1403-1407。

[6]D.L.Snyder.A State Space Approach to Analog Communication System.Cambridge，MA：MIT Press。

[7]L.Griffiths.Rapid Measurement of Digital Instantaneous Frequency.IEEE Trans.on ASSP，23(2)，1975：207-222。

[8]S.Haykin.Adaptive Filter Theory.Englwood Cliffs，NJ：Prentice Hall，1991 2nded。

[9]T.A.C.M.Classen，W.F.G.Mecklenbrauker.The Wigner Distribution.Part II.PhilipsJournal of Research，Vol.35，1980：276-300。

[10]L.B.White，B.Boashash.Estimating the Instantaneous Frequency of a Gaussian RandomProcess.IEEE Trans.on ASSP，Vol.36，1988：1518-1521。

[11]A.M.Haimovich，et al.SAR Imagery of Moving targets：Application of Time-FrequencyDistributions for Estimating Motion Parameters.Proc.1994 SPIE′s International Symposium onAerospace and Sensing，1994，Vol.2238：238-247。

[12]B.Boashash.Representation of temps-frequency.Ph.D Thesis，Univ.Grenoble，France，1982。

[13]Rao，F.Taylor.Estimation of Instantaneous Frequency Using the Discrete WignerDistribution.Electron.Lett，26(4)，1990：246-248。

[14]B.Boashash，P.O.Shea.Time-varying High Order Spectra.Proc.SPIE，Conf.AdvancedAlgorithms and Architectures of Signal Processing VI，1991。

[15]S.Barbarossa，V.Petrone.Analysis of Polynomial-Phase Signals by the IntegratedGeneralized Ambiguity Function.IEEE Trans.on SP，1997，45(2)：316-327。

[16]S.Peleg，B.Porat.Estimation and Classification of Polynomial-Phase Signal.IEEE Trans.on Info.Theory，37(3)，1991：422-430。

[17]H.Choi，W.J.Williams.Improved Time-frequency Representation of MulticomponentSignals Using Exponential Kernels.IEEE Trans.on SP，1988，37(6)：862-871。

[18]S.Barbarossa.Analysis of Multicomponent LFM Signals by a Combined Wigner-HoughTransform.IEEE Trans.on SP，1995，43(6)：1511-1515。

[19]田孝华，廖桂生，吴云韬.LFM脉冲雷达回波Doppler与多径时延的联合估计.电子学报，2002，30(6)：857-860。

[20]T.J.Abatzoglou.Fast maximum likelihood joint estimation of frequency and frequencyrate.IEEE Trans.on AES，1986，22(6)：708-715。

[21]R.Kuaresen，S.Verma.On Estimating the Parameters of Chirp Using Rank ReductionTechniques.Proc.21st Asilomar Conference on Signal，Systems and Computers，1987：555-558。

[22]陶然，邓兵，王越.分数阶傅里叶变换及其应用[M].北京：清华大学出版社，2009.9。

[23]L.B.Almeida.The fractional Fourier transform and time-frequency representations[J].IEEE Trans on SP，1994，42(11)：3084～3091。

[24]H.M.Ozaktas，Orhan Arikan.Digital Computation of the Fractional Fourier Transform[J]IEEE Trans on SP，1996，44(9)：2141～2150。

[25]齐林，陶然，周思永，王越.基于分数阶Fourier变换的多分量LFM信号的检测和参数估计[J].中国科学(E辑)，2003，33(8)：749～759。

[26]Akay O.Fractional convolution and correlation via operator methods and an application todetection of linear FM signals[J].IEEE Trans on SP，2001，49(5)：979-993。

[27]Qi Lin，Tao Ran，Zhou Siyong，et al.Detection and parameter estimation ofmulticomponent LFM signal based on the fractional Fourier transform[J].Science in China，Ser.F，2004，47(2)：184-198。

[28]丁康，江利旗.离散频谱的能量重心校正法[J].振动工程学报，2001，14(3)：354～359。

Claims (3)

1.一种线性调频信号参数估计方法，其特征是，包括下列步骤：

借助基于能量重心插值估计和Radon模糊变换的调频率估计方法，具体为：对于给定输入信号s(n)，首先计算s(n)的模糊函数AF_s(τ，ξ)；然后从初始调频率μ=0开始以指定的Δμ为调频率扫描步长逐次等间隔地在模糊域(τ,ξ)平面内构造过原点的扫描直线,则第k次扫描直线可表示为ξ＝(kΔμ)τ，τ、ξ分别为模糊参数，再将模糊函数AF_s(τ，ξ)在该直线上进行投影得到投影值η(k)，并记录所有投影值η(k)以构造全景离散谱，接着搜索该离散谱峰位置

记录峰值谱k＝k^*处及其前后各p个投影值，再借助能量重心插值公式

μ_k＝k^*+kΔμ，即得出信号调频率的估计值

然后，用得到的调频率估计值计算出分数阶傅立叶变换的角度

以α₀作为参数对信号做分数阶傅立叶变换，而得到分数阶傅立叶变换的离散谱S(u_i),其中u_i=iΔu，Δu为离散谱间隔，i为分数阶傅立叶域的谱序号，进而找到其峰值谱位置u=mΔu，记录峰值谱u=mΔu处及其前后各q个投影值，借助

u_i=iΔu计算出谱能量重心位置

再借助即可得到信号中心频率的估计。

2.一种线性调频信号参数估计实施装置，其特征是，包括：

信号调理电路，用于对输入信号进行模拟预处理，以对信号幅度范围进行必要调整，并去除外干扰噪声；

A/D模数转化器，用于采样得到样本序列s(n)，以并行数字输入的形式送入数字信号处理器；经过数字信号处理器的内部的算法处理，而得到信号参数的估计，

数字信号处理器，用于：

(1)调用核心算法，其为如权利要求1所述的线性调频信号参数估计方法，完成接收信号的参数估计处理；

(2)根据实际需要调整采样率f_s，使得在该采样率下，尽量高精度地估计出信号参数；

(3)内部RAM存储数据不足时，将处理数据与外部RAM进行数据交换，以配合核心算法处理；

(4)将相位估计结果实时输出到显示模块及其输出驱动，用于显示出调频率和中心频率的估计值。

3.根据权利要求2所述的一种线性调频信号参数估计实施装置，其特征是，数字信号处理器进一步包括：①数据预处理模块；②RAT变换处理模块；③RAT谱峰的搜索和校正模块；

④分数阶Fourier变换FRFT变换处理模块；⑤分数阶Fourier变换FRFT谱峰的搜索和校正模块。

Images