CN105785324B - 基于mgcstft的线性调频信号参数估计方法 - Google Patents

Abstract

本发明公开了一种基于MGCSTFT的线性调频信号参数估计方法，利用广义柯西分布，构造了一类适用于α稳定分布噪声环境的损失函数，通过最大似然估计理论准确实现了线性调频信号的参数估计。具体步骤包括：1、采集线性调频信号，2、确定判别门限，3、判断是否含有脉冲噪声，4、设置损失函数，5、做最大似然广义柯西短时傅里叶变换MGCSTFT，6、提取线性调频信号的参数。本发明克服了传统的时频分析技术在α稳定分布噪声环境下估计信号参数失效的问题，提高了参数估计精度。

Description

基于MGCSTFT的线性调频信号参数估计方法

技术领域

本发明属于通信信号处理领域，更进一步涉及到雷达信号处理技术领域中的一种基于最大似然广义柯西短时傅里叶变换MGCSTFT(maxmium-likehood generalized CauchyShort Time Fourier Transform MGCSTFT)的线性调频信号参数估计方法。本发明利用时频分析技术，在时频域经过迭代去除脉冲噪声，可以准确实现在α稳定分布噪声环境中线性调频信号的参数提取。

背景技术

线性调频信号LFM(Linear Frequency Modulation LFM)，又称chirp信号，作为非平稳信号的典型代表，因为保留连续信号和脉冲的特性，广泛应用于雷达，语音，声呐和生物医学等方面。在数字接收的情况下，对线性调频信号初始频率和调频斜率参数进行精确估计，可以实现电子侦察***中的目标检测和识别。目前，研究者针对非平稳信号的参数估计先后提出了多种有效的处理方法，其中包括以短时傅里叶变换(STFT)为代表的线性时频分析方法和Wigner-Ville分布(WVD)为代表的双线性时频分析方法。短时傅里叶变换作为最早提出的一类时频分析方法，因为可以结合各种窗函数的优点，且对于多分量信号的处理没有交叉项干扰，因而被广泛应用于非平稳信号的处理中。但是STFT的时频聚集性不够良好，容易受到噪声的干扰，尤其在α稳定分布噪声影响下，STFT方法对非平稳信号参数的估计性能大幅下降。

南京信息工程大学申请的专利“基于分数低阶统计量的正交小波盲均衡方法”(申请号CN201110208437A，公开号CN102355435A)中公开了一种分数低阶统计量的方法。该方法利用分数低阶统计量来抑制α稳定分布噪声，利用信源的先验信息，在迭代过程中自适应修正模值，并且对均衡器输入信号进行了正交小波变换，在减小输入信号自相关性的同时，提高了***的均衡性。虽然该方法利用的降阶技术在一定程度上抑制了较大脉冲，但是该方法存在的不足之处是，当脉冲噪声的强度增强，该方法性能退化，且分数低阶算子的取值没有严格的选取标准。

朱敏等人提出的L-DFT估计方法(Alpha稳定分布噪声下基于L-DFT的LFM-BPSK复合调制信号参数估计，振动与冲击，2015.9)，利用平方倍频法消除编码调相，L-DFT方法抑制脉冲噪声的同时，估计线性调频信号的起始频率和调制斜率，且效果明显优于分数低阶统计量，但此类方法存在的不足是，L估计求解过程复杂，对强脉冲噪声较为敏感，当广义信噪比降低时，此方法性能下降，难以满足非平稳信号参数估计的要求。

综上所述，对于线性调频信号在脉冲噪声下的参数估计方法，目前已有的时频分析技术在脉冲噪声的环境下不具有稳健性，容易受到脉冲噪声的干扰，从而使信号的参数估计达不到实际要求的精度，尤其当广义信噪比降低时，基于非线性滤波的方法性能退化，甚至失效。

发明内容

本发明的目的在于克服上述已有时频分析技术在脉冲噪声环境下无法实现信号参数估计的不足，在传统的短时傅里叶变换时频分析技术的基础上，利用广义柯西分布，构造了一类适用于α稳定分布噪声环境的损失函数，通过最大似然估计理论，得到了一种改进的短时傅里叶变换，即最大似然广义柯西短时傅里叶变换MGCSTFT。此类方法估计线性调频信号的参数，可以提高α稳定分布噪声下线性调频信号参数估计的准确性和稳定性，从而提高线性调频信号的参数估计性能。

