CN113515725B - 一种基于参数预估计的改进径向高斯核时频分析方法 - Google Patents

Abstract

本发明公开了一种基于参数预估计的改进径向高斯核时频分析方法，该方法包括如下步骤：1、获取待处理脉冲信号；2、计算待处理脉冲信号直角坐标系下的模糊函数；3、预估信号模糊域自项径向角；4、将直角坐标系下模糊函数转换为极坐标系下模糊函数；5、利用Step‑Project算法计算极坐标系下信号最优核函数扩展矢量；6、计算直角坐标系下信号最优核函数；7、利用信号最优核函数对信号模糊函数滤波，将滤波后的结果变换到时频域得到信号的时频分析结果。本发明是对传统的基于径向高斯核时频分析方法的改进，可在保持其算法性能的基础上，提高其运算效率，适合对未知分量数的线性调频和双曲调频信号等进行高时效时频分析。

Description

一种基于参数预估计的改进径向高斯核时频分析方法

技术领域

本发明属于信号处理技术领域，尤其涉及一种基于参数预估计的改进径向高斯核时频分析方法。

背景技术

大多数现实生活中的信号具有非平稳特性，对于时变或频变的非平稳信号，传统的。一维时间或频率分析无法全面地描述出信号的本质，即时间和频率无法同时联合描述信号，不能反映出非平稳信号的统计量时间变化特征。因此，能够同时在时间和频率这两个维度上表示信号能量的时频分析方法在信号处理应用中有着重要意义。目前，时频分析方法主要分为线性时频分析方法和二次型时频分析方法。

最为经典的线性时频分析方法是短时傅里叶变换，它的优点是计算简单，运算量小，且不存在交叉项干扰问题，但是其窗函数的形状和宽度大大地限制了该方法的时频分辨率，其时频分辨率无法同时提高，因此存在较大的应用局限性。

二次型时频分析方法建立在Cohen的双线性时频分布理论的基础上，即可以认为一般的二次型时频分布是利用核函数对于Wigner-Ville分布的二位滤波结果。Wigner-Ville分析方法在一定程度上弥补了短时傅里叶变换在分辨率上的缺点，该方法是单分量信号时频分析的最佳方法，但是当它被用于多分量信号的分析上时，会产生严重的交叉项干扰问题。对于使用核函数的时频分析方法，选用不同的核函数可以得到不同性能的时频性能。一种基于信号自适应的径向高斯核函数时频分析方法选用了径向高斯函数作为其核函数，根据信号的模糊域特征自适应地调整核函数的分布形状，可以有效抑制交叉项干扰，保留自项的能量集中，同时具有良好的时频分辨率。

径向高斯核函数时频分析方法是一种性能优越的时频分析方法，尤其适用于分析多分量的非平稳信号。但是，由于它在解决关于径向高斯函数中的扩展函数σ(ψ)的优化问题时，采取了在整个模糊域上半平面即(0,π)角度范围内迭代搜索的方法，因此在处理实际的大量观测信号数据时会面临过大的运算量问题，导致计算时间过久而不适用于实时的处理。

发明内容

本发明目的在于提供一种基于参数预估计的改进径向高斯核时频分析方法,以解决模糊函数过原点且集中于原点的线性调频、双曲调频信号，减少优化过程中的扩展矢量迭代搜索的维数，提高优化过程迭代搜索扩展矢量的效率的技术问题。

为解决上述技术问题，本发明的具体技术方案如下：

一种基于参数预估计的改进径向高斯核时频分析方法，包括如下步骤：

步骤1、获取待处理脉冲信号的采样数据序列s(a)；

步骤2、计算待处理脉冲信号采样数据序列s(a)直角坐标系下的模糊函数A_z(m,n)；

步骤3、在模糊域内检测信号自项并估计信号自项的径向角度φ_k；

步骤4、对待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)进行插值，得到待处理脉冲信号极坐标系下的模糊函数A_p(p,q)；

步骤5、利用Step-Project算法计算待处理脉冲信号的采样数据序列在极坐标系下对应模糊域上半平面的最优径向高斯核函数扩展矢量σ_opt；

步骤6、对极坐标系下的扩展矢量σ进行插值，转换为直角坐标系下的扩展矢量，并计算直角坐标系下的核函数Φ_d(m,n)；

步骤7、将最优核函数Φ_d(m,n)与信号的模糊函数A_z(m,n)相乘，得到模糊域滤波后的特征函数A_F(m,n)，并将A_F(m,n)变换到时频域得到信号的时频分析结果ρ(a,b)。

进一步的，步骤1中，采用如下方法获取处理的脉冲信号采样数据序列：

从传感器接收L个采样点的实时采集数据作为待处理脉冲信号的采样数据序列s(a)，已知采样频率f_s，所述的L为包含整个多脉冲信号脉宽长度所对应的采样点个数，取值为2的正整数次幂，且L≥4。

进一步的，步骤1中，采用如下方法获取处理的脉冲信号采样数据序列：从存储器中提取包含整个多脉冲信号的L个采样点数据作为待处理的数据序列s(a)，已知采样频率f_s，所述的L为包含整个多脉冲信号脉宽长度所对应的采样点个数，取值为2的正整数次幂，且L≥4。

