一种线性调频信号调频率和起始频率估计方法

摘要

本发明公开了一种线性调频信号调频率和起始频率估计方法，该方法包括以下步骤：第一步：获取数据序列x(n),n＝0,1,2,…,N-1；第二步：参数初始化；第三步：计算第i个短时窗内的数据序列xi(m)的功率谱Yi(l2)；第四步：采用Rife插值算法估计出本短时窗内数据序列的瞬时频率估计值fi；第五步：判断是否处理完所有短时窗的数据序列：如果未处理完,返回第三步，否则转入第六步；第六步：计算第k次迭代权重第七步：判断是否满足迭代加权最小二乘线性拟合停止条件：如不满足，返回第六步，否则转入第八步；第八步：计算出调频率和起始频率参数。该方法无需进行复杂的计算和参数搜索，在保证快速估计前提下，可以提高参数估计的精度。

说明书

技术领域

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

背景技术

线性调频信号的低峰值功率、大时宽带宽积特性，使其具有较好的可压缩 性，这使得线性调频信号在声呐和雷达等领域得到了广泛的应用。初始频率和 调频率是表征线性调频信号频率特性的两个基本参数，其估计问题是一直是信 号处理领域的重要研究内容。线性调频信号是一种典型的非平稳信号，使用简 单的傅里叶变换只能获得信号的整体频谱，难以获得其局部特性。因此，估计 线性调频信号的参数常需使用时频分析方法。目前常用于估计线性调频信号参 数的时频方法有Wigner-Hough变换(WHT)、Radon-Ambiguity变换(RAT)和分 数阶傅里叶变换(FRFT)，这几种方法都利用了线性调频信号在时频域的聚集特 性。其中，WHT是将线性调频信号的参数的估计问题转换成时频图的直线搜 索问题，RAT是将其转换成模糊图中搜索过原点直线的问题，这两种方法都 需要先进行复杂的运算再进行直线搜索，运算速度非常慢，限制了其工程应用， 而且RAT法丢失了信号的初始频率信息，只适合仅对调频率感兴趣的场合。 FRFT充分利用线性调频信号在时频域的聚集特性，可以同时估计出信号的调 频率和起始频率参数，但是该方法需要对线性调频信号的调频率和起始频率参 数进行二维搜索，计算量大，而且参数估计的精度受到FRFT固有分辨率的限 制。

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

发明内容

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

技术方案：本发明的线性调频信号调频率和起始频率估计方法，包括以下 步骤：

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

第二步：参数初始化：设置短时窗长M、短时窗移动步进L、最大迭代 次数门限K和精度控制指标ε，计算出总的短时窗个数表示向下取整运算，初始化短时窗序号i＝1，所述短时窗长M取值为为2的 整数次幂且满足M＜N/2，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\left(l_1\right)=\underset{m=0}{\overset{M-1}{\Sigma }}{x}_i\left(m\right)e^{-j\frac{2\pi }{M}\mathrm{ml}},$l₁＝0,1…M-1  式1

其中x_i(l₁)表示离散傅里叶变换的结果，j表示虚数单位，即则功率谱Y_i(l₂)为

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

第四步：采用Rife插值算法估计出第i个短时窗内数据序列的瞬时频率估 计值f_i：Rife插值算法的公式为

$f_i=(L_i-1+{\delta }_0^i)\mathrm{\Delta f}$式3

其中L_i为功率谱Y_i(l)最大值对应的离散频率索引值，Δf＝f_s/M为短时窗长 度为M的离散傅里叶变换的频率分辨率，f_s为采样频率，为该短时窗内数 据序列频率与功率谱Y_i(l)最大值对应的离散频率索引值的相对偏差，其表达式 为

${\delta }_0^i=\left(\begin{array}{cc}\frac{Y_i(L_i+1)}{Y_i(L_i+1)+Y_i\left(L_i\right)}& Y_i(L_i+1)>Y_i(L_i-1)\\ \frac{{-Y}_i(L_i-1)}{Y_i(L_i-1)+Y_i\left(L_i\right)}& Y_i(L_i+1)<Y_i(L_i-1)\\ 0& Y_i(L_i+1)=Y_i(L_i-1)\end{array}\right)$式4

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

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

