CN103675758A

CN103675758A - 一种双曲调频信号周期斜率和起始频率估计方法

Info

Publication number: CN103675758A
Application number: CN201310652510.2A
Authority: CN
Inventors: 方世良; 姚帅; 王晓燕; 韩宁; 王莉
Original assignee: Southeast University
Current assignee: Southeast University
Priority date: 2013-12-05
Filing date: 2013-12-05
Publication date: 2014-03-26
Anticipated expiration: 2033-12-05
Also published as: CN103675758B

Abstract

本发明公开了一种双曲调频信号周期斜率和起始频率估计方法，该方法包括以下步骤：第一步：获取数据序列；第二步：参数初始化；第三步：计算第i个短时窗内数据的功率谱；第四步：采用最大线谱法估计第i个短时窗内数据的瞬时频率f_i；第五步：估计第i个短时窗内数据的过零点间隔g_i=1/f_i；第六步：判断是否处理完所有短时窗的数据：如未处理完,返回第三步，否则转入第七步；第七步：对过零点间隔估计序列{g_i,i＝1,2,…,I}进行滑动中值滤波得第八步：计算第k次迭代权重第九步：判断是否满足迭代停止条件：如不满足，返回第八步，否则转入第十步；第十步：计算出周期斜率和起始频率。该方法无需复杂的计算和参数搜索，稳健性强，可以实现参数的快速、高精度估计。

Description

一种双曲调频信号周期斜率和起始频率估计方法

技术领域

本发明属于信号处理领域，具体涉及一种双曲调频信号周期斜率和起始频率估计方法。

背景技术

双曲调频信号具有多普勒不变性，这一性质使其特别适用于检测高速小目标，这使得双曲调频信号在水声和雷达等领域得到了广泛的应用。周期斜率和起始频率是表征线双曲调频信号频率特性的两个基本参数，如果能够估计得到这两个参数，在已知信号脉宽的条件下，即可恢复出获取的双曲调频信号，这对水声对抗和雷达对抗有着十分重要的意义，因此其估计问题是水声和雷达信号处理领域的重要研究内容。

目前用于估计双曲调频信号参数的方法主要有最大似然法，非线性最小二乘法和时频分析与图像处理相结合的方法。最大似然法和非线性最小二乘法都需要解非线性方程，而对于双曲调频信号来说，这两种方法需要求解的非线性方程都不存在解析解，因此需要使用数值解法，这就加大了求解的难度和计算的复杂度。时频分析与图像处理相结合的方法为了达到较好的估计结果，需要重复使用时频滤波器，这使得该算法的计算量大大增加，限制了其工程实用性。

在进行双曲调频信号周期斜率和起始频率估计前，需要先进行数据采集和信号检测，由数据采集部分可以得到采样频率f_s，通过信号检测可以得到信号的起始和终止时刻，检测到的信号的终止和起始时刻的差值即为信号脉宽，信号脉宽所对应的采样点数即为N，在检测到信号起始时刻并得信号脉宽所对应的采样点数N后，即可进行双曲调频信号周期斜率和起始频率估计。

发明内容

技术问题：本发明提供了一种无需进行复杂的计算和参数搜索，可在保证参数快速估计前提下，提高参数估计精度的双曲调频信号周期斜率和起始频率估计方法。

技术方案：本发明的双曲调频信号周期斜率和起始频率估计方法，包括以下步骤：

第一步，获取待处理的数据序列x(n),n＝0,1,…,N-1：从传感器接收N个采样点的实时采集数据作为待处理的数据序列x(n),n＝0,1,…,N-1，或从存储器中提取从检测到信号时刻起始的N个采样点的数据作为待处理的数据序列x(n),n＝0,1,…,N-1，所述的N为检测到的双曲调频信号脉宽长度所对应的采样点个数；

第二步，参数初始化：设置短时窗长M、短时窗移动步进L、最大迭代次数门限K和精度控制指标ε，计算出总的短时窗个数

