CN104360322B - 基于四阶非对称乘积型核函数的qfm信号参数估计方法 - Google Patents

Abstract

本发明属于雷达信号检测与估计技术领域，特别涉及基于四阶非对称乘积型核函数的QFM信号参数估计方法，其步骤为：计算初始信号能量，设定循环初始值以及建立存放后续各分量参数估计值的矩阵；将函数的时延函数与其共轭时延函数相乘，再对其相位匹配变换，在频域峰值处可得到信号分量第二个参数的估计值；对补偿第二个参数估计值后的信号构造新的四阶非对称乘积型核函数再做快速傅里叶变换得到频率‑调频率变化率的二维分布图，从峰值处可得到信号分量的第一个和第三个参数估计值；步骤4，将通过上述步骤估计出的信号分量从信号中滤除；步骤5，重置循环变量和信号，计算信号总能量与初始信号能量比值。

Description

基于四阶非对称乘积型核函数的QFM信号参数估计方法

技术领域

本发明属于雷达信号检测与估计技术领域，特别涉及基于四阶非对称乘积型核函数的QFM信号参数估计方法。

背景技术

调频信号是一种瞬时频率随着时间而变化的信号，因此是典型的非平稳信号。根据调频率与时间的线性关系，调频信号可以分为线性调频和非线性调频两类信号。瞬时频率与时间具有二次函数关系的QFM(Quadratic frequency modulated,二次调频)信号是一种普遍存在于自然界和人工应用中的信号，在声呐、高分辨雷达、地震信号分析、遥感遥测、通信等工程和科学领域有着广泛的应用。

时频分析方法能够同时描述信号在不同时间不同频率的能量强度或密度分布，是分析非平稳信号的一种直接和有效的方法。时频分析方法可以分为线性变换和非线性变换两类，典型的线性变换有短时傅立叶变换(Short time Fourier Transform,STFT)、Gabor变换、小波变换等。线性变换由于受到测不准原理的限制，不可能同时得到很高的时间分辨率和频率分辨率。典型的非线性变换方法有：Wigner-Ville分布(WVD)、模糊函数、Cohen类时频分布等双线性时频分布和多项式Wigner-Ville分布(Polynomial WVD，PWVD)等高阶时频分布。双线性时频分布能够较好地分析线性调频(Linear frequency modulated,LFM)信号，但在分析QFM信号时会产生严重的交叉项(包括单分量时的内部交叉项和多分量时的外部交叉项)，从而影响分析结果。多项式Wigner-Ville分布通过选取合适的参数能够对单分量的QFM信号具有最佳的时频聚集性，但对于多分量的信号会产生严重的外部交叉项。时频分布是非参数化的分析方法，直接由时频分布可以描述信号的瞬时频率信息，但不能完成对信号参数的直接估计以及信号的重构。对于QFM信号，最大似然估计(MaximumLikelihood，ML)无疑是具有最优检测和估计性能的参数估计方法，但它在实际应用中需要进行3维的最大值搜索，运算量极大。为了避免3维目标函数和相应的3维的最大值搜索，三次相位函数(Cubic Phase Function，CPF)被提出并运用，但对于多分量的QFM信号却无法处理；PHMT(Product High-Order Matched-Phase Transform)和PGCPF(ProductGeneralized Cubic Phase Function)方法可以应用于多分量QFM信号的参数估计，但这两种方法经过高阶非线性变换后均只能估计最高阶的参数，无法同时估计其他参数，需要通过迭代补偿不断搜索估计不同阶数的参数，补偿精度影响估计效果。

发明内容

本发明的目的在于提出基于四阶非对称乘积型核函数的QFM信号参数估计方法，本发明先估计信号的调频率，再同时估计信号的调频率变化率和中心频率，最后估计信号的幅度，从而重构原信号。

为实现上述技术目的，本发明采用如下技术方案予以实现。

基于四阶非对称乘积型核函数的QFM信号参数估计方法包括以下步骤：

步骤1，获取初始多分量二次调频信号x(n)，n表示采样序列时刻点，N表示总采样点数；得出初始多分量二次调频信号的总能量E₀；设置循环次数变量r，r＝0,1,2,...，当r＝0时，执行步骤2；

步骤2，构造信号x(n)的时延函数x(n+m)，m表示信号x(n)时延的时刻点数，得出函数y(n,m)，y(n,m)＝x(n+m)x(-n-m)；将信号x(n)中具有最大幅度值的分量记为信号x(n)的第j个分量，根据函数y(n,m)，得出信号x(n)中具有最大幅度值的分量的调频率估计值