实现本发明的具体思路是：将广义柯西分布用来拟合α稳定分布噪声，由其概率密度函数得到最大似然估计理论中的损失函数，此类损失函数可以较好适应α稳定分布噪声环境。在短时傅里叶变换的时频平面，线性调频信号的最大似然估计由一个迭代的过程来实现。通过迭代，得到线性调频信号的最大似然广义柯西短时傅里叶变换MGCSTFT时频分布图，对时频分布图进行霍夫变换直线检测，最终得到线性调频信号的初始频率和调频斜率。

本发明的具体步骤包括如下：

(1)采集线性调频信号：

信号采集***通过脉压雷达的接收机设备，采集雷达天线中任意一段含有实际噪声的线性调频信号；

(2)确定判别门限：

采用局部幅值特征方法，得到局部阈值，将该阈值作为判别门限；

(3)判断幅值统计模块采集的线性调频信号局部幅值是否大于判别门限，若是，则认为线性调频信号中含有α稳定分布噪声，执行步骤(4)，否则，对采集的线性调频信号做短时傅里叶变换后，执行步骤(6b)；

(4)按照下式，设置最大似然估计的损失函数：

F(e)＝-logp_v(t)

其中，F(e)表示最大似然估计的损失函数，e表示估计误差，t表示采集信号的采样时间，v(t)表示采集信号中的背景噪声，p_v(t)表示背景噪声v(t)的概率密度函数；

(5)做最大似然广义柯西短时傅里叶变换MGCSTFT：

(5a)对步骤(3)中得到的脉冲指示信号做短时傅里叶变换，将短时傅里叶变换的结果，作为最大似然广义柯西短时傅里叶变换MGCSTFT迭代的初始值；

(5b)按照下式，依次迭代得到最大似然广义柯西短时傅里叶变换MGCSTFT值：

其中，M⁽ⁱ⁾表示第i次迭代得到的最大似然广义柯西短时傅里叶变换MGCSTFT值，i表示迭代次数，Q表示窗长函数在离散时间点的累加值，N表示采集信号的总采样点数，N的大小由采集信号的长度和采样间隔确定，m表示离散时间初始点，ω表示窗函数，n表示离散时间采样点，Δt表示离散时间段，t表示采集信号的采样时间，x表示采集信号，f表示采集信号的采样频率，e表示误差函数，k表示广义柯西分布的尺度参数，k的大小由α稳定分布噪声的系数α确定，p表示广义柯西分布的阶数，取值范围为[1,2]；

(5c)判断相邻两次迭代值的相对误差是否小于门限值η，若是，则执行步骤(6)，否则，执行步骤(5b)；

(6)提取线性调频信号的参数：

(6a)在时间-频率平面画出与最大似然广义柯西短时傅里叶变换MGCSTFT值对应的时频分布图；

(6b)采用霍夫变换直线检测方法，在时频分布图上，提取线性调频信号的参数。

与现有技术相比，本发明具有如下优点：

第一，由于本发明采用的做最大似然广义柯西短时傅里叶变换MGCSTFT，在时频平面通过迭代抑制α稳定分布噪声，克服了现有技术中采用的时频分析技术在α稳定分布噪声中性能退化的缺点，使得本发明可以在较低的广义信噪比下，有效地提高线性调频信号参数估计的准确率，获得更小的均方根误差，从而提高线性调频信号参数估计的性能。

第二，由于本发明采用的判断相邻两次迭代值的相对误差是否小于门限值η，η取值越小，信号参数的估计精度越高，克服了现有技术中无法灵活地调整信号参数估计精度的局限性，使得本发明对线性调频信号的参数估计达到更小的精度，从而有效地凸显线性调频信号的局部信息。

附图说明

图1是本发明的流程图；

图2是本发明的仿真图；

图3是本发明的仿真效果对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1，采集线性调频信号。

信号采集***通过脉压雷达的接收机设备，采集雷达天线中任意一段含有实际噪声的线性调频信号。

信号采集***通过接收天线，选取任意一段含有原始线性调频信号s(t)和脉冲噪声v(t)的接收信号，将所选取的接收信号作为采集信号，即接收信号：x(t)＝s(t)+v(t)，其中线性调频信号s(t)的模型可表示为：