进一步的，步骤2中，计算信号待处理脉冲信号的采样数据序列s(a)直角坐标系下的模糊函数A_z(m,n)的方法，包括以下步骤：

步骤2.1、计算待处理脉冲信号采样数据序列s(a)的离散傅里叶变换S(α)，即：

其中，α为S(α)的离散频率索引，j表示虚数单位，即

利用待处理脉冲信号的离散傅里叶变换S(α)生成s(a)的离散解析信号z(a)的傅里叶变换Z(β)，β＝0,1,…L-1，即：

其中，β为Z(β)的离散频率索引；

对Z(β)进行离散逆傅里叶变换得到待处理脉冲信号的离散解析信号z(a)，即：

步骤2.2，计算离散解析信号z(a)的自相关函数K_z(a,n)，即：

其中，

表示

的共轭形式，n代表离散延时量，n＝0,1,…L-1，a＝0,1,…L-1；

步骤2.3，依据离散解析信号z(a)的自相关函数K_z(a,n)计算待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)，即：

其中，m代表离散频移，n代表离散延时量。

进一步的，步骤3中，在模糊域检测信号自项并估计信号自项的径向角度具体包括以下步骤：

步骤3.1、计算待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)在径向距离为0处的离散Radon变换R₀(θ_x)，即：

其中，R₀(θ_x)代表径向距离为0处的Radon变换，θ_x代表Radon变换的离散径向角度，其取值为

x为整数，且0≤x≤X-1，X为总的离散径向角度，X取值为≥2的正整数；δ(·)代表冲激函数；

步骤3.2、估计信号自项在模糊域的径向角度φ_k：

对R₀(θ_x)进行去趋势项处理，得到其去趋势项分量R_d,0(θ_x)，即：

R_d,0(θ_x)＝R₀(θ_x)-R_t,0(θ_x)，

其中，R_t,0(θ_x)代表R₀(θ_x)的趋势项分量，其取值为

其中，B为观测矩阵，为单位矩阵，即：

为采用正则化最小二乘法估计得到的参数，即：

其中，λ为正则化参数，且为大于0的实数，正则化参数λ取值为10，D₂为二阶导数算子的离散形式，D₂是一个(L-2)×(L-1)的二维矩阵，即：

步骤3.3、在得到去趋势项分量R_d,0(θ_x)后，从去趋势项分量R_d,0(θ_x)中搜索所有峰值点(φ_e,R(φ_e))，e＝1,2,...,E，其中φ_e为峰值点所对应的径向角，E为总的峰值点数，为正整数，R(φ_e)表示径向角φ_e所对应的去趋势项分量R_d,0(φ_e)的值；

定义自项的检测门限值R_th为

其中，μ是所有搜索峰值的平均值，即

从R_d,0(θ_x)中搜索出所有大于R_th的径向角，记为φ_k，k＝1,2,...,K，其中K为搜索出的大于R_th的径向角的总个数，即检测到的自项分量的总个数，K≤E。

进一步的，步骤4中，采用如下方法对待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)进行插值，转换待处理脉冲信号为极坐标系下的模糊函数A_p(p,q)：

步骤4.1、确定直角坐标系的两个坐标轴对应矢量的取值，横轴对应的矢量为v，纵轴对应的矢量为τ，它们均为长度为L的矢量，分别为

τ＝[-LU,…,(L-1)U]^T，

其中，上标中的T代表矩阵转置符号，矩阵中的U为与离散间隔相关的值，

步骤4.2，生成直角坐标系对应的坐标区域，横坐标区域为V，纵坐标区域为Τ，它们均为L×L维矩阵，分别为

Τ_L×L＝[τ … τ]；

步骤4.3，确定极坐标系的两个坐标轴对应矢量的取值，极角对应的矢量为ψ，极径对应的矢量为ρ；

极角矢量为

极径矢量为

其中，ψ的采样间隔为

ρ的采样间隔为

Q为ψ的采样点数，Q＝L，P为ρ的采样点数，

round(·)为取整函数；

步骤4.4、生成极坐标系对应的坐标区域，横坐标区域为X，纵坐标区域为Y，它们均为P×Q维矩阵，分别为

X＝ρ^T×cos(ψ)，

Y＝ρ^T×sin(ψ)；

步骤4.5、运用双线性插值运算将待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)转换为极坐标系下的模糊函数A_p(p,q)，即：

其中，A_p(p,q)为待插值点，{A_z(m,n),A_z(m+1,n),A_z(m,n+1),A_z(m+1,n+1)}为A_z(m,n)上最小包围点A_p(p,q)的相邻的四个模糊函数点值，d_v,min(p,q)为点A_p(p,q)到点A_z(m,n)和点A_z(m,n+1)所在直线的距离，d_h,min(p,q)为点A_p(p,q)到点A_z(m,n)和点A_z(m+1,n)所在直线的距离。

进一步的，步骤5中，利用Step-Project算法计算待处理脉冲信号直角坐标系下对应模糊域上半平面的最优径向高斯核函数的扩展矢量σ_opt，具体包括以下步骤：

步骤5.1、确定扩展矢量σ_q的搜索维数R，信号自项径向角对应的极角矢量的索引为

其中，

表示使得括号内取最小值的q值；