步骤3，利用信号x(n)中具有最大幅度值的分量的调频率估计值对信号x(n)进行补偿，得到补偿后信号x_d1(n)；利用补偿后信号x_d1(n)，构造四阶非对称乘积型核函数y₁(n,m)：

$y_1(n,m)=x_{d1}(n+\frac{2}{3}m)x_{d1}(n+m){x_{d1}}^{*}(n+\frac{2+\sqrt{22}}{6}m){x_{d1}}^{*}(n+\frac{2-\sqrt{22}}{6}m)$

其中，上标*表示取共轭；

对四阶非对称乘积型核函数y₁(n,m)进行相位匹配变换，得到峰值函数Y₁(f₃,m)，f₃表示对应于时刻n的频域频率变量；对峰值函数Y₁(f₃,m)进行由m时域至f₁频域的快速傅里叶变换，得到峰值函数Y₂(f₃,f₁)，f₁表示对应于时延m的频域频率变化率；将峰值函数Y₂(f₃,f₁)的频谱图的最大幅度值在f₁轴上的投影作为信号x(n)中具有最大幅度值的分量的中心频率的估计值将峰值函数Y₂(f₃,f₁)的频谱图的最大幅度值在f₃轴上的投影作为信号x(n)中具有最大幅度值的分量的调频率变化率的估计值

步骤4，得出信号x(n)中第j个分量的幅度值的估计值利用信号x(n)中第j个分量的幅度值的估计值、中心频率的估计值、调频率的估计值、调频率变化率的估计值，重构出信号x(n)中第j个分量s_j(n)；在信号x(n)中滤除第j个分量s_j(n)，得出信号x_r(n)；

步骤5，得出信号x_r(n)的总能量E，计算出信号x_r(n)的总能量E和初始多分量二次调频信号的总能量E₀的比值ε，如果ε大于或等于设定的门限值，则令x(n)＝x_r(n)，并令r的值自增1，返回至步骤2；否则，说明二次调频信号参数估计过程完成。

本发明的有益效果为：

1)本发明采用的四阶非对称乘积型核函数，与常规的核函数(非线性阶数过高的核函数，或者虽然含有较低的非线性但却含有复数时延的核函数)相比，不会产生由于过多的外部交叉项影响积分结果，也不会产生由于复数时延对信号的频率和幅度调制而导致积分后的最大值位置与信号的自身参数不对应的现象，很好的解决过高的非线性阶数和复数时延的问题。

2)现有的如PHMT(Product High order Matched-phase Transform，高次乘积相位匹配变换)和PGCPF(Product Generalized Cubic Phase Function，三次乘积型相位核函数)方法都是顺序估计信号的三次项，二次项和一次项系数。由于低次项系数的估计是在高次项估计并补偿之后进行的，因次高次项系数的估计精度会影响低次项的估计精度，并且三次项系数的估计值是在构造高阶核函数(PHMT构造6阶核函数，PGCPF构造4阶核函数)得到的，信噪比门限会比较高。与PHMT和PGCPF方法不同，本发明先通过构造2阶核函数估计二次项的系数，信噪比门限更低，在多分量信号情况下可以更好的检测并估计信号的参数，本发明中新的乘积型的变换核函数保证了三次项的系数和一次项的系数分别沿时间轴和时延轴彼此独立，因此可以同时估计三次项和一次项的系数，这样可以避免高次项的估计精度对低次项的估计精度的影响。

附图说明

图1为本发明的基于四阶非对称乘积型核函数的QFM信号参数估计方法的流程图；

图2为仿真实验1中初始多分量二次调频信号x(n)的两个分量在调频率变化率-频率平面的分布示意图；

图3为仿真实验1中初始多分量二次调频信号x(n)的两个分量在调频率变化率-频率-信号幅度值三维空间的分布示意图；

图4为仿真实验2中初始多分量二次调频信号中加入0dB的高斯白噪声时得出的初始多分量二次调频信号的各分量在调频率变化率-频率-信号幅度值三维空间的分布示意图；

图5为仿真实验3中采用PHMT方法得出的初始多分量二次调频信号各分量的峰值位置的估计结果示意图；

图6为仿真实验4中采用PHMT方法得出的初始多分量二次调频信号各分量的峰值位置的估计结果示意图；