表示向下取整运算，初始化短时窗序号i＝1，所述短时窗长M取值为2的整数次幂且满足M＜N/4，L取值为

K取值为大于等于2的正整数，ε取值为小于等于0.1的正数；

第三步，对第i个短时窗内的数据序列x_i(m)做离散傅里叶变换并计算其功率谱Y_i(l₂)：第i个短时窗内的数据序列为x_i(m)=x(n_i),m＝0,1,…M-1，n_i＝(i-1)L,(i-1)L+1,…,(i-1)L+M-1，用下列式1对x_i(m)做离散傅里叶变换：

X_{i} (l_{1}) = Σ_{m = 0}^{M - 1} x_{i} (m) e^{- j \frac{2 π}{M} ml}, l_{1} = 0,1 \cdot \cdot \cdot, M

式1

其中X_i(l₁)表示离散傅里叶变换的结果，j表示虚数单位，即

l₁为X_i(l₁)的离散频率序号，则第i个短时窗内的数据序列x_i(m)的功率谱Y_i(l₂)为：

Y_{i} (l_{2}) = \frac{1}{M} {| X_{i} (l_{1}) |}^{2},

l₁＝l₂且l₂＝0,1,2…M/2-1 式2

其中l₂为Y_i(l₂)的离散频率序号；

第四步，采用最大线谱法，按照下式估计出第i个短时窗内数据序列的瞬时频率估计值f_i：

f_i＝(L_i-1)Δf 式3

其中L_i为所有功率谱Y_i(l₂)，l₂＝0,1,2…M/2-1中的最大值对应的离散频率序号，Δf为短时窗长度为M的离散傅里叶变换的频率分辨率，Δf＝f_s/M，f_s为采样频率；

第五步，对所述第四步得到的瞬时频率估计值f_i取倒数，作为第i个短时窗内数据序列的过零点间隔估计值g_i：

g_i＝1/f_i 式4；

第六步，判断是否处理完所有短时窗的数据序列：如果i≤I-1，则令i＝i+1，并返回到第三步，否则令迭代次数k＝0，并进入第七步；

第七步，对过零点间隔估计值的序列{g_i,i＝1,2,…,I}进行滑动中值滤波得到

{g_{i}^{m}, i = 1,2 \cdot \cdot \cdot, I},

的表达式为：

g_{i}^{m} = \{\begin{matrix} g_{1} & i = 1 \\ median (g_{i - 1}, g_{i}, g_{i + 1}) & 2 < i < I \\ g_{I} & i = I \end{matrix}

式5

其中median(g_i-1,g_i,g_i+1)是取g_i-1,g_i和g_i+1的中值；

第八步，计算第k次迭代中第i个短时窗内数据序列的瞬时频率估计值所对应的加权权重

的计算过程如下：

d_{i}^{k} = | g_{i}^{m} - a^{k - 1} - b^{k - 1} i |, k &GreaterEqual; 1

式6

λ_{i}^{k} = (d_{i}^{k} - d_{\min}^{k}) / d_{mean}^{k}, k &GreaterEqual; 1

式7

w_{i}^{k} = \{\begin{matrix} 1 & k = 0 \\ \frac{1}{λ_{i}^{k} + δ_{1}} & k &GreaterEqual; 1 \end{matrix}

式8

式6中，a^k-1和b^k-1的表达式为：

a^{k - 1} = \frac{Σ_{i = 1}^{I} i^{2} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} g_{p}^{m} - Σ_{i = 1}^{I} {iw}_{i}^{k - 1} Σ_{p = 1}^{I} {pw}_{p}^{k - 1} g_{p}^{m}}{Σ_{i = 1}^{I} i^{2} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} - Σ_{i = 1}^{I} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1}}, k &GreaterEqual; 1

式9