扩展矢量σ_q中，对于每一个自项径向角对应的极角矢量索引，定义极角矢量索引的单边搜索索引宽度为N_s，

其中round()代表四舍五入取整数运算，则需要搜索的索引范围为

S_k＝{l_k-N_s,l_k-N_s+1,…,l_k,…l_k+N_s-1,l_k+N_s}。

扩展矢量σ_q中需要搜索的总索引范围为

记S的元素个数为R，即扩展矢量σ_q的搜索维数；

步骤5.2、确定Step-Project算法中要使用的各元素的初值：

μ(0)＝1，μ_m＝1001，

ε₁＝0.1，ε₂＝1-ε₁，

Δ＝10，ε₃＝1×10^-5，

κ＝0；

其中，γ为核函数的体积限制参数，σ(0)为扩展矢量的初始值，μ(0)为初始步长，μ_m为最大步长，ε₁和ε₂为测试步长大小的标准值，ε₃为迭代终止误差中设置的参数，Δ为步长倍增或倍减因子，κ为迭代次数；

步骤5.3、开始进行Step-Project算法，子步骤如下：

步骤5.3.1、求梯度向量

计算目标函数f(σ_q(κ))的值，即

其中，σ_q(κ)为第κ次迭代时的扩展矢量的值，

代表求梯度，

代表f对σ_q(κ)求偏导数；

步骤5.3.2、更新下一次迭代的扩展矢量和优化目标，即

步骤5.3.3、根据f(σ_q(κ))与f(σ_q(κ-1))测试步长μ(κ)，若过小则令步长倍增，若超过最大步长μ_m则停止迭代，转步骤5.4，若在μ_m以内过大则令步长倍减；定义判断测试步长大小的依据g为优化目标实际上升值与预期上升值的比值，即

若g＜ε₁，则令μ(κ+1)＝μ(κ)/Δ，σ_q(κ+1)＝σ_q(κ)；若g＞ε₂，则令

若ε₁≤g≤ε₂，则设置两个迭代终止误差为

∈₂＝(f(σ_q(κ+1))-f(σ_q(κ)))-ε₃·(1+f(σ_q(κ)))，

且迭代次数递增，即令κ＝κ+1；

步骤5.3.4、判断若迭代终止误差值均为非正，或迭代次数κ超过最大迭代次数50，则停止迭代，输出σ_q(κ)，记为σ_q'，否则转步骤5.3.1；

步骤5.4：确定完整的最优径向高斯核函数的扩展矢量σ_opt，首先预置

再代入有值部分，即

σ_opt(S(i))＝σ_q'(i),i＝1,2,…,R。

进一步的，所述步骤6中，对极坐标系下的扩展矢量插值转换为直角坐标下的扩展矢量，具体包括以下步骤：

步骤6.1：将ψ和σ_opt补充为整个平面上的表示，即

ψ＝[ψ,ψ(i)+π],i＝2,3,…,Q，

σ＝[σ_opt,σ_opt(i)],i＝2,3,…,Q；

步骤6.2：依据直角坐标轴V和Τ计算它们对应于极坐标中的角度，即

且需满足0＜Ψ_d＜2π；

步骤6.3：首先将直角坐标点的角度转换为一个列向量，即Ψ_d'，再对极坐标系中的扩展矢量进行线性插值得到每个直角坐标点对应的扩展矢量值，即

其中，j满足如下约束关系：

j＝argmaxΨ(j)≤Ψ_d'(i)。

对σ_d进行按列重排，变换为L×L维的矩阵，即得到直角坐标系下的扩展矢量，仍记为σ_d；

步骤6.4：计算直角坐标系下的最优核函数Φ_opt，即

Φ_opt记为Φ_opt(m,n),m＝0,1,…,L-1,n＝0,1,…,L-1。

进一步的，步骤7中，模糊域滤波后的特征函数A_F(m,n)为

A_F(m,n)＝A_z(m,n)·Φ_opt(m,n)，m＝0,1,…,L-1，n＝0,1,…,L-1；

信号的时频分析结果ρ(a,b)为

其中，a代表离散频率索引，b代表离散频率索引。

进一步的，步骤7中，模糊域到时频域的转换用快速傅里叶变换实现。

本发明的一种基于参数预估计的改进径向高斯核时频分析方法，具有以下优点：

(1)本发明由于是建立在传统的径向高斯核函数时频分析方法的基础上的，因此可以具备与其完全一致的时频分析性能，即主要适用于模糊函数过原点且集中于原点的线性调频、双曲调频等信号，能够不遗漏地提取多分量信号的分量信息，完整地表示信号的时频分布，在保持自项能量集中的同时，有效地抑制交叉项的干扰。特别地，在对实际试验得到的非平稳信号进行分析时，虽然引入了传统径向高斯核函数时频分析结果中没有的噪声干扰，但是由于相对能量集中部分可以忽略不计，因此认为本发明的时频分析性能较好。

(2)本发明在核函数计算过程中，首先对信号的分量和自项径向角度进行估计，确认分量角度估计值后大大缩小了扩展矢量的求解维数，从而减少了迭代搜索的运算时间。经过验证得出，迭代搜索的效率得到了有效的提高，因此是一种高效的时频分析方法。

附图说明