图7，为仿真实验4中采用PGCPF方法得出的初始多分量二次调频信号各分量在调频率变化率-频率-信号幅度值三维空间的分布示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的基于四阶非对称乘积型核函数的QFM信号参数估计方法的流程图。该基于四阶非对称乘积型核函数的QFM信号参数估计方法包括以下步骤：

步骤1，获取初始多分量二次调频信号x(n)，n表示采样序列时刻点，N表示总采样点数；得出初始多分量二次调频信号的总能量E₀；设置循环次数变量r，r＝0,1,2,...，当r＝0时，执行步骤2；

其具体步骤为：

获取初始多分量二次调频信号(QFM信号)x(n)，初始多分量二次调频信号x(n)的表达式为：

$x\left(n\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\σ}}_i\exp \left(j2\mathrm{\π}(a_{i,1}n+\frac{1}{2}{a}_{i,2}{n}^2+\frac{1}{6}{a}_{i,3}{n}^3)\right),n\∈[-\frac{N-1}{2},\frac{N-1}{2}]$

其中，n表示采样序列时刻点，N表示总采样点数，N为自然数，i取1至k，k表示初始多分量二次调频信号的分量个数；σ_i表示初始多分量二次调频信号的第i个分量的幅度值，a_i,1表示初始多分量二次调频信号的第i个分量的中心频率，a_i,2表示初始多分量二次调频信号的第i个分量的调频率，a_i,3表示初始多分量二次调频信号的第i个分量的调频率变化率。

得出初始多分量二次调频信号x(n)的总能量E₀：

$E_0=\underset{n=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\mathrm{\Σ}}}{\left|x\left(n\right)\right|}^2.$

其中，|x(n)|表示初始多分量二次调频信号x(n)的模值。

构建大小为4×k的矩阵signal，矩阵signal每个元素的初始值为任意值。

设置循环次数变量(估计出的信号分量的个数)r，r＝0,1,2,...，当r＝0时，执行步骤2。

步骤2，构造信号x(n)的时延函数x(n+m)，m表示信号x(n)时延的时刻点数，得出函数y(n,m)，y(n,m)＝x(n+m)x(-n-m)；将信号x(n)中具有最大幅度值的分量记为信号x(n)的第j个分量，根据函数y(n,m)，得出信号x(n)中具有最大幅度值的分量的调频率估计值

其具体子步骤为：

(2.1)构造信号x(n)的时延函数x(n+m)，m表示信号x(n)时延的时刻点数，m为自然数。得出函数y(n,m)，

$y(n,m)=x(n+m)x(-n-m)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\σ}}_i\exp \left[j2\mathrm{\π}{a}_{i,2}{(n+m)}^2\right]+B(n,m).$

其中，B(n,m)表示信号分量相乘的交叉项。

(2.2)将信号x(n)中具有最大幅度值的分量记为信号x(n)的第j个分量。根据函数y(n,m)，得出信号x(n)中具有最大幅度值的分量的调频率估计值

${\hat{a}}_{j,2}=\underset{f_2}{\arg \max }\left|\mathrm{PY}\left(f_2\right)\right|$

$\mathrm{PY}\left(f_2\right)=\underset{n_1=-\frac{N}{4}}{\overset{\frac{N}{4}}{\mathrm{\Π}}}Y(n_1,f_2)$

$\begin{array}{c}Y(n,f_2)=\underset{m=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\mathrm{\Σ}}}y(n,m)e^{-j2\mathrm{\π}{f}_2(m^2+2\mathrm{nm})}\\ =\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\σ}}_i\exp \left(j2\mathrm{\π}{a}_{i,2}{n}^2\right)\·\mathrm{\δ}(f_2-a_{i,2}){+B}^{\′}(n,f_2)\end{array}$

其中，f₂表示时延m对应的频域频率变量，|PY(f₂)|表示PY(f₂)的模值，

$n_1\∈[-\frac{N}{4},\frac{N}{4}],n\∈[-\frac{N-1}{2},\frac{N-1}{2}],m\∈[-\frac{N-1}{2},\frac{N-1}{2}],$ B′(n,f₂)为经过B(n,m)相位匹配变换后的交叉项，δ(f₂-a_i,2)为在f₂频域a_i,2位置处输出尖峰的δ函数。Y(n,f₂)为信号x(n)各分量相乘的交叉项在调频率变量处的峰值函数，由于m从至将函数值累加，因此该函数为时间-频率函数。PY(f₂)为累积相乘后函数，本发明实施例中，通过累积相乘后函数PY(f₂)的频谱图幅值的最大值在频率f₂轴上的投影可以得到信号x(n)中具有最大幅度值的分量的调频率估计值