$w_i^k=\left(\begin{array}{cc}1& k=1\\ \frac{1}{|f_i-f_i^{k-1}|+{\sigma }_1}& k>1\end{array}\right)$式5

其中，σ₁为权重修正因子，σ₁为任一大于0的数，为第k-1次迭代过 程中第i个时间窗内数据序列的瞬时频率拟合值，通过对所述第四步中得到的 瞬时频率估计值f_i进行加权最小二乘线性拟合得到，即：

$k_i^{k-1}=b^{k}i+a^k$

其中a^k,b^k通过在表达式中分别对a^k,b^k求导并令导数为 0，并写成矩阵形式为$\left(\begin{array}{cc}\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k& \underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k\\ \underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k& \underset{i=0}{\overset{I-1}{\Sigma }}{w}_i^k\end{array}\right)\left(\begin{array}{c}{b}^k\\ a^k\end{array}\right)=\left(\begin{array}{c}\underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k{f}_i\\ \underset{i=1}{\overset{I}{\Sigma }}{w}_i^k{f}_i\end{array}\right)$

对上式求解即可得到

$a^k=\frac{\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k{f}_p-\underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k\underset{p=1}{\overset{I}{\Sigma }}{\mathrm{pw}}_p^w{f}_p}{\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k-\underset{i=1}{\overset{I}{\Sigma }}{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k}$式6

$b^k=\frac{\underset{i=1}{\overset{I}{\Sigma }}{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{\mathrm{pw}}_p^k{f}_p-\underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k{f}_p}{\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k-\underset{i=1}{\overset{I}{\Sigma }}{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k}$式7

其中p＝1,2,…I；

第七步：判断是否满足迭代加权最小二乘线性拟合停止条件：用 计算得到第k次迭代的加权残差平方和如果k≤K-1 且则令k=k+1，并回到第六步，否则进入第八步；

第八步：计算出信号调频率和起始频率参数：分别通过计算得到起始频率f_l的估计值和调频率μ的估计值

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

本发明方法的第四步中，首先搜索功率谱值Y_i(l₂)最大值，通过所述搜索 功率谱值Y_i(l₂)最大值的位置确定离散频率索引值L_i，然后根据公式4确定相对 偏差并代入式3，估计出本短时窗内数据序列的瞬时频率估计值f_i。其中， 搜索功率谱值Y_i(l₂)最大值时取l₂＝0,1,2…M/2-1。

本发明方法的第六步中，当σ₁=1时，估计效果较好。

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

1．目前常用于估计线性调频信号参数的有Wigner-Hough变换(WHT)、 Radon-Ambiguity变换(RAT)和分数阶傅里叶变换(FRFT)。其中，WHT是将线 性调频信号的参数的估计问题转换成时频图中直线搜索的问题，RAT是将其 转换成模糊图中搜索过原点直线的问题，这两种方法都需要先进行复杂的运算 再进行直线搜索，而且RAT法丢失了信号的起始频率信息；FRFT充分利用线 性调频信号在分数阶域Fourier的聚集特性，在已知信号起始时间的条件下， 可以同时估计出信号的调频率和起始频率参数，但该方法需要对线性调频信号 的调频率和起始频率参数进行二维搜索，而且参数估计精度受到FRFT固有分 辨率的限制。以上三种方法虽然可以获得高精度的调频率估计，但都需要复杂 的计算和参数搜索，计算量大，运算速度慢，而在实际工程中常需处理数据流， 要求较快的运算速度，这限制了上述三种算法的工程应用。本发明的估计方法 无需进行复杂的计算和参数搜索，可以快速实现，具有工程实用性。

2．本发明的估计方法将Rife插值算法应用于短时傅里叶变换，能够提高 短时傅里叶变换中每段短时窗内频率估计的精度，相对于仅用短时傅里叶变换， 能够更准确地估计线性调频信号的瞬时频率。

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

附图说明

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

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

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

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

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

具体实施方式

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

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

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

第二步：参数初始化：设置短时窗长M、短时窗移动步进L、最大迭代 次数门限K和精度控制指标ε，计算出总的短时窗个数表示向下取整运算，初始化短时窗序号i＝1，所述短时窗长M取值为为2的 整数次幂且满足M＜N/2，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\left(l_1\right)=\underset{m=0}{\overset{M-1}{\Sigma }}{x}_i\left(m\right)e^{-j\frac{2\pi }{M}\mathrm{ml}},$l₁＝0,1…M-1  式1

其中x_i(l1)表示离散傅里叶变换的结果，j表示虚数单位，即则功率谱Y_i(l₂)为

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

其中，离散傅里叶变换可以通过快速傅里叶变换实现；

第四步：采用Rife插值算法估计出第i个短时窗内数据序列的瞬时频率估 计值f_i：Rife插值算法的公式为

$f_i=(L_i-1+{\delta }_0^i)\mathrm{\Delta f}$式3

其中L_i为功率谱Y_i(l)最大值对应的离散频率索引值，Δf＝f_s/M为短时窗长 度为M的离散傅里叶变换的频率分辨率，f_s为采样频率，为该短时窗内数 据序列频率与功率谱Y_i(l)最大值对应的离散频率索引值的相对偏差，其表达式 为

${\delta }_0^i=\left(\begin{array}{cc}\frac{Y_i(L_i+1)}{Y_i(L_i+1)+Y_i\left(L_i\right)}& Y_i(L_i+1)>Y_i(L_i-1)\\ \frac{{-Y}_i(L_i-1)}{Y_i(L_i-1)+Y_i\left(L_i\right)}& Y_i(L_i+1)<Y_i(L_i-1)\\ 0& Y_i(L_i+1)=Y_i(L_i-1)\end{array}\right)$式4

在第四步中，Rife插值算法属于现有技术，例如《电子学报》2004年， 32(4)的第625页至第628页中公开的内容，本发明方法将其应用于短时傅里 叶变换，可以提高每一短时窗内数据序列瞬时频率估计值的精度。

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

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

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

$w_i^k=\left(\begin{array}{cc}1& k=1\\ \frac{1}{|f_i-f_i^{k-1}|+{\sigma }_1}& k>1\end{array}\right)$式5

其中，σ₁为权重修正因子，σ₁为任一大于0的数，为第k-1次迭代过 程中第i个时间窗内数据序列的瞬时频率拟合值，通过对所述第四步中得到的 瞬时频率估计值f_i进行加权最小二乘线性拟合得到，即：

$k_i^{k-1}=b^{k}i+a^k$

其中a^k,b^k通过在表达式中分别对a^k,b^k求导并令导数为 0，并写成矩阵形式为$\left(\begin{array}{cc}\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k& \underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k\\ \underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k& \underset{i=0}{\overset{I-1}{\Sigma }}{w}_i^k\end{array}\right)\left(\begin{array}{c}{b}^k\\ a^k\end{array}\right)=\left(\begin{array}{c}\underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k{f}_i\\ \underset{i=1}{\overset{I}{\Sigma }}{w}_i^k{f}_i\end{array}\right)$

对上式求解即可得到

$a^k=\frac{\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k{f}_p-\underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k\underset{p=1}{\overset{I}{\Sigma }}{\mathrm{pw}}_p^w{f}_p}{\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k-\underset{i=1}{\overset{I}{\Sigma }}{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k}$式6

$b^k=\frac{\underset{i=1}{\overset{I}{\Sigma }}{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{\mathrm{pw}}_p^k{f}_p-\underset{i=1}{\overset{I}{\Sigma }}{\mathrm{iw}}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k{f}_p}{\underset{i=1}{\overset{I}{\Sigma }}{i}^2{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k-\underset{i=1}{\overset{I}{\Sigma }}{w}_i^k\underset{p=1}{\overset{I}{\Sigma }}{w}_p^k}$式7

其中p＝1,2,…I；

在第六步中，加权重修正因子σ₁是为了避免当时，出现权重无限 大的情况，取σ₁＝1时，效果最好，物理意义最明确。当时当 且与f_i差值越大越小，通过采用上述权重实现了拟 合过程中正常样本点与异常样本点的自动区分，降低异常值对参数估计结果的 影响。式6和式7中p与i均为窗序号，用不同的字母表示，表达不同短时窗 序号的组合。

第七步：判断是否满足迭代加权最小二乘线性拟合停止条件：用 计算得到第k次迭代的加权残差平方和如果k≤K-1 且则令k=k+1，并回到第六步，否则进入第八步；

第八步：计算出信号调频率和起始频率参数：分别通过计算得到起始频率f_l的估计值和调频率μ的估计值

本发明的线性调频信号调频率和起始频率估计方法，首先使用基于Rife 插值算法改进的短时傅里叶变换估计线性调频信号的瞬时频率估计值，然后采 取变权最小二乘线性拟合的方法进行迭代计算，从而估计出线性调频信号的起 始频率和调频率参数。

本发明的原理是利用线性调频信号的瞬时频率是时间的线性函数。因此， 如果能获得线性调频信号的瞬时频率，则可以结合线性拟合估计出信号参数。

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

其中A为信号幅度，为初始相位，τ为脉冲宽度，f_l为起始频率，μ＝B/τ 为调频率,B为脉冲带宽。仿真信号参数分别设置为：信号幅度A=1，初始相位 脉宽τ＝20ms，采样频率f_s＝70kHz，对应的采样点数N＝1400，起始频率 f_l＝5kHz，调频率μ＝500KHz/s，叠加零均值高斯白噪声，方差σ²的大小由信噪 比SNR决定：SNR＝10lg(A²/2σ²)。

仿真线性调频信号波形示意图如图2所示，从图2可以看出，信号波形越 来越密，这是因为仿真线性调频信号是上调频信号，信号频率随着时间的增大 而线性增大，因此波形越来越密。图3所示是仿真线性调频信号的瞬时频率曲 线，从图3也可看出，仿真线性调频信号的瞬时频率是时间的线性函数。图4 是在信噪比为-9dB的情形下，叠加了背景噪声后的仿真接收信号波形图。图5 是每个短时窗瞬时频率估计值，拟合直线和瞬时频率实际值的比较图，从图5 中可以看出瞬时频率估计值存在异常值，本文发明的方法得到的拟合直线与瞬 时频率实际值吻合较好。以该仿真信号模拟接收到的受噪声污染后的采样信号 x(n),n＝0,1,...,N-1，N＝1400。下面对x(n)进行调频率和起始频率的估计。

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

然后移动时间窗，利用Rife插值算法估计每个短时窗内的数据序列的瞬 时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2，…,20，如表1所示

表1利用Rife插值算法估计得到的瞬时频率序列

其中f₃，f₄，f₁₁，f₁₅，f₁₈，f₁₉是瞬时频率估计值中存在的异常值。

最后通过加权最小二乘线性拟合的方法进行迭代计算，估计出线性调频信 号的调频率估计值相对误差为起始频率估计值 ${\hat{f}}_l=5112.6\mathrm{Hz},$相对误差为$|{\hat{f}}_l-f_l|/f_l=0.0225.$

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

然后移动时间窗，利用Rife插值算法估计每个短时窗内的数据序列的瞬 时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2，…,20；

最后通过加权最小二乘线性拟合的方法进行迭代计算，估计出线性调频信 号的调频率估计值相对误差为起始频率估计值 ${\hat{f}}_l=5097.6\mathrm{Hz},$相对误差为$|{\hat{f}}_l-f_l|/f_l=0.01952.$

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

然后移动时间窗，利用Rife插值算法估计每个短时窗内的数据序列的瞬 时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2，…,9；

最后通过加权最小二乘线性拟合的方法进行迭代计算，估计出线性调频信 号的调频率估计值相对误差为起始频率估计值 ${\hat{f}}_l=5121.3\mathrm{Hz},$相对误差为$|{\hat{f}}_l-f_l|/f_l=0.02426.$

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

然后移动时间窗，利用Rife插值算法估计每个短时窗内的数据序列的瞬 时频率估计值，得到瞬时频率估计值序列f_i,i＝1,2，…,9；

最后通过加权最小二乘线性拟合的方法进行迭代计算，估计出线性调频信 号的调频率估计值相对误差为起始频率估计值 ${\hat{f}}_l=5081.1\mathrm{Hz},$相对误差为$|{\hat{f}}_l-f_l|/f_l=0.01622.$

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