图1为本发明的一种基于参数预估计的改进径向高斯核时频分析方法的流程示意图；

图2为本发明第一实施例计算得到的最优核函数的模糊域图像；

图3为本发明第一实施例仿真信号的时频图像；

图4为本发明第二实施例计算得到的最优核函数的模糊域图像；

图5为本发明第二实施例仿真信号的时频图像；

图6为本发明第三实施例计算得到的最优核函数的模糊域图像；

图7为本发明第三实施例仿真信号的时频图像；

具体实施方式

为了更好地了解本发明的目的、结构及功能，下面结合附图，对本发明一种基于参数预估计的改进径向高斯核时频分析方法做进一步详细的描述。

如图1所示，一种基于参数预估计的改进径向高斯核时频分析方法，包括如下步骤：

步骤1、获取待处理脉冲信号的采样数据序列：从传感器接收L个采样点的实时采集数据作为待处理脉冲信号的采样数据序列s(a)，a＝0,1,…L-1；或从存储器中提取包含整个多脉冲信号的L个采样点数据作为待处理脉冲信号的采样数据序列s(a)，已知采样频率f_s，所述的L为包含整个多脉冲信号脉宽长度所对应的采样点个数，取值为2的正整数次幂，且L≥4。

步骤2、计算待处理脉冲信号采样数据序列s(a)直角坐标系下的模糊函数A_z(m,n)，具体包括以下步骤：

步骤2.1、计算待处理脉冲信号采样数据序列s(a)的离散傅里叶变换S(α)，即：

其中，α为S(α)的离散频率索引，j表示虚数单位，即

利用离散傅里叶变换S(α)生成s(a)的离散解析信号z(a)的傅里叶变换Z(β)，β＝0,1,…L-1，即：

其中，β为Z(β)的离散频率索引；

对Z(β)进行离散逆傅里叶变换得到待处理脉冲信号的离散解析信号z(a)，即：

步骤2.2、计算离散解析信号z(a)的自相关函数K_z(a,n)，即：

其中，

表示

的共轭形式，n代表离散延时量，n＝0,1,…L-1，a＝0,1,…L-1；

步骤2.3、依据离散解析信号z(a)的自相关函数K_z(a,n)计算待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)，即：

其中，m代表离散频移，n代表离散延时量。

步骤3、在模糊域内检测信号自项并估计信号自项的径向角度φ_k，具体包括以下步骤：

步骤3.1、计算待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)在径向距离为0处的离散Radon变换R₀(θ_x)，即：

其中，R₀(θ_x)代表径向距离为0处的Radon变换，θ_x代表Radon变换的离散径向角度，其取值为

x为整数，且0≤x≤X-1，X为总的离散径向角度，X取值为大于等于2的正整数；δ(·)代表冲激函数；

步骤3.2、估计信号自项在模糊域的径向角度φ_k：

对R₀(θ_x)进行去趋势项处理，得到其去趋势项分量R_d,0(θ_x)，即：

R_d,0(θ_x)＝R₀(θ_x)-R_t,0(θ_x)，

其中，R_t,0(θ_x)代表R₀(θ_x)的趋势项分量，其取值为

其中，B为观测矩阵，为单位矩阵，即：

为采用正则化最小二乘法估计得到的参数，即：

其中，λ为正则化参数，且为大于0的实数，正则化参数λ取值为10，D₂为二阶导数算子的离散形式，D₂是一个(L-2)×(L-1)的二维矩阵，即：

步骤3.3、在得到去趋势项分量R_d,0(θ_x)后，从去趋势项分量R_d,0(θ_x)中搜索所有峰值点(φ_e,R(φ_e))，e＝1,2,...,E，其中，φ_e为峰值点所对应的径向角，E为总的峰值点数，为正整数，R(φ_e)表示径向角φ_e所对应的去趋势项分量R_d,0(φ_e)的值；

定义自项的检测门限值R_th为

其中，μ是所有搜索峰值的平均值，即

从R_d,0(θ_x)中搜索出所有大于R_th的径向角，记为φ_k，k＝1,2,...,K，其中K为搜索出的大于R_th的径向角的总个数，即检测到的自项分量的总个数，K≤E。

步骤4、对待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)进行插值，转换待处理脉冲信号为极坐标系下的模糊函数A_p(p,q)，具体包括以下步骤：

步骤4.1、确定直角坐标系的两个坐标轴对应矢量的取值，横轴对应的矢量为v，纵轴对应的矢量为τ，它们均为长度为L的矢量，分别为

τ＝[-LU,…,(L-1)U]^T，

其中，上标中的T代表矩阵转置符号，矩阵中的U为与离散间隔相关的值，

步骤4.2、生成直角坐标系对应的坐标区域，横坐标区域为V，纵坐标区域为Τ，它们均为L×L维矩阵，分别为

Τ_L×L＝[τ … τ]；

步骤4.3、确定极坐标系的两个坐标轴对应矢量的取值，极角对应的矢量为ψ，极径对应的矢量为ρ。

极角矢量和极径矢量分别为

其中，ψ的采样间隔为

ρ的采样间隔为

Q为ψ的采样点数，Q＝L，P为ρ的采样点数，

round(·)为取整函数；