步骤3，利用信号x(n)中具有最大幅度值的分量的调频率估计值对信号x(n)进行补偿，得到补偿后信号x_d1(n)；利用补偿后信号x_d1(n)，构造四阶非对称乘积型核函数y₁(n,m)：

$y_1(n,m)=x_{d1}(n+\frac{2}{3}m)x_{d1}(n+m){x_{d1}}^{*}(n+\frac{2+\sqrt{22}}{6}m){x_{d1}}^{*}(n+\frac{2-\sqrt{22}}{6}m)$

其中， $m\∈[-\frac{N-1}{2},\frac{N-1}{2}],$ 上标*表示取共轭；

对四阶非对称乘积型核函数y₁(n,m)进行相位匹配变换，得到峰值函数Y₁(f₃,m)，f₃表示对应于时刻n的频域频率变量；对峰值函数Y₁(f₃,m)进行由m时域至f₁频域的快速傅里叶变换，得到峰值函数Y₂(f₃,f₁)，f₁表示对应于时延m的频域频率变化率；将峰值函数Y₂(f₃,f₁)的频谱图的最大幅度值在f₁轴上的投影作为信号x(n)中具有最大幅度值的分量的中心频率的估计值将峰值函数Y₂(f₃,f₁)的频谱图的最大幅度值在f₃轴上的投影作为信号x(n)中具有最大幅度值的分量的调频率变化率的估计值

其具体子步骤为：

(3.1)利用信号x(n)中具有最大幅度值的分量的调频率估计值对信号x(n)进行补偿，得到补偿后信号x_d1(n)； $n\∈[-\frac{N-1}{2},\frac{N-1}{2}],$ N表示总采样点数。

利用补偿后信号x_d1(n)，构造四阶非对称乘积型核函数y₁(n,m)：

$y_1(n,m)=x_{d1}(n+\frac{2}{3}m)x_{d1}(n+m){x_{d1}}^{*}(n+\frac{2+\sqrt{22}}{6}m){x_{d1}}^{*}(n+\frac{2-\sqrt{22}}{6}m)$

其中，上标*表示取共轭。经过推导可以得出：

$y_1(n,m)={\mathrm{\σ}}_j\exp \left(j2\mathrm{\π}(a_{j,1}+\frac{1}{2}{a}_{j,3}{\mathrm{mn}}^2)\right)+A(n,m),$

其中，σ_j表示信号x(n)中第j个分量的幅度值，a_j,1表示信号x(n)中具有最大幅度值的分量的中心频率，a_j,3表示信号x(n)中具有最大幅度值的分量的调频率变化率。A(n,m)为含有信号各分量相乘的交叉项以及除能量最强信号分量之外所有信号各分量自身项的函数。

(3.2)对四阶非对称乘积型核函数y₁(n,m)进行相位匹配变换，得到峰值函数Y₁(f₃,m)，

$\begin{array}{c}{Y}_1(f_3,m)=\underset{n=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\mathrm{\Σ}}}{y}_1(n,m)e^{-j2\mathrm{\π}{f}_3{n}^{2}m}\left|n\right|\\ ={\mathrm{\σ}}_j\exp \left(j2\mathrm{\π}{a}_{j,1}m\right)\·\mathrm{\δ}(f_3-\frac{1}{2}{a}_{j,3})+A^{\′}(f_3,m)\end{array}$

其中，f₃表示对应于时刻n的频域频率变量，|n|表示n的绝对值，为在f₃频域a_j,3位置处输出尖峰的δ函数，A′(f₃,m)为A(n,m)的相位匹配变换后的交叉项以及除最强分量外所有信号各分量自身项函数。

对峰值函数Y₁(f₃,m)进行由m时域至f₁频域的快速傅里叶变换，得到在位置a_j,1和位置a_j,3同时输出峰值的峰值函数Y₂(f₃,f₁)，f₁表示对应于时延m的频域频率变化率；