b^{k - 1} = \frac{Σ_{i = 1}^{I} w_{i}^{k - 1} Σ_{p = 1}^{I} p w_{p}^{k - 1} g_{p}^{m} - Σ_{i = 1}^{I} {iw}_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} g_{p}^{m}}{Σ_{i = 1}^{I} i^{2} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} - Σ_{i = 1}^{I} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1}}, k &GreaterEqual; 1

式10

其中p＝1,2…,I；

式7中，

是过零点间隔序列的最小值，

是过零点间隔序列

{d_{i}^{k}, i = 1,2, \cdot \cdot \cdot, I}

的平均值；

式8中，δ₁为权重修正因子，δ₁为任一大于0的数；

第九步，判断是否满足迭代加权最小二乘线性拟合停止条件：用

计算得到第k次迭代的加权残差平方和如果k≤K-1且

则令k＝k+1，并回到第八步，否则进入第十步；

第十步，分别通过和

计算得到信号周期斜率k₀的估计值

和起始频率参数f_l的估计值

本发明方法的第三步中，x_i(m)的离散傅里叶变换是通过快速傅里叶变换实现。

本发明方法的优选方案中，第八步中，当权重修正因子δ₁=1时，估计效果较好。

本发明方法的第四步中，通过搜索功率谱值Y_i(l₂)最大值所对应的离散频率序号L_i，然后将L_i代入式3，估计出本短时窗内数据序列的瞬时频率估计值f_i。

有益效果：与现有技术相比，本发明具有以下优点：

1．目前用于估计双曲调频信号参数的方法主要有最大似然法，非线性最小二乘法和时频分析与图像处理相结合的方法。最大似然法和非线性最小二乘法都需要解非线性方程，而对于双曲调频信号来说，这两种方法需要求解的非线性方程都不存在解析解，因此需要使用数值解法，这就加大了求解的难度和计算的复杂度。时频分析与图像处理相结合的方法为了达到较好的估计结果，需要重复使用时频滤波器，这使得该算法的计算量大大增加，限制了其工程实用性。以上三种方法都需要复杂的计算和参数搜索，计算量大，运算速度慢，而在实际工程中常需处理数据流，要求较快的运算速度，这限制了上述三种算法的工程应用。本发明的估计方法通过对瞬时频率取倒数得到过零点间隔，然后进行拟合，将复杂的非线性问题转化为线性问题，大大降低了参数估计的计算复杂度，无需进行参数搜索，可以快速实现，具有工程实用性。

2.计算每个窗内的功率谱时需要计算该短时窗内数据序列的离散傅里叶变换，而离散傅里叶变换可以通过快速傅里叶变换实现，快速傅里叶变换是一种使用十分广泛而且有效的计算工具，这大大提高了估计方法的运算效率和实用性，使其能够较好的应用于工程实际。

3.通过对{g_i,i＝1,2,…,I}进行滑动中值滤波，可以滤除{g_i,i＝1,2,…,I}中非连续出现的异常值，大大降低了该类异常值对参数估计结果的影响，提高了方法的稳健性。

4.传统的最小二乘线性拟合方法，对所有的样本点同等对待，在样本中无异常值时是最优估计方法，但是当信噪比较低，导致提取的瞬时频率中存在异常值时，会受异常值影响，参数估计性能急剧下降，而本发明的估计方法使用改进的迭代变权最小二乘线性拟合方法：一方面，估计方法不受拟合值绝对大小的影响，适用性强；另一方面，能够有效地降低估计的瞬时频率中存在的异常值对参数估计结果的影响，提高参数估计的精度，降低对信噪比的要求。

附图说明

图1为本发明方法的流程图。

图2为本发明方法实施例中双曲调频仿真信号波形图。

图3为本发明方法实施例中双曲调频仿真信号的瞬时频率曲线。

图4为本发明方法实施例中双曲调频仿真信号的过零点间隔曲线。

图5为本发明方法实施例中在信噪比为-9dB的情形下，叠加了背景噪声后的接收信号仿真波形图。