步骤4.4、生成极坐标系对应的坐标区域，横坐标区域为X，纵坐标区域为Y，它们均为P×Q维矩阵，分别为

X＝ρ^T×cos(ψ)，

Y＝ρ^T×sin(ψ)；

步骤4.5、运用双线性插值运算将直角坐标系对应的坐标区域中的信号模糊函数采样A_z(m,n)转换为极坐标系对应的坐标区域中的信号模糊函数A_p(p,q)，即：

其中，A_p(p,q)为待插值点，{A_z(m,n),A_z(m+1,n),A_z(m,n+1),A_z(m+1,n+1)}为A_z(m,n)上最小包围点A_p(p,q)的相邻的四个模糊函数点值，d_v,min(p,q)为点A_p(p,q)到点A_z(m,n)和点A_z(m,n+1)所在直线的距离，d_h,min(p,q)为点A_p(p,q)到点A_z(m,n)和点A_z(m+1,n)所在直线的距离。

步骤5、利用Step-Project算法计算待处理脉冲信号直角坐标系下对应模糊域上半平面的最优径向高斯核函数的扩展矢量σ_opt，具体包括如下步骤：

步骤5.1、确定扩展矢量σ_q的搜索维数R。信号自项径向角对应的极角矢量的索引为

其中，

表示使得括号内取最小值的q值。

扩展矢量σ_q中，对于每一个自项径向角对应的极角矢量索引，定义极角矢量索引的单边搜索索引宽度为N_s，

其中round()代表四舍五入取整数运算，则需要搜索的索引范围为

S_k＝{l_k-N_s,l_k-N_s+1,…,l_k,…l_k+N_s-1,l_k+N_s}。

扩展矢量σ_q中需要搜索的总索引范围为

记S的元素个数为R，即扩展矢量σ_q的搜索维数；

步骤5.2、确定Step-Project算法中要使用的各元素的初值。

μ(0)＝1，μ_m＝1001，

ε₁＝0.1，ε₂＝1-ε₁，

Δ＝10，ε₃＝1×10^-5，

κ＝0；

其中，γ为核函数的体积限制参数，σ(0)为扩展矢量的初始值，μ(0)为初始步长，μ_m为最大步长，ε₁和ε₂为测试步长大小的标准值，ε₃为迭代终止误差中设置的参数，Δ为步长倍增或倍减因子，κ为迭代次数。α的理想取值区间为α∈[2,5]。若信号在模糊域中的交叉项距离原点特别远，则可以考虑α＞5以获得更好的抑制交叉项的效果；

步骤5.3、开始进行Step-Project算法，子步骤如下：

步骤5.3.1、求梯度向量

计算目标函数f(σ(κ))的值，即

其中，σ(κ)为第κ次迭代时的扩展矢量的值，

代表求梯度，

代表f对σ_q(κ)求偏导数。

步骤5.3.2、更新下一次迭代的扩展矢量和优化目标，即

步骤5.3.3、根据f(σ(κ))与f(σ(κ-1))测试步长μ(κ)是否合适，若过小则令步长倍增，若超过最大步长μ_m则停止迭代，转步骤5.4，若在μ_m以内过大则令步长倍减。定义判断测试步长大小的依据g为优化目标实际上升值与预期上升值的比值，即

若g＜ε₁，则令μ(κ+1)＝μ(κ)/Δ，σ(κ+1)＝σ(κ)；若g＞ε₂，则令

若ε₁≤g≤ε₂，则设置两个迭代终止误差为

∈₂＝(f(σ(κ+1))-f(σ(κ)))-ε₃·(1+f(σ(κ)))，

且迭代次数递增，即令κ＝κ+1；

步骤5.3.4、判断若迭代终止误差值均为非正，或迭代次数κ超过最大迭代次数50，则停止迭代，输出σ(κ)，记为σ_q'，否则转步骤5.3.1。

步骤5.4：确定完整的最优径向高斯核函数的扩展矢量σ_opt，首先预置

再代入有值部分，即

σ_opt(S(i))＝σ_q'(i),i＝1,2,…,R。

步骤6、对极坐标系下的扩展矢量σ进行插值，转换为直角坐标系下的扩展矢量，并计算直角坐标系下的核函数Φ_d(m,n)，具体包括以下步骤：

步骤6.1、将ψ和σ_opt补充为整个平面上的表示，即

ψ＝[ψ,ψ(i)+π],i＝2,3,…,Q，

σ＝[σ_opt,σ_opt(i)],i＝2,3,…,Q；

步骤6.2、依据直角坐标轴V和Τ计算它们对应于极坐标中的角度，即

且需满足0＜Ψ_d＜2π；

步骤6.3、首先将直角坐标点的角度转换为一个列向量，即Ψ_d'，再对极坐标系中的扩展矢量进行线性插值得到每个直角坐标点对应的扩展矢量值，即

其中，j满足如下约束关系：

j＝arg maxΨ(j)≤Ψ_d'(i)。

对σ_d进行按列重排，变换为L×L维的矩阵，即得到直角坐标系下的扩展矢量，仍记为σ_d；

步骤6.4、计算直角坐标系下的最优核函数Φ_opt，即

Φ_opt可以记为Φ_opt(m,n),m＝0,1,…,L-1,n＝0,1,…,L-1。