$\begin{array}{c}{Y}_2\left(f_3{,f}_1\right)=\underset{m\→f_1}{\mathrm{FFT}}\left[Y_1(f_3,m)\right]\\ ={\mathrm{\σ}}_j\mathrm{\δ}(f_1-a_{j,1})\·\mathrm{\δ}(f_3-\frac{1}{2}{a}_{j,3})+A^{\′\′}(f_3,f_1)\end{array}$

其中，δ(f₁-a_j,1)为在f₁频域a_j,1位置处输出尖峰的δ函数，A″(f₃,f₁)表示A′(f₃,m)由m时域至f₁频域的快速傅里叶变换。

得出信号x(n)中具有最大幅度值的分量的中心频率的估计值和信号x(n)中具有最大幅度值的分量的调频率变化率的估计值信号x(n)中具有最大幅度值的分量的中心频率的估计值为：峰值函数Y₂(f₃,f₁)的频谱图的最大幅度值在f₁轴上的投影，信号x(n)中具有最大幅度值的分量的调频率变化率的估计值为：峰值函数Y₂(f₃,f₁)的频谱图的最大幅度值在f₃轴上的投影。

步骤4，得出信号x(n)中第j个分量的幅度值的估计值利用信号x(n)中第j个分量的幅度值的估计值、中心频率的估计值、调频率的估计值、调频率变化率的估计值，重构出信号x(n)中第j个分量s_j(n)；在信号x(n)中滤除第j个分量s_j(n)，得出信号x_r(n)。

其具体子步骤为：

(4.1)得出信号x(n)中第j个分量的幅度值的估计值

${\widehat{\mathrm{\σ}}}_j=|\underset{n=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\mathrm{\Σ}}}x\left(n\right)\·\exp [-j2\mathrm{\π}({\hat{a}}_{j,1}n+\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3)\left]\right|/N$

其中，N表示总采样点数；

(4.2)利用信号x(n)中第j个分量的幅度值的估计值、中心频率的估计值、调频率的估计值、调频率变化率的估计值，重构出信号x(n)中第j个分量s_j(n)，则有：

$s_j\left(n\right)={\widehat{\mathrm{\σ}}}_j\exp \left[j2\mathrm{\π}({\hat{a}}_{j,1}n+\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3)\right]$

在重构出信号x(n)中第j个分量s_j(n)之后，将矩阵signal的第1行第r+1列的元素替换为将矩阵signal的第2行第r+1列的元素替换为将矩阵signal的第3行第r+1列的元素替换为将矩阵signal的第4行第r+1列的元素替换为

(4.3)在信号x(n)中滤除第j个分量s_j(n)，得出信号x_r(n)。

具体地说，用信号x(n)乘以相位再对n快速傅里叶变换至频域f得到在处形成窄带脉冲的频域函数，用窄带的带阻滤波器把附近的窄带频谱滤除，把最强分量j的频谱的主瓣滤掉，将滤除最强分量的信号乘以构造新的信号x_r(n)。