图6为本发明方法实施例中瞬时频率实际值、瞬时频率估计值和拟合曲线的比较。

具体实施方式

下面结合附图，对本发明的技术方案进行详细的说明。

如图1所示，本发明的双曲调频信号周期斜率和起始频率估计方法，包括以下步骤：

K取值为大于等于2的正整数，ε取值为小于等于0.1的正数；

X_{i} (l_{1}) = Σ_{m = 0}^{M - 1} x_{i} (m) e^{- j \frac{2 π}{M} ml}, l_{1} = 0,1 \cdot \cdot \cdot, M

式1

其中X_i(l₁)表示离散傅里叶变换的结果，j表示虚数单位，即

Y_{i} (l_{2}) = \frac{1}{M} {| X_{i} (l_{1}) |}^{2},

l₁＝l₂且l₂0,1,₂…M/2-1 式2

其中l₂为Y_i(l₂)的离散频率序号；

在第三步中，x_i(m)的离散傅里叶变换即式1，是通过快速傅里叶变换实现的，利用快速傅里叶变换可以降低算法的运算量，提高算法的计算效率；

f_i＝(L_i-1)Δf 式3

在第四步中，搜索Y_i(l₂)最大值对应的离散频率序号L_i时，取l₂＝0,1,2…M/2-1，这是因为实数据序列的离散傅里叶变换关于中心对称，因此搜索功率谱值Y_i(l₂)最大值所对应的离散频率序号L_i时，l₂可以只取前M/2个点；

g_i＝1/f_i 式4；

通过式4，将具有双曲线形式的瞬时频率f_i转换为具有直线形式的过零点间隔g_i，通过这一步处理，大大降低了双曲调频信号参数估计的复杂度；

{g_{i}^{m}, i = 1,2 \cdot \cdot \cdot, I},

的表达式为：

g_{i}^{m} = \{\begin{matrix} g_{1} & i = 1 \\ median (g_{i - 1}, g_{i}, g_{i + 1}) & 2 < i < I \\ g_{I} & i = I \end{matrix}

式5

其中median(g_i-1,g_i,g_i+1)是取g_i-1,g_i和g_i+1的中值，通过滑动中值滤波，可以滤除序列{g_i,i＝1,2,…,I}中非连续出现的异常值，可以降低该类异常值对参数估计结果的影响。

的计算过程如下：

d_{i}^{k} = | g_{i}^{m} - a^{k - 1} - b^{k - 1} i |, k &GreaterEqual; 1

式6

λ_{i}^{k} = (d_{i}^{k} - d_{\min}^{k}) / d_{mean}^{k}, k &GreaterEqual; 1

式7

w_{i}^{k} = \{\begin{matrix} 1 & k = 0 \\ \frac{1}{λ_{i}^{k} + δ_{1}} & k &GreaterEqual; 1 \end{matrix}

式8

式6中，

和

的表达式为：

a^{k - 1} = \frac{Σ_{i = 1}^{I} i^{2} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} g_{p}^{m} - Σ_{i = 1}^{I} {iw}_{i}^{k - 1} Σ_{p = 1}^{I} {pw}_{p}^{k - 1} g_{p}^{m}}{Σ_{i = 1}^{I} i^{2} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} - Σ_{i = 1}^{I} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1}}, k &GreaterEqual; 1

式9

b^{k - 1} = \frac{Σ_{i = 1}^{I} w_{i}^{k - 1} Σ_{p = 1}^{I} p w_{p}^{k - 1} g_{p}^{m} - Σ_{i = 1}^{I} {iw}_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} g_{p}^{m}}{Σ_{i = 1}^{I} i^{2} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1} - Σ_{i = 1}^{I} w_{i}^{k - 1} Σ_{p = 1}^{I} w_{p}^{k - 1}}, k &GreaterEqual; 1

式10

其中p＝1,2…,I；

式7中，

是过零点间隔序列

的最小值，

是过零点间隔序列