步骤7、最优核函数Φ_d(m,n)与信号的模糊函数A_z(m,n)相乘，得到模糊域滤波后的特征函数A_F(m,n)，并将A_F(m,n)变换到时频域得到信号的时频分析结果ρ(a,b)：

模糊域滤波后的特征函数A_F(m,n)为

A_F(m,n)＝A_z(m,n)·Φ_opt(m,n)，m＝0,1,…,L-1，n＝0,1,…,L-1。

信号的时频分析结果ρ(a,b)为

其中，a代表离散频率索引，b代表离散频率索引。

本发明第一实施例：

仿真信号为单分量线性调频信号，信号参数设置为：信号幅度A＝1，初始相位

脉宽τ＝0.4s，信号起始频率f₁＝300Hz，信号终止频率f₂＝600Hz，采样频率f_s＝2000Hz，观测数据序列点数L＝τf_s＝800。

依据步骤2，计算直角坐标系下信号采样数据序列s(a)的模糊函数A_z(m,n)；

在步骤3，在模糊域内检测信号自项并估计信号自项的径向角度，结果为1.2566，是弧度制的角度表示；

依据步骤4、步骤5和步骤6，计算出直角坐标系下的模糊域核函数如图2所示；

依据步骤7，将最优核函数Φ_d(m,n)与信号的模糊函数A_z(m,n)相乘，得到模糊域滤波后的特征函数A_F(m,n)，并将A_F(m,n)变换到时频域得到信号的时频分析结果ρ(a,b)，如图3所示。

本发明第二实施例：

仿真信号为线性调频和双曲调频组合信号，信号参数设置为：信号幅度均为A＝1，初始相位均为

双曲调频分量脉宽τ₁＝0.4s，起始时间和终止时间分别为0s和0.4s，起始频率f₁＝300Hz，信号终止频率f₂＝600Hz，采样频率f_s＝2000Hz，观测数据序列点数L₁＝τ₁f_s＝800；线性调频分量脉宽τ₂＝0.25s，起始时间和终止时间分别为0.05s和0.3s，信号起始频率f₃＝125Hz，信号终止频率f₄＝250Hz，观测数据序列点数L₂＝τ₂f_s＝500。

依据步骤2，计算直角坐标系下信号采样数据序列s(a)的模糊函数A_z(m,n)；

在步骤3，在模糊域内检测信号自项并估计信号自项的径向角度，结果为1.3265，是弧度制的角度表示；

依据步骤4、步骤5和步骤6，计算出直角坐标系下的模糊域核函数如图4所示；

依据步骤7，将最优核函数Φ_d(m,n)与信号的模糊函数A_z(m,n)相乘，得到模糊域滤波后的特征函数A_F(m,n)，并将A_F(m,n)变换到时频域得到信号的时频分析结果ρ(a,b)，如图5所示。

本发明第三实施例：

仿真信号为两个平行分量和一个非平行分量组合线性调频信号，信号参数设置为：信号幅度均为A＝1，初始相位均为

分量一脉宽τ₁＝0.6s，起始时间和终止时间分别为0s和0.6s，起始频率f₁＝100Hz，信号终止频率f₂＝62.5Hz，采样频率f_s＝2000Hz，观测数据序列点数L₁＝τ₁f_s＝1200；分量二脉宽τ₂＝0.8s，起始时间和终止时间分别为0s和0.8s，信号起始频率f₃＝200Hz，信号终止频率f₄＝400Hz，观测数据序列点数L₂＝τ₂f_s＝1600；分量三脉宽τ₃＝0.6s，起始时间和终止时间分别为0.2s和0.8s，信号起始频率f₃＝350Hz，信号终止频率f₄＝500Hz，观测数据序列点数L₂＝τ₃f_s＝1200。

依据步骤2，计算直角坐标系下信号采样数据序列s(a)的模糊函数A_z(m,n)；

在步骤3，在模糊域内检测信号自项并估计信号自项的径向角度，结果为[1.6057,1.3963]，是弧度制的角度表示；

依据步骤4、步骤5和步骤6，计算出直角坐标系下的模糊域核函数如图6所示；

依据步骤7，将最优核函数Φ_d(m,n)与信号的模糊函数A_z(m,n)相乘，得到模糊域滤波后的特征函数A_F(m,n)，并将A_F(m,n)变换到时频域得到信号的时频分析结果ρ(a,b)，如图7所示。

可以理解，本发明是通过一些实施例进行描述的，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明所保护的范围内。

Claims (10)

1.一种基于参数预估计的改进径向高斯核时频分析方法，其特征在于，包括如下步骤：

步骤1、获取待处理脉冲信号的采样数据序列s(a)；

步骤2、计算待处理脉冲信号采样数据序列s(a)直角坐标系下的模糊函数A_z(m,n)；

步骤3、在模糊域内检测信号自项并估计信号自项的径向角度φ_k；

步骤4、对待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)进行插值，得到待处理脉冲信号极坐标系下的模糊函数A_p(p,q)；

步骤5、利用Step-Project算法计算待处理脉冲信号的采样数据序列在极坐标系下对应模糊域上半平面的最优径向高斯核函数扩展矢量σ_opt；