$x_r\left(n\right)=x_{c,i}\left(n\right)\·\exp \left(j2\mathrm{\π}(\frac{1}{2}{\hat{a}}_{j,,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3)\right)$

$x_{c,j}\left(n\right)=\mathrm{IFFT}\left\{{\mathrm{Win}}_j\left({\hat{a}}_{j,1}\right)\·\mathrm{FFT}[x\left(n\right)\·\exp (-j2\mathrm{\π}(\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3))]\right\}$

其中，为加在信号x(n)中第j个分量上的频域矩形窗函数，FFT[·]表示进行傅里叶变换，IFFT{·}表示逆傅里叶变换，f_j，L为设定的频率值，f_j,R为设定的频率值，f_j,L和f_j,R的数值是根据左右的频谱宽度设定的。

步骤5，得出信号x_r(n)的总能量E，|x_r(n)|表示x_r(n)的模值；计算出信号x_r(n)的总能量E和初始多分量二次调频信号的总能量E₀的比值ε，如果ε大于或等于设定的门限值，则令x(n)＝x_r(n)，并令r的值自增1，返回至步骤2，重复执行步骤2至步骤5；否则，如果ε小于设定的门限值，则二次调频信号参数估计过程完成。

下面结合仿真实验对本发明的效果做进一步说明。

仿真实验1：初始多分量二次调频信号x(n)的信号形式为：

$x\left(n\right)=\exp \left(j2\mathrm{\π}(-\frac{1}{32}n+\frac{1}{10N^2}{n}^3)\right)+\exp \left(j2\mathrm{\π}(\frac{1}{16}n+\frac{7}{20N^2}{n}^3)\right),$

其中，n＝-256,-255,…,256；在仿真实验1中，共采样513个点(N为513)，两点之间的采样间隔为1s，采样频率为F_s＝1Hz。对初始多分量二次调频信号x(n)采用本发明进行参数估计，得出x(n)中的两个分量。参照图2，为仿真实验1中初始多分量二次调频信号x(n)的两个分量在调频率变化率-频率平面的分布示意图。图2中，横轴表示频率，单位为赫兹，纵轴表示调频率变化率，不同的灰度值表示不同的信号幅度值。参照图3，为仿真实验1中初始多分量二次调频信号x(n)的两个分量在调频率变化率-频率-信号幅度值三维空间的分布示意图。图3中，水平的两个轴分别表示频率(单位为赫兹)和调频率变化率，竖直的轴表示信号幅度值。由图2和图3可知，通过峰值点的位置即可估计原信号的调频率变化率和中心频率，可以看出，两分量能够很好地分离，很好地抑制了外部交叉项的干扰，能实现信号各分量的参数估计。

仿真实验2：了验证在噪声环境下，多分量QFM信号参数的估计能力，在仿真实验2中，设初始多分量二次调频信号为三个幅度相同的二次调频信号之和(即初始多分量二次调频信号有三个分量)。在初始多分量二次调频信号分别加入10dB，5dB和0dB的高斯白噪声，采样时间点n＝-256,-255,…,256，采样频率F_s＝1Hz。仿真实验2中初始多分量二次调频信号各分量对应的参数如表1所示，其中，σ代表初始多分量二次调频信号各分量的幅度值，a₁代表初始多分量二次调频信号各分量的三次项系数(调频率变化率)，a₂代表初始多分量二次调频信号各分量的二次项系数(调频率)，a₃代表初始多分量二次调频信号各分量的一次项系数(中心频率)。

表1 QFM信号真实值

在仿真实验2中，针对初始多分量二次调频信号中分别加入10dB、5dB和0dB的高斯白噪声得出的三组信号，采用本发明进行参数估计，得出初始多分量二次调频信号的各分量在调频率变化率-频率-信号幅度值三维空间的分布示意图。参照图4，为仿真实验2中初始多分量二次调频信号中加入0dB的高斯白噪声时得出的初始多分量二次调频信号的各分量在调频率变化率-频率-信号幅度值三维空间的分布示意图，图4中，水平的两个轴分别表示频率(单位为赫兹)和调频率变化率，竖直的轴表示信号幅度值。从图4中可以看出，得出的初始多分量二次调频信号的各分量在调频率变化率-频率-信号幅度值三维空间能够很好的区分。表2为仿真实验2中初始多分量二次调频信号中加入10dB的高斯白噪声时得出的各分量对应的参数估计表，表3为仿真实验2中初始多分量二次调频信号中加入5dB的高斯白噪声时得出的各分量对应的参数估计表，表4为仿真实验2中初始多分量二次调频信号中加入0dB的高斯白噪声时得出的各分量对应的参数估计表。表2至表4中，代表初始多分量二次调频信号各分量的幅度值估计值，代表初始多分量二次调频信号各分量的三次项系数(调频率变化率)估计值，代表初始多分量二次调频信号各分量的二次项系数(调频率)估计值，代表初始多分量二次调频信号各分量的一次项系数(中心频率)估计值。对表1至表4进行对比可以看出，在加入不同信噪比的高斯白噪声后，采用本发明得出的参数的估计结果相差很小，说明了本发明可以在较低的信噪比情况下检测QFM信号并估计其参数，可见本文所提出的方法具有很高的估计精度和较低的信噪比门限。

表2 10dB信噪比条件下QFM信号估计值

表3 5dB信噪比条件下QFM信号估计值

表4 0dB信噪比条件下QFM信号估计值

当信号分量个数较少且信噪比较高时，本发明与PHMT方法和PGCPF方法的估计精度相差不大，但是当分量个数较多或者多分量信号的二次项系数差别较小时，本发明具有更为准确的估计结果，为此对两种情况下参数估计性能分别进行验证，参见仿真实验3和仿真实验4。

仿真实验3：为了验证本发明在多分量QFM信号情况下有较好的参数估计能力，在仿真实验3中，设初始多分量二次调频信号为三个幅度相同的二次调频信号之和(即初始多分量二次调频信号有三个分量)。采样时间点n＝-256,-255,…,256采样频率F_s＝1Hz。仿真实验3中初始多分量二次调频信号各分量对应的参数如表1所示，在仿真实验3中，采用PHMT方法得出初始多分量二次调频信号各分量的峰值位置。参照图5，为仿真实验3中采用PHMT方法得出的初始多分量二次调频信号各分量的峰值位置的估计结果示意图。图5中，横轴表示初始多分量二次调频信号各分量的峰值位置，纵轴表示归一化相对幅度。从图5可以看出此时三个峰值的位置分别出现在128，149和369，但三个峰值正确的位置应分别出现在122，347和437。而由仿真实验2的结果(图4)可以看出本发明相比PHMT方法具有较高的估计精度。

仿真实验4：为了验证本发明在多分量QFM信号的二次项系数差别较小情况下有较好的参数估计能力。设初始多分量二次调频信号为两个幅度相同的二次调频信号之和(即初始多分量二次调频信号有两个分量)，n＝-256,-255,…,256采样频率F_s＝1Hz。仿真实验4中初始多分量二次调频信号各分量对应的参数如表5所示，表5中，σ代表初始多分量二次调频信号各分量的幅度值，a₁代表初始多分量二次调频信号各分量的三次项系数(调频率变化率)，a₂代表初始多分量二次调频信号各分量的二次项系数(调频率)，a₃代表初始多分量二次调频信号各分量的一次项系数(中心频率)。

表5 QFM信号真实值

在仿真实验4中，分别采用本发明和PGCPF方法得出初始多分量二次调频信号各分量的峰值位置。参照图6，为仿真实验4中采用PHMT方法得出的初始多分量二次调频信号各分量的峰值位置的估计结果示意图。图6中，横轴表示初始多分量二次调频信号各分量的峰值位置，纵轴表示归一化相对幅度。从图6中可以看出此时只有一个较高的峰值(此时的估计值为具有相同二次项系数的信号的三次项系数的平均值)，但正确的结果应该有两个较高的峰值(峰值位置对应于三次项系数)。参照图7，为仿真实验4中采用PGCPF方法得出的初始多分量二次调频信号各分量在调频率变化率-频率-信号幅度值三维空间的分布示意图。图7中，水平的两个轴分别表示频率(单位为赫兹)和调频率变化率，竖直的轴表示信号幅度值。从图7可以看出，此时有两个较高的峰值，且峰值位置反映参数大小。因此，当多分量QFM信号的二次项系数差别较小时，说明本发明相对于PGCPF方法具有更好的参照估计效果。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.基于四阶非对称乘积型核函数的QFM信号参数估计方法，其特征在于，包括以下步骤：

步骤1，获取初始多分量二次调频信号x(n)，n表示采样序列时刻点，N表示总采样点数；得出初始多分量二次调频信号的总能量E₀；设置循环次数变量r，r＝0，1，2，...，当r＝0时，执行步骤2；

步骤2，构造信号x(n)的时延函数x(n+m)，m表示信号x(n)时延的时刻点数，得出函数y(n，m)，y(n，m)＝x(n+m)x(-n-m)；将信号x(n)中具有最大幅度值的分量记为信号x(n)的第j个分量，根据函数y(n，m)，得出信号x(n)中具有最大幅度值的分量的调频率估计值

根据以下公式得出信号x(n)中具有最大幅度值的分量的调频率估计值

${\hat{a}}_{j,2}=\underset{f_2}{\mathrm{argmax}}\left|PY\left(f_2\right)\right|$

$PY\left(f_2\right)=\underset{n_1=-\frac{N}{4}}{\overset{\frac{N}{4}}{\Π}}Y(n_1,f_2)$

$Y(n,f_2)=\underset{m=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\Σ}}y(n,m)e^{-j2{\mathrm{\πf}}_2(m^2+2nm)}$