{d_{i}^{k}, i = 1,2, \cdot \cdot \cdot, I}

的平均值；

式8中，δ₁为权重修正因子，δ₁为任一大于0的数；

在第八步中，加权重修正因子δ₁是为了避免当

时，出现权重无限大的情况，取δ₁=1时，效果最好，物理意义最明确。当

时，

当且

差值越大

越小，通过采用上述权重：一方面，使得权重的大小与

的绝对大小无关，仅与其相对大小有关，提高了算法的适用性；另一方面，实现了拟合过程中正常样本点与异常样本点的自动区分，降低异常值对参数估计结果的影响。式9和式10中p与i均为短时窗序号，用不同的字母表示，表达不同短时窗序号的组合；

计算得到第k次迭代的加权残差平方和

如果k≤K-1且

则令k＝k+1，并回到第八步，否则进入第十步；

第十步，分别通过

和

计算得到信号周期斜率k₀的估计值

和起始频率参数f_l的估计值

本发明的双曲调频信号周期斜率和起始频率估计方法，首先使用基于短时傅里叶变换峰值搜索的方法估计双曲调频信号的瞬时频率估计值，然后对瞬时频率值取倒数得到过零点间隔，最后采取改进的变权最小二乘线性拟合的方法进行迭代计算，从而估计出双曲调频信号的周期斜率和起始频率。

本发明的原理是利用双曲调频信号的过零点间隔是时间的线性函数。因此，如果能获得双曲调频信号的瞬时频率，并对其取倒数获得双曲调频信号的过零点间隔，则可以结合线性拟合估计出信号参数。

本发明的实施例中，仿真接收信号模型为：

其中A为信号幅度，

为初始相位，τ为脉冲宽度，f_l为起始频率，k₀=(f_h-f_l)/(τf_lf_h)为周期斜率,f_h为终止频率，其值为f_h＝1/(-k₀τ+1/f_l)。仿真信号参数分别设置为：信号幅度A＝1，初始相位

脉宽τ=0.8s，采样频率f_s＝2048Hz，对应的采样点数N＝1638，起始频率f_l＝230Hz，周期斜率k₀＝2.5×10^-3，叠加零均值高斯白噪声，方差σ²的大小由信噪比SNR决定：SNR＝10log₁₀[A²/(2σ²)]。

仿真双曲调频信号波形示意图如图2所示，从图2可以看出，信号波形越来越密，这是因为仿真双曲调频信号是上调频信号，信号频率随着时间的增大而线性增大，因此波形越来越密。图3所示是仿真双曲调频信号的瞬时频率曲线，从图3可以看出，仿真双曲调频信号的瞬时频率是时间的双曲函数。图4所示是仿真双曲调频信号的过零点间隔曲线，从图4可以看出，仿真双曲调频信号的过零点间隔是时间的线性函数。图5是在信噪比为-9dB的情形下，叠加了背景噪声后的仿真接收信号波形图，从图5可以看出，信号已经完全淹没在背景噪声中。图6是每个短时窗瞬时频率实际值，瞬时频率估计值和拟合曲线比较图，从图6可以看出，瞬时频率估计值存在多个异常值，本文发明的方法得到的拟合曲线与瞬时频率实际值吻合较好。

以该仿真信号模拟接收到的受噪声污染后的采样信号x(n),n＝0,1,…,N-1，N＝1638。下面对x(n)进行周期斜率和起始频率的估计。

实施例1：首先，进行参数初始化，设置短时窗长M＝128，短时窗移动步进L＝32，权重修正因子δ₁=1，最大迭代次数门限k=10和精度控制指标ε=10^-3，计算出总的短时窗个数初始化窗移动次数i＝1,迭代次数k＝1。

然后，移动时间窗，利用最大线谱法估计每个短时窗内的数据序列的瞬时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2,…,50，如表1所示：

f₁	f₂	f₃	f₄	f₅	f₆	f₇	f₈	f₉	f₁₀
										304	128	240	288	656	656	240	240	256	256