步骤6、对极坐标系下的扩展矢量σ进行插值，转换为直角坐标系下的扩展矢量，并计算直角坐标系下的核函数Φ_d(m,n)；

步骤7、将最优核函数Φ_d(m,n)与信号的模糊函数A_z(m,n)相乘，得到模糊域滤波后的特征函数A_F(m,n)，并将A_F(m,n)变换到时频域得到信号的时频分析结果ρ(a,b)。

2.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤1中，采用如下方法获取处理的脉冲信号采样数据序列：

从传感器接收L个采样点的实时采集数据作为待处理脉冲信号的采样数据序列s(a)，已知采样频率f_s，所述的L为包含整个多脉冲信号脉宽长度所对应的采样点个数，取值为2的正整数次幂，且L≥4。

3.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤1中，采用如下方法获取处理的脉冲信号采样数据序列：从存储器中提取包含整个多脉冲信号的L个采样点数据作为待处理的数据序列s(a)，已知采样频率f_s，所述的L为包含整个多脉冲信号脉宽长度所对应的采样点个数，取值为2的正整数次幂，且L≥4。

4.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤2中，计算信号待处理脉冲信号的采样数据序列s(a)直角坐标系下的模糊函数A_z(m,n)的方法，包括以下步骤：

步骤2.1、计算待处理脉冲信号采样数据序列s(a)的离散傅里叶变换S(α)，即：

其中，α为S(α)的离散频率索引，j表示虚数单位，即

利用待处理脉冲信号的离散傅里叶变换S(α)生成s(a)的离散解析信号z(a)的傅里叶变换Z(β)，β＝0,1,…L-1，即：

其中，β为Z(β)的离散频率索引；

对Z(β)进行离散逆傅里叶变换得到待处理脉冲信号的离散解析信号z(a)，即：

步骤2.2、计算离散解析信号z(a)的自相关函数K_z(a,n)，即：

其中，

表示

的共轭形式，n代表离散延时量，n＝0,1,…L-1，a＝0,1,…L-1；

步骤2.3、依据离散解析信号z(a)的自相关函数K_z(a,n)计算待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)，即：

其中，m代表离散频移，n代表离散延时量。

5.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤3中，在模糊域检测信号自项并估计信号自项的径向角度具体包括以下步骤：

步骤3.1、计算待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)在径向距离为0处的离散Radon变换R₀(θ_x)，即：

其中，R₀(θ_x)代表径向距离为0处的Radon变换，θ_x代表Radon变换的离散径向角度，其取值为

x为整数，且0≤x≤X-1，X为总的离散径向角度，X取值为≥2的正整数；δ(·)代表冲激函数；

步骤3.2、估计信号自项在模糊域的径向角度φ_k：

对R₀(θ_x)进行去趋势项处理，得到其去趋势项分量R_d,0(θ_x)，即：

R_d,0(θ_x)＝R₀(θ_x)-R_t,0(θ_x)，

其中，R_t,0(θ_x)代表R₀(θ_x)的趋势项分量，其取值为

其中，B为观测矩阵，为单位矩阵，即：

为采用正则化最小二乘法估计得到的参数，即：

其中，λ为正则化参数，且为大于0的实数，正则化参数λ取值为10，D₂为二阶导数算子的离散形式，D₂是一个(L-2)×(L-1)的二维矩阵，即：

步骤3.3、在得到去趋势项分量R_d,0(θ_x)后，从去趋势项分量R_d,0(θ_x)中搜索所有峰值点(φ_e,R(φ_e))，e＝1,2,...,E，其中φ_e为峰值点所对应的径向角，E为总的峰值点数，为正整数，R(φ_e)表示径向角φ_e所对应的去趋势项分量R_d,0(φ_e)的值；

定义自项的检测门限值R_th为

其中，μ是所有搜索峰值的平均值，即

从R_d,0(θ_x)中搜索出所有大于R_th的径向角，记为φ_k，k＝1,2,...,K，其中K为搜索出的大于R_th的径向角的总个数，即检测到的自项分量的总个数，K≤E。

6.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤4中，采用如下方法对待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)进行插值，转换待处理脉冲信号为极坐标系下的模糊函数A_p(p,q)：

步骤4.1、确定直角坐标系的两个坐标轴对应矢量的取值，横轴对应的矢量为v，纵轴对应的矢量为τ，它们均为长度为L的矢量，分别为

其中，上标中的T代表矩阵转置符号，矩阵中的U为与离散间隔相关的值，

步骤4.2，生成直角坐标系对应的坐标区域，横坐标区域为V，纵坐标区域为Τ，它们均为L×L维矩阵，分别为

Τ_L×L＝[τ…τ]；

步骤4.3，确定极坐标系的两个坐标轴对应矢量的取值，极角对应的矢量为ψ，极径对应的矢量为ρ；

极角矢量为

极径矢量为

其中，ψ的采样间隔为

ρ的采样间隔为

Q为ψ的采样点数，Q＝L，P为ρ的采样点数，

round(·)为取整函数；

步骤4.4、生成极坐标系对应的坐标区域，横坐标区域为X，纵坐标区域为Y，它们均为P×Q维矩阵，分别为

X＝ρ^T×cos(ψ)，

Y＝ρ^T×sin(ψ)；