其中，f₂表示时延m对应的频域频率变量，|PY(f₂)|表示PY(f₂)的模值，

步骤3，利用信号x(n)中具有最大幅度值的分量的调频率估计值对信号x(n)进行补偿，得到补偿后信号x_d1(n)；利用补偿后信号x_d1(n)，构造四阶非对称乘积型核函数y₁(n，m)：

$y_1(n,m)=x_{d1}(n+\frac{2}{3}m)x_{d1}(n+m){x_{d1}}^{*}(n+\frac{2+\sqrt{22}}{6}m){x_{d1}}^{*}(n+\frac{2-\sqrt{22}}{6}m)$

其中，上标*表示取共轭；

对四阶非对称乘积型核函数y₁(n，m)进行相位匹配变换，得到峰值函数Y₁(f₃，m)，f₃表示对应于时刻n的频域频率变量；对峰值函数Y₁(f₃，m)进行由m时域至f₁频域的快速傅里叶变换，得到峰值函数Y₂(f₃，f₁)，f₁表示对应于时延m的频域频率变化率；将峰值函数Y₂(f₃，f₁)的频谱图的最大幅度值在f₁轴上的投影作为信号x(n)中具有最大幅度值的分量的中心频率的估计值将峰值函数Y₂(f₃，f₁)的频谱图的最大幅度值在f₃轴上的投影作为信号x(n)中具有最大幅度值的分量的调频率变化率的估计值