s(t)＝Aexp(j2πft+(jπkt²),t∈[0,T]

其中，A，f，k分别表示信号的幅值、初始频率和调频斜率，T表示时间区间。LFM信号频率随时间呈线性变化。

步骤2，确定判别门限：

采用局部幅值特征方法，得到局部阈值，将该阈值作为判别门限。

局部阈值的确定步骤如下：

第一步，将采集信号输入判别模块，作为判别信号。

第二步，设置一个固定长度为W的检测窗口，长度W的取值范围为大于0小于线性调频信号采样点总数的值。

第三步，利用检测窗口的时域平滑，将线性调频信号按时间段截断，划分为多个等长度时间段、互不重叠的子信号，采用局部幅值特征方法计算子信号的幅度均值，得到局部阈值。

步骤3，判断幅值统计模块采集的线性调频信号局部幅值是否大于判别门限，若是，则认为线性调频信号中含有α稳定分布噪声，执行步骤4，否则，直接对采集的信号做短时傅里叶变换后，执行步骤6的第二步；

步骤4，按照下式，设置最大似然估计的损失函数：

F(e)＝-logp_v(t)

其中，F(e)表示最大似然估计的损失函数，e表示估计误差，t表示采集信号的采样时间，v(t)表示采集信号中的背景噪声，p_v(t)表示背景噪声v(t)的概率密度函数。

在估计含有噪声的信号参数时，基于最大似然估计理论的方法是最优的，其损失函数F(e)＝-logp_v(t)，其中，p_v(t)表示噪声v(t)的概率密度函数。在高斯背影噪声中，由其概率密度函数正态分布可得损失函数F(e)＝|e|²，最大似然理论表明，可以得到信号的最小均方误差估计，信号s(n)由均值滤波的输出得到，即：

其中，s(n)表示待估计的信号，表示信号s(n)的最大似然估计，M表示离散采样点的移动幅度，n表示离散采样点，x(k)表示含有噪声的信号，mean表示均值滤波函数。

但是，最大似然估计对偏离假定噪声模型的干扰非常敏感，这就导致了在α稳定分布噪声下，基于最大似然估计理论的常规分析方法性能退化严重，因此本发明提出了一类新的损失函数，此类损失函数基于广义柯西分布，其概率密度函数为：

f(x)＝a(k^p+|x|^p)^-2/p

其中，k表示尺度参数，在区间(0,2]上取值；Γ(·)表示伽马函数，广义柯西分布包含两个特例，即p＝1时的Merid分布和p＝2时的柯西分布，由最大似然估计理论，可得损失函数函数为：F(e)＝log(k^p+|e|^p)，其中k＞0，0＜p≤2，e表示误差。

步骤5，做最大似然广义柯西短时傅里叶变换MGCSTFT：

第一步，令i＝0，对采集信号x(t)做短时傅里叶变换STFT，作为迭代的初始值，即M⁽⁰⁾＝STFT_x，其中，M⁽⁰⁾表示最大似然广义柯西短时傅里叶变换MGCSTFT的初始值，STFT_x表示采集信号x(t)的短时傅里叶变换。

第二步，令i＝i+1，迭代的步骤如下：

其中，M⁽ⁱ⁾表示i次迭代时得到的最大似然广义柯西短时傅里叶变换MGCSTFT值，i表示迭代次数，Q表示窗长函数在离散时间点的累加值，Q的取值由确定，N表示采集信号的总采样点数，N的大小由采集信号的长度和采样间隔确定，m表示离散时间初始点，ω表示窗函数，n表示离散时间采样点，Δt表示离散时间段，t表示采集信号的采样时间，x表示采集信号，f表示采集信号的采样频率，e表示误差函数，e^(i-1)＝x(nΔt)e^-j2πfnΔt-M^(i-1)，k表示广义柯西分布的尺度参数，k的大小由α稳定分布噪声的系数α确定，即p表示广义柯西分布的阶数，取值范围为[1,2]；

第三步，判断相邻两次迭代值的相对误差是否小于门限值η，若其中，M⁽ⁱ⁾表示i次迭代时得到的最大似然广义柯西短时傅里叶变换MGCSTFT值，则执行步骤6，否则，执行步骤5的第二步；

门限值η，由信号参数的估计精度决定，取为0.1；

步骤6，提取线性调频信号的参数：

第一步：对最大似然广义柯西短时傅里叶变换MGCSTFT得到的时频图进行霍夫变换，检索到峰值点的极坐标(ρ,θ)，其中，ρ表示极坐标的极径，θ表示极坐标的极角；

第二步：由坐标变换公式得到线性调频信号的参数，其中，f₀表示线性调频信号的调频斜率，k₀表示线性调频信号的初始频率，ρ表示极坐标的极径，θ表示极坐标的极角。

下面结合仿真实验对本发明做进一步的描述。

1、仿真条件：

本发明的仿真信号为线性调频信号，其参数为：采样频率f_s＝10kHz，初始频率f₀＝1kHz，调频斜率k₀＝10kHz/s，采样点数N＝1600。仿真软件环境为Intel(R)Pentium(TM)G3260 [email protected]，Windows 7 32bit操作***下的Matlab7.8。

本发明采用归一化均方根误差NRMSE来评价参数估计方法的性能，假设MonteCarlo仿真实验的次数为R，所估计的参数为θ，且在第i次实验中的估计值为，则参数θ估计的归一化均方根误差为

本发明采用广义信噪比来度量信号和脉冲噪声能量的大小，由于α稳定分布不存在有限的二阶矩，在常规信噪比下噪声方差失去意义，故重新定义广义信噪比：式中，表示信号能量；γ_v为α稳定分布噪声的分散系数。

2、仿真内容与结果分析：

仿真1：采用MGCSTFT方法估计线性调频信号的参数。

在特征指数α＝0.8的α稳定分布噪声下，线性调频信号采用MGCSTFT方法的计算机仿真如图，取p＝1，η＝0.1，图2(a)为不含噪声的线性调频原始信号的短时傅里叶变换图,由图2(a)可知，线性调频信号在时频域清晰明显，频率随时间线性变化，呈现直线的特征。

图2(b)是线性调频信号中加入高斯白噪声后，对混合信号做短时傅里叶变换的时频分布图，信噪比SNR＝3dB。STFT方法对于高斯白噪声具有一定的稳健性，在时频分布图上线性调频信号轮廓清晰，可见常规的时频分析技术可以处理含有高斯白噪声的信号。

图2(c)是加入广义信噪比GSNR＝3dB，脉冲强度α＝0.8的α稳定分布噪声后，线性调频信号的时频图，可见线性调频信号被噪声完全淹没，时频平面出现混乱。线性调频信号受脉冲噪声影响严重，与高斯白噪声影响下的时频图完全不同，因此传统的时频分析方法STFT难以对含有脉冲噪声的线性调频信号进行参数估计，验证了基于高斯模型的传统处理方法在脉冲噪声环境下失效。

图2(d)是α稳定分布噪声的线性调频信号应用MGCSTFT方法，在经过4次迭代后，得到的时频分布图。从图2(d)可以看到，经过MGCSTFT方法处理后，脉冲噪声已被去除，线性调频信号在时频平面清晰明确。

图2(e)是在图2(c)时频分布图的基础上，利用霍夫直线检测方法，提取线性调频信号的参数，从图2(e)可以看出，霍夫变换已经搜索到峰值点，其顶点对应的坐标(ρ,θ)经过坐标变换公式后，得到了线性调频信号的初始频率和调频斜率。

仿真2：各种方法的误差比较。

改变广义信噪比GSNR，在广义信噪比从-2dB到8dB变化时，将MGCSTFT方法得出的线性调频信号初始频率和调频斜率的归一化均方根误差。将MGCSTFT方法分别与基于Myriad滤波器的STFT(MYRSTFT)，Merid滤波器的STFT(MERSTFT)，传统的STFT处理方法得到的归一化均方根误差做对比，经过100次Monte Carlo仿真实验，其仿真结果如图3(a)，图3(b)。

从图3(a)，图3(b)可以看到，当GSNR≥1dB时，MGCSTFT方法可以准确估计线性调频信号的调频斜率和初始频率；当GSNR≥6dB时，基于Merid滤波器的STFT可准确获得线性调频信号的调频斜率，GSNR≥5dB时，可准确估计线性调频信号的初始频率；当GSNR≥8dB时，基于Myriad滤波器的STFT才可以得到线性调频信号调频斜率和初始频率的准确估计；传统的STFT方法受强脉冲噪声影响严重，性能严重退化，难以对信号的参数做出准确估计。基于Myriad滤波器和Merid滤波器等多种方法虽然在信噪比较高的情况下均可较为准确进行线性调频信号参数估计，但在强脉冲噪声和低信噪比下，估计性能都有所下降。

综上所述，传统的时频分析方法STFT在强脉冲噪声的影响下无法估计线性调频信号的参数，通过Myriad滤波器和Merid滤波器抑制脉冲噪声后，虽然MYRSTFT和MERSTFT都可以在一定程度上抑制脉冲噪声，但是在较低的广义信噪比下，两者的估计仍然出现较大的误差。MGCSTFT方法不需要在时域进行滤波，避免了应用滤波器造成的线性调频信号参数估计性能下降的缺点，在时频平面通过迭代后，准确地实现了信号参数的估计。因此，本发明提出的MGCSTFT新方法适用于α稳定分布噪声下线性调频信号的参数估计，且估计性能和精度良好，优于其他常规的处理方法。

Claims (3)

1.基于MGCSTFT的线性调频信号参数估计方法，包括如下步骤：

(1)采集线性调频信号：

信号采集***通过脉压雷达的接收机设备，采集雷达天线中任意一段含有实际噪声的线性调频信号；

(2)确定判别门限：

采用局部幅值特征方法，得到局部阈值，将该阈值作为判别门限；

(3)判断幅值统计模块采集的线性调频信号局部幅值是否大于判别门限，若是，则认为线性调频信号中含有α稳定分布噪声，执行步骤(4)，否则，执行步骤(6b)；

(4)按照下式，设置最大似然估计的损失函数：

F(e)＝-logp_v(t)

其中，F(e)表示最大似然估计的损失函数，e表示估计误差，t表示采集信号的采样时间，v(t)表示采集信号中的背景噪声，p_v(t)表示背景噪声v(t)的概率密度函数；

(5)做最大似然广义柯西短时傅里叶变换MGCSTFT：

(5a)对步骤(3)中得到的脉冲指示信号做短时傅里叶变换，将短时傅里叶变换的结果，作为最大似然广义柯西短时傅里叶变换MGCSTFT迭代的初始值；

(5b)按照下式，依次迭代得到最大似然广义柯西短时傅里叶变换MGCSTFT值：

<mrow> <msup> <mi>M</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msup> <mo>=</mo> <mfrac> <mn>1</mn> <mi>Q</mi> </mfrac> <mo>&amp;CenterDot;</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mfrac> <mrow> <mi>&amp;omega;</mi> <mrow> <mo>(</mo> <mi>n</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>-</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>n</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mi>n</mi> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> </mrow> <mrow> <mo>|</mo> <msup> <mi>e</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msup> <msup> <mo>|</mo> <mi>p</mi> </msup> <mo>+</mo> <msup> <mi>k</mi> <mi>p</mi> </msup> </mrow> </mfrac> </mrow>

其中，M⁽ⁱ⁾表示第i次迭代得到的最大似然广义柯西短时傅里叶变换MGCSTFT值，i表示迭代次数，Q表示窗长函数在离散时间点的累加值，N表示采集信号的总采样点数，N的大小由采集信号的长度和采样间隔确定，m表示离散时间初始点，ω表示窗函数，n表示离散时间采样点，Δt表示离散时间段，t表示采集信号的采样时间，x表示采集信号，f表示采集信号的采样频率，e表示误差函数，k表示广义柯西分布的尺度参数，k的大小由α稳定分布噪声的系数α确定，p表示广义柯西分布的阶数，取值范围为[1,2]；

(5c)判断相邻两次迭代值的相对误差是否小于门限值η，若是，则执行步骤(6)，否则，执行步骤(5b)；

(6)提取线性调频信号的参数：

(6a)在时间-频率平面画出与最大似然广义柯西短时傅里叶变换MGCSTFT值对应的时频分布图；

(6b)采用霍夫变换直线检测方法，在时频分布图上，提取线性调频信号的参数。

2.根据权利要求1所述的基于MGCSTFT的线性调频信号参数估计方法，其特征在于：步骤(5c)所述的门限值η是由信号参数的估计精度决定，取为0.1。

3.根据权利要求1所述的基于MGCSTFT的线性调频信号参数估计方法，其特征在于：步骤(6b)所述的霍夫变换直线检测方法的具体步骤如下：

第一步：对最大似然广义柯西短时傅里叶变换MGCSTFT得到的时频图进行霍夫变换，检索到峰值点的极坐标(ρ,θ)，其中，ρ表示极坐标的极径，θ表示极坐标的极角；

第二步：由坐标变换公式得到线性调频信号的参数，其中，f₀表示线性调频信号的调频斜率，k₀表示线性调频信号的初始频率，ρ表示极坐标的极径，θ表示极坐标的极角。