步骤4.5、运用双线性插值运算将待处理脉冲信号直角坐标系下的模糊函数A_z(m,n)转换为极坐标系下的模糊函数A_p(p,q)，即：

其中，A_p(p,q)为待插值点，{A_z(m,n),A_z(m+1,n),A_z(m,n+1),A_z(m+1,n+1)}为A_z(m,n)上最小包围点A_p(p,q)的相邻的四个模糊函数点值，d_v,min(p,q)为点A_p(p,q)到点A_z(m,n)和点A_z(m,n+1)所在直线的距离，d_h,min(p,q)为点A_p(p,q)到点A_z(m,n)和点A_z(m+1,n)所在直线的距离。

7.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，步骤5中，利用Step-Project算法计算待处理脉冲信号直角坐标系下对应模糊域上半平面的最优径向高斯核函数的扩展矢量σ_opt，具体包括以下步骤：

步骤5.1、确定扩展矢量σ_q的搜索维数R，信号自项径向角对应的极角矢量的索引为

其中，

表示使得括号内取最小值的q值；

扩展矢量σ_q中，对于每一个自项径向角对应的极角矢量索引，定义极角矢量索引的单边搜索索引宽度为N_s，

其中round()代表四舍五入取整数运算，则需要搜索的索引范围为

S_k＝{l_k-N_s,l_k-N_s+1,…,l_k,…l_k+N_s-1,l_k+N_s}；

扩展矢量σ_q中需要搜索的总索引范围为

记S的元素个数为R，即扩展矢量σ_q的搜索维数；

步骤5.2、确定Step-Project算法中要使用的各元素的初值：

μ(0)＝1，μ_m＝1001，

ε₁＝0.1，ε₂＝1-ε₁，

Δ＝10，ε₃＝1×10^-5，

κ＝0；

其中，γ为核函数的体积限制参数，σ(0)为扩展矢量的初始值，μ(0)为初始步长，μ_m为最大步长，ε₁和ε₂为测试步长大小的标准值，ε₃为迭代终止误差中设置的参数，Δ为步长倍增或倍减因子，κ为迭代次数；

步骤5.3、开始进行Step-Project算法，子步骤如下：

步骤5.3.1、求梯度向量▽f(κ)，计算目标函数f(σ_q(κ))的值，即

其中，σ_q(κ)为第κ次迭代时的扩展矢量的值，▽代表求梯度，

代表f对σ_q(κ)求偏导数；

步骤5.3.2、更新下一次迭代的扩展矢量和优化目标，即

步骤5.3.3、根据f(σ_q(κ))与f(σ_q(κ-1))测试步长μ(κ)，若过小则令步长倍增，若超过最大步长μ_m则停止迭代，转步骤5.4，若在μ_m以内过大则令步长倍减；定义判断测试步长大小的依据g为优化目标实际上升值与预期上升值的比值，即

若g＜ε₁，则令μ(κ+1)＝μ(κ)/Δ，σ_q(κ+1)＝σ_q(κ)；若g＞ε₂，则令

若ε₁≤g≤ε₂，则设置两个迭代终止误差为

且迭代次数递增，即令κ＝κ+1；

步骤5.3.4、判断若迭代终止误差值均为非正，或迭代次数κ超过最大迭代次数50，则停止迭代，输出σ_q(κ)，记为σ_q'，否则转步骤5.3.1；

步骤5.4：确定完整的最优径向高斯核函数的扩展矢量σ_opt，首先预置

再代入有值部分，即

σ_opt(S(i))＝σ_q'(i),i＝1,2,…,R。

8.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤6中，对极坐标系下的扩展矢量插值转换为直角坐标下的扩展矢量，具体包括以下步骤：

步骤6.1：将ψ和σ_opt补充为整个平面上的表示，即

ψ＝[ψ,ψ(i)+π],i＝2,3,…,Q，

σ＝[σ_opt,σ_opt(i)],i＝2,3,…,Q；

步骤6.2：依据直角坐标轴V和Τ计算它们对应于极坐标中的角度，即

且需满足0＜Ψ_d＜2π；

步骤6.3：首先将直角坐标点的角度转换为一个列向量，即Ψ_d'，再对极坐标系中的扩展矢量进行线性插值得到每个直角坐标点对应的扩展矢量值，即

其中，j满足如下约束关系：

j＝arg maxΨ(j)≤Ψ_d'(i)；

对σ_d进行按列重排，变换为L×L维的矩阵，即得到直角坐标系下的扩展矢量，仍记为σ_d；

步骤6.4：计算直角坐标系下的最优核函数Φ_opt，即

Φ_opt记为Φ_opt(m,n),m＝0,1,…,L-1,n＝0,1,…,L-1。

9.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤7中，模糊域滤波后的特征函数A_F(m,n)为

A_F(m,n)＝A_z(m,n)·Φ_opt(m,n)，m＝0,1,…,L-1，n＝0,1,…,L-1；

信号的时频分析结果ρ(a,b)为

其中，a代表离散频率索引，b代表离散频率索引。

10.根据权利要求1所述的基于参数预估计的改进径向高斯核时频分析方法，其特征在于，所述步骤7中，模糊域到时频域的转换用快速傅里叶变换实现。

Images