步骤4，得出信号x(n)中第j个分量的幅度值的估计值利用信号x(n)中第j个分量的幅度值的估计值、中心频率的估计值、调频率的估计值、调频率变化率的估计值，重构出信号x(n)中第j个分量s_j(n)；在信号x(n)中滤除第j个分量s_j(n)，得出信号x_r(n)；

步骤5，得出信号x_r(n)的总能量E，计算出信号x_r(n)的总能量E和初始多分量二次调频信号的总能量E₀的比值ε，如果ε大于或等于设定的门限值，则令x(n)＝x_r(n)，并令r的值自增1，返回至步骤2；否则，说明二次调频信号参数估计过程完成。

2.如权利要求1所述的基于四阶非对称乘积型核函数的QFM信号参数估计方法，其特征在于，在步骤1中，初始多分量二次调频信号的总能量E₀为：

$E_0=\underset{n=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\Σ}}|x\left(n\right){|}^2$

其中，|x(n)|表示初始多分量二次调频信号的模值。

3.如权利要求1所述的基于四阶非对称乘积型核函数的QFM信号参数估计方法，其特征在于，在步骤3中，补偿后信号x_d1(n)为：

$x_{d1}\left(n\right)=x\left(n\right)\·\exp (-j2\mathrm{\π}(\frac{1}{2}{\hat{a}}_{j,2}{n}^2\left)\right).$

4.如权利要求1所述的基于四阶非对称乘积型核函数的QFM信号参数估计方法，其特征在于，在步骤3中，峰值函数Y₁(f₃，m)为：

$Y_1(f_3,m)=\underset{n=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\Σ}}{y}_1(n,m)e^{-j2{\mathrm{\πf}}_3{n}^{2}m}\left|n\right|$

其中，f₃表示对应于时刻n的频域频率变量，|n|表示n的绝对值。

5.如权利要求1所述的基于四阶非对称乘积型核函数的QFM信号参数估计方法，其特征在于，所述步骤4的具体子步骤为：

(4.1)得出信号x(n)中第j个分量的幅度值的估计值

${\widehat{\mathrm{\σ}}}_j=|\underset{n=-\frac{N-1}{2}}{\overset{\frac{N-1}{2}}{\Σ}}x\left(n\right)\·\exp \[-j2\mathrm{\π}({\hat{a}}_{j,1}n+\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3)\]|/N$

其中，N表示总采样点数；

(4.2)利用信号x(n)中第j个分量的幅度值的估计值、中心频率的估计值、调频率的估计值、调频率变化率的估计值，重构出信号x(n)中第j个分量s_j(n)，则有：

$s_j\left(n\right)={\widehat{\mathrm{\σ}}}_j\exp \[j2\mathrm{\π}({\hat{a}}_{j,1}n+\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3)\];$

(4.3)在信号x(n)中滤除第j个分量s_j(n)，得出信号x_r(n)；

$x_r\left(n\right)=x_{c,i}\left(n\right)\·\exp \left(j2\mathrm{\π}\right(\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3\left)\right)$

$x_{c,j}\left(n\right)=IFFT\{{\mathrm{Win}}_j\left({\hat{a}}_{j,1}\right)\·FFT\[x\left(n\right)\·\exp (-j2\mathrm{\π}(\frac{1}{2}{\hat{a}}_{j,2}{n}^2+\frac{1}{6}{\hat{a}}_{j,3}{n}^3))\]\}$

其中，FFT[·]表示进行傅里叶变换，IFFT{·}表示逆傅里叶变换，f_j，L为设定的频率值，f_j，R为设定的频率值。