f₁₁	f₁₂	f₁₃	f₁₄	f₁₅	f₁₆	f₁₇	f₁₈	f₁₉	f₂₀
										256	256	256	256	256	256	272	272	272	272
f₂₁	f₂₂	f₂₃	f₂₄	f₂₅	f₂₆	f₂₇	f₂₈	f₂₉	f₃₀
										272	288	288	288	288	288	816	304	304	240
f₃₁	f₃₂	f₃₃	f₃₄	f₃₅	f₃₆	f₃₇	f₃₈	f₃₉	f₄₀
										320	976	240	320	336	336	336	336	48	352
f₄₁	f₄₂	f₄₃	f₄₄	f₄₅	f₄₆	f₄₇	f₄₈	f₄₉	f₅₀
										640	976	464	368	384	384	384	400	512	416

表1利用最大线谱法估计得到的瞬时频率序列

其中f₁,f₂,f₄,f₅,f₆,f₂₇,f₃₀,f₃₂,f₃₃,f₃₉,f₄₁,f₄₂,f₄₃,f₄₉是瞬时频率估计值中的异常值。

下一步，对f_i,i＝1,2,…,50取倒数得到g_i,i＝1,2,…,50，并对g_i,i＝1,2,…,50进行滑动中值滤波得到

最后，通过改进的加权最小二乘线性拟合的方法进行迭代计算，估计出双曲调频信号的周期斜率估计值

相对误差为

起始频率估计值

{\hat{f}}_{l} = 232.42 Hz,

相对误差为

| {\hat{f}}_{l} - f_{l} | / f = 0.0105 .

实施例2：首先，进行参数初始化，设置短时窗长M＝256，短时窗移动步进L＝64，权重修正因子δ₁=1，最大迭代次数门限k=10和精度控制指标ε=10^-3，计算出总的短时窗个数

初始化窗移动次数i＝1,迭代次数k＝1。

然后，移动时间窗，利用最大线谱法估计每个短时窗内的数据序列的瞬时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2,…,25；下一步，对f_i,i＝1,2,…,25取倒数得到g_i,i＝1,2,…,25，并对g_i,i＝1,2,…,25进行滑动中值滤波得到

相对误差为

起始频率估计值

{\hat{f}}_{l} = 230.07 Hz,

相对误差为

| {\hat{f}}_{l} - f_{l} | / f_{l} = 0.0003 .

实施例3：首先进行参数初始化，设置短时窗长M＝256，短时窗移动步进L＝64，权重修正因子δ₁=1，最大迭代次数门限k=10000和精度控制指标ε=10^-6，计算出总的短时窗个数

初始化窗移动次数i＝1,迭代次数k＝1。

相对误差为

起始频率估计值

{\hat{f}}_{l} = 231.54 Hz,

相对误差为

| {\hat{f}}_{l} - f_{l} | / f = 0.0067 .

实施例4：首先，进行参数初始化，设置短时窗长M＝128，短时窗移动步进L＝32，权重修正因子δ₁=1000，最大迭代次数门限k=100和精度控制指标ε=10^-3，计算出总的短时窗个数

初始化窗移动次数i＝1,迭代次数k＝1。

然后，移动时间窗，利用最大线谱法估计每个短时窗内的数据序列的瞬时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2,…,50；下一步，对f_i,i＝1,2,…,50取倒数得到g_i,i＝1,2,…,50，并对g_i,i＝1,2,…,50进行滑动中值滤波得到

相对误差为起始频率估计值

{\hat{f}}_{l} = 224.76 Hz,

相对误差为

| {\hat{f}}_{l} - f_{l} | / f = 0.0228 .

从实施例1、实施例2、实施例3和实施例4的结果可以看出，本发明估计方法可以获得良好的估计精度，而且计算简单，计算量小，适用于高精度快速估计双曲调频信号的周期斜率和起始频率参数的场合。