一种多分量多项式相位信号的分离方法

摘要

本发明公开一种多分量多项式相位信号的分离方法，步骤是：步骤A，准备待检测处理的分量个数为NUM的多分量多项式相位信号，该信号以采样间隔1/f

说明书

技术领域

本发明涉及一种多项式相位信号模型，特别涉及一种基于时频分布的多分量多项式相位信号分离方法。

背景技术

多分量多项式相位信号的分离方法能够很好地实现在时-频平面上存在交点信号之间的分离，也能实现只在一段时间内存在的相位信号的分离。

多项式相位信号具有非平稳的特性，在自然界或无线通信、雷达、地震、生物医学等领域有广泛的应用。多项式相位信号可分为单分量、多分量、高信噪比、低信噪比等情况。单分量和高信噪比前提下的多项式相位信号分析具有较完善的理论，但是这样的传统分析理论不适用于多分量、低信噪比和非高斯噪声环境下的相位信号分离。因此这方面的研究具有现实意义。

B.Barkat和K.Abed-Meraim首先在《Algorithms for Blind ComponentsSeparation and Extraction from the Time-Frequency Distribution of TheirMixture》中估计了分量的个数、进行分量分离，然后进行脊提取和脊路径重组，并使用一个给定的混合信号的无交叉项TFD来分离所有的分量。李宏坤在《基于时频分布的欠定信号盲分离与微弱特征提取》中基于时频分析的盲源分离方法，结合时频分析对非平稳信号分析的优势进行分离。向强等在《基于线性正则变换的时频信号分离方法》中把在时频面互不重叠，但在时域和频域均存在较强耦合的多分量合成信号进行分离。上述方法中，当信号的相位多项式阶数较高或信号分量数目较多时，会存在过多的交叉项干扰，因此具有一定的局限性。

发明内容

本发明的目的，在于提供一种多分量多项式相位信号的分离方法，能够克服传统的分离方法不能在多分量下有效分离多项式相位信号的问题。

为了达成上述目的，本发明的解决方案是：

一种多分量多项式相位信号的分离方法，包括如下步骤：

步骤A，准备待检测处理的分量个数为NUM的多分量多项式相位信号，该信号以采样间隔1/f_s进行离散化，得到s(n),1≤n≤N，其中，f_s为采样频率，N为数据长度；

步骤B，令ss(n)＝s(n)，number＝1；

步骤C，用时频分析方法计算信号ss(n)的时频分布TFR(t,f)，其中，t为时频分布时间方向上的矢量，长度为N；f为时频分布频率方向上的矢量，长度为M，M为进行时频分析时FFT变换的频率点数；

步骤D，对信号的时频分布通过二维峰值搜索算法进行脊线提取，获得相位多项式信号的瞬时频率向量IF；

步骤E，对步骤D获得的相位多项式信号瞬时频率向量IF进行解调频和滤波；

步骤F，修正平滑后的瞬时频率向量F_smooth(n)，得到最终的输出瞬时频率向量IF_Out(n)；

步骤G，将信号ss(n)减去已提取出的分量Out_ss(n)，得到新的ss(n)；number＝number+1；如果number≤NUM，跳转至步骤C；否则，分离程序结束。

上述步骤D的具体步骤是：

步骤D-1，设置频率搜索范围delta，搜索脊线峰值时的幅度门限值thres，求脊线斜率时两点之间的时间间隔tdis；

步骤D-2，对所得到的时频分布取绝对值，并寻找TFR(t,f)最大值所对应的时间t₀和频率f₀；

步骤D-3，对时频分布进行归一化，得到归一化后的时频分布TFD(t,f)；

步骤D-4，设置时频信号分量的起始时刻为tstart，结束时刻tend，设置出现时频断点时的频率搜索范围初始值delta0，令估计的相位多项式信号分量瞬时频率向量为IF，长度为N，并令IF(t₀)＝f₀；

步骤D-5，以时频分布最大值为起点，先在时间轴方向向右搜索；设时间搜索变量j＝t₀+1，频率搜索变量f_R为f₀；

步骤D-6，如果j>N，跳转至步骤D-8，否则，确定频率搜索范围的上下界，其中low为频率搜索范围的下界，low＝max{f_R-delta,1}，up为频率搜索范围的上界，up＝min{f_R+delta,M}；在上下界范围内以设定的门限thres寻找脊线峰值，提取其在j时刻脊线峰值对应的频率向量F_j(i)，其中i为峰值个数；

步骤D-7，对提取的脊线峰值频率向量F_j(i)进行分析处理；

步骤D-8，以时频分布最大值为起点，再在时间轴方向向左搜索，设时间搜索变量k＝t₀-1，频率搜索变量f_L为f₀；

步骤D-9，如果k<1，跳转至步骤E；否则，确定频率搜索范围的上下界，其中low为频率搜索范围的下界，low＝max{f_L-delta,1}，up为频率搜索范围的上界，up＝min{f_L+delta,M}；在上下界范围内以设定的门限thres寻找脊线峰值，提取其在k时刻脊线峰值对应的频率向量F_k(i)，其中i为峰值个数；

步骤D-10，对提取的脊线峰值频率向量F_k(i)进行分析处理。

上述步骤D-7中，对提取的脊线峰值频率向量F_j(i)进行分析处理的具体过程是：

当没有出现脊线峰值，即i＝0时，将结束时刻修改为tend＝j-1，令jj＝j；

步骤D-7-1，如果jj>N，跳转至步骤D-8，否则，重新确定频率搜索范围的上下界，其中频率下界low＝max{f_R-delta0,1}，频率上界up＝min{f_R+delta0,M}；在此上下界频率范围内求TFD(t,f)最大值，并提取其对应的频率f_R，令IF(j)＝f_R，j＝j+1，jj＝jj+1；

步骤D-7-2，如果TFD(t,f)最大值大于门限值thres，将结束时刻修改为tend＝N，并跳出当前循环至步骤D-6，否则，跳转至步骤D-7-1；

当出现一个脊线峰值，即i＝1时，令IF(j)＝F_j(1)，j＝j+1，跳转至步骤D-6；

当出现多个脊线峰值，即i>1时，如果j-t₀≥2·tdis，则计算频率向量的斜率差dis1(i)：

取dis1(i)最小值对应的频率f_R，并令IF(j)＝f_R，j＝j+1，返回步骤D-6；

如果j-t₀<2·tdis，则计算频率向量的距离dis2(i)：

dis2(i)＝abs[F_j(i)-IF(j-1)]

取dis2(i)最小值对应的频率f_R，并令IF(j)＝f_R，j＝j+1，返回步骤D-6。

上述步骤D-10中，对提取的脊线峰值频率向量F_k(i)进行分析处理进行分析处理的具体过程是：

当没有出现脊线峰值，即i＝0时，将开始时刻修改为tstart＝k+1，kk＝k；

步骤D-10-1，如果kk<1，跳转至步骤E，否则，重新确定频率搜索范围的上下界，其中频率下界low＝max{f_L-delta0,1}，频率上界up＝min{f_L+delta0,M}；在此频率上下界范围内求TFD(t,f)最大值，并提取其对应的频率f_L，令IF(k)＝f_L，k＝k-1，kk＝kk-1；

步骤D-10-2，如果最大值大于门限值thres，将开始时刻修改为tstart＝1，并跳出当前循环至步骤D-9，否则，跳转至步骤D-10-1；

当出现一个脊线峰值，即i＝1时，令IF(k)＝F_k(1)，k＝k-1，跳转至步骤D-9；

当出现多个脊线峰值，即i>1时，如果N-k≥2·tdis，则计算频率向量的斜率差dis3(i)：

取dis3(i)最小值对应的频率f_L，并令IF(k)＝f_L，k＝k-1，跳转至步骤D-9；

如果N-k<2·tdis，则计算频率向量的距离dis4(i)：

dis4(i)＝abs[F_k(i)-IF(k+1)]

取dis4(i)最小值对应的频率f_L，并令IF(k)＝f_L，k＝k-1，跳转至步骤D-9。

上述步骤E的具体步骤是：

步骤E-1，将瞬时频率向量IF转换至标准频率得到新的瞬时频率向量IF_New：

IF_New＝(IF-N/2)/f_s

步骤E-2，对瞬时频率向量IF_New进行平滑滤波得到平滑后的频率向量F_smooth；

步骤E-3，对平滑后的频率向量进行积分得到相位多项式信号瞬时相位向量phase(n)；

步骤E-4，对多分量相位多项式信号进行解调频：

S_comp(n)＝ss(n)·exp(-j·phase(n))(tstart≤n≤tend)

步骤E-5，对解调频后的S_comp(n)信号进行低通滤波获得该分量的瞬时幅度向量IA(n)；

步骤E-6，根据瞬时相位向量phase(n)和瞬时幅度向量IA(n)重构信号，得到该分量的时域形式Out_ss(n)。

上述步骤E-2的具体内容是：用修正二阶差分矩阵Ω，通过下式进行平滑：

其中，

其中，I是单位矩阵，T是转置运算，beta＝10^-4为平滑度，其值越小，则频率曲线越平滑。

上述步骤E-5中，滤波器采用线性相位FIR滤波器，窗函数选择为汉明窗，截止频率为CutFreq＝10Hz，滤波器长度为[1000×0.8]，其中，运算符号[x]表示取不大于x的最大整数运算。

上述步骤F的具体步骤是：

步骤F-1，计算标准频率向量：

步骤F-2，令kk＝1；

步骤F-3，如果kk<1或kk>N跳转至步骤G，否则，dis5(kk,n)＝abs[F_smooth(kk)-F_std(n)]，求出dis5(kk,n)的最小值所对应的F_std(n)中的频率值f_kk；

步骤F-4，令IF_Out(kk)＝f_kk，kk＝kk+1，跳转至步骤F-3。

采用上述方案后，本发明利用对多分量多项式相位信号的时频分布进行脊线提取，通过确定频率搜索的上下界和峰值幅度的门限值来统计峰值个数。在此基础上根据峰值个数的不同通过不同的方法筛选频率，得到瞬时频率向量；而后对瞬时频率向量进行平滑得到平滑后的瞬时频率向量，对其积分得到瞬时相位向量，在此基础上通过解调频的方法重构信号，对重构的信号滤波得到该分量对应的时域形式，将该分量从原信号中减去，重复进行以上步骤直到将所有多项式相位信号分量全部提取出来。

本发明能够克服传统的分离方法不能在多分量下有效分离多项式相位信号的问题，本发明描述了多项式相位信号的时频分布特性，并在此基础上通过二维峰值搜索算法提取信号的时频图像脊线，通过设置频率查找范围和判断峰值的门限值对信号进行分离处理，最终将所有分量提取出来。由于多分量相位信号是非平稳信号处理的重要研究对象，本发明能够很好地实现在时-频平面上存在交点的多分量相位多项式信号之间的分离，也能实现只在一时间段内存在的相位多项式信号的分离。

附图说明

图1是待检测的多分量相位多项式信号的时域波形图s(n)；

图2是待检测信号的时-频图；

图3是提取第一个多项式信号分量后的剩余信号的时-频图；

图4是提取前两个多项式信号分量后的剩余信号的时-频图；

图5是提取前三个多项式信号分量后的剩余信号的时-频图；

图6是提取出的各个相位多项式信号分量的瞬时频率。

具体实施方式

以下以多项式相位信号为例，该信号在时-频平面上存在交点，且部分分量在某段时间内不存在。结合附图，详细说明本发明的实施方式。

本发明是一种多分量多项式相位信号的分离方法，包括如下步骤：

步骤A：准备待检测处理的分量个数为NUM＝4的多分量多项式相位信号，该信号以采样间隔1/f_s进行离散化，得到s(n),1≤n≤N，其中，f_s＝1000Hz为采样频率，N＝1000为数据长度。时域波形见图1。

s₁(n)＝A₁(n)exp(2πj·(200·(n/fs)²))；

s₂(n)＝A₂(n)exp(j·(140π·(n/fs)+150·sin(2π·(n/fs))))；

s₃(n)＝A₃(n)exp(j·(-140π·(n/fs)-150·sin(2π·(n/fs))))；

s(n)＝s₁(n)+s₂(n)+s₃(n)+s₄(n)

式中，信号幅度A₁(n)＝A₂(n)＝A₃(n)＝A₄(n)＝1。

步骤B：令ss(n)＝s(n)，number＝1。

步骤C：用时频分析方法计算信号ss(n)的时频分布TFR(t,f)，这里采用短时傅里叶方法。其中，t为时频分布时间方向上的矢量，长度为N；f为时频分布频率方向上的矢量，长度为M，M＝1024为进行时频分析时FFT变换的频率点数。

TFR(t,f)＝STFT_ss(t,f)

步骤D：对信号的时频分布通过二维峰值搜索算法进行脊线提取，获得相位多项式信号的瞬时频率向量IF。其过程如下：

步骤D-1：设置频率搜索范围delta＝20，搜索脊线峰值时的幅度门限值thres＝0.1，求脊线斜率时两点之间的时间间隔tdis＝35。

步骤D-2：对所得到的时频分布取绝对值，并寻找TFR(t,f)最大值所对应的时间t₀和频率f₀，其算法实现为：

TFR0(t,f)＝abs[TFR(t,f)]

[t₀,f₀]＝argmax_t,f[TRF0(t,f)]

步骤D3：对时频分布进行归一化，得到归一化后的时频分布TFD(t,f)：

TFD(t,f)＝TFR0(t,f)/max[TFR0(t,f)]

步骤D-4：设置时频信号分量的起始时刻为tstart＝1，结束时刻tend＝N，设置出现时频断点时的频率搜索范围初始值delta0＝2，令估计的相位多项式信号分量瞬时频率向量为IF，长度为N，并令IF(t₀)＝f₀。

步骤D-5：以时频分布最大值为起点，先在时间轴方向向右搜索。设时间搜索变量j＝t₀+1，频率搜索变量f_R为f₀。

步骤D-6：如果j>N，跳转至步骤D-8，否则，确定频率搜索范围的上下界，其中low为频率搜索范围的下界，low＝max{f_R-delta,1}，up为频率搜索范围的上界，up＝min{f_R+delta,M}。在上下界范围内以设定的门限thres寻找脊线峰值，提取其在j时刻脊线峰值对应的频率向量F_j(i)，其中i为峰值个数。

步骤D-7：对提取的脊线峰值频率向量F_j(i)进行分析处理。

当没有出现脊线峰值，即i＝0时，将结束时刻修改为tend＝j-1，令jj＝j；

步骤D-7-1：如果jj>N，跳转至步骤D-8，否则，重新确定频率搜索范围的上下界，其中频率下界low＝max{f_R-delta0,1}，频率上界up＝min{f_R+delta0,M}。在此上下界频率范围内求TFD(t,f)最大值，并提取其对应的频率f_R，令IF(j)＝f_R，j＝j+1，jj＝jj+1。

步骤D-7-2：如果TFD(t,f)最大值大于门限值thres，将结束时刻修改为tend＝N，并跳出当前循环至步骤D-6，否则，跳转至步骤D-7-1。

当出现一个脊线峰值，即i＝1时，令IF(j)＝F_j(1)，j＝j+1，跳转至步骤D-6。

当出现多个脊线峰值，即i>1时，如果j-t₀≥2·tdis，则计算频率向量的斜率差dis1(i)：

取dis1(i)最小值对应的频率f_R，并令IF(j)＝f_R，j＝j+1，返回步骤D-6。

如果j-t₀<2·tdis，则计算频率向量的距离dis2(i)：

dis2(i)＝abs[F_j(i)-IF(j-1)]

取dis2(i)最小值对应的频率f_R，并令IF(j)＝f_R，j＝j+1，返回步骤D-6。

步骤D-8：以时频分布最大值为起点，再在时间轴方向向左搜索。设时间搜索变量k＝t₀-1，频率搜索变量f_L为f₀。

步骤D-9：如果k<1，跳转至步骤E；否则，确定频率搜索范围的上下界，其中low为频率搜索范围的下界，low＝max{f_L-delta,1}，up为频率搜索范围的上界，up＝min{f_L+delta,M}。在上下界范围内以设定的门限thres寻找脊线峰值，提取其在k时刻脊线峰值对应的频率向量F_k(i)，其中i为峰值个数。

步骤D-10：对提取的脊线峰值频率向量F_k(i)进行分析处理。

当没有出现脊线峰值，即i＝0时，将开始时刻修改为tstart＝k+1，kk＝k；

步骤D-10-1：如果kk<1，跳转至步骤E，否则，重新确定频率搜索范围的上下界，其中频率下界low＝max{f_L-delta0,1}，频率上界up＝min{f_L+delta0,M}。在此频率上下界范围内求TFD(t,f)最大值，并提取其对应的频率f_L，令IF(k)＝f_L，k＝k-1，kk＝kk-1。

步骤D-10-2：如果最大值大于门限值thres，将开始时刻修改为tstart＝1，并跳出当前循环至步骤D-9，否则，跳转至步骤D-10-1。

当出现一个脊线峰值，即i＝1时，令IF(k)＝F_k(1)，k＝k-1，跳转至步骤D-9。

当出现多个脊线峰值，即i>1时，如果N-k≥2·tdis，则计算频率向量的斜率差dis3(i)：

取dis3(i)最小值对应的频率f_L，并令IF(k)＝f_L，k＝k-1，跳转至步骤D-9；

如果N-k<2·tdis，则计算频率向量的距离dis4(i)：

dis4(i)＝abs[F_k(i)-IF(k+1)]

取dis4(i)最小值对应的频率f_L，并令IF(k)＝f_L，k＝k-1，跳转至步骤D-9。

步骤E：对步骤D获得的相位多项式信号瞬时频率向量IF进行解调频和滤波，其过程如下：

步骤E-1：将瞬时频率向量IF转换至标准频率得到新的瞬时频率向量IF_New：

IF_New＝(IF-N/2)/f_s

步骤E-2：对瞬时频率向量IFNew进行平滑滤波得到平滑后的频率向量F_smooth。平滑的具体过程为：用修正二阶差分矩阵Ω，通过下式进行平滑：

这里，

其中，I是单位矩阵，T是转置运算，beta＝10^-4为平滑度，其值越小，则频率曲线越平滑。

步骤E-3：对平滑后的频率向量进行积分得到相位多项式信号瞬时相位向量phase(n)。

步骤E-4：对多分量相位多项式信号进行解调频：

S_comp(n)＝ss(n)·exp(-j·phase(n))(tstart≤n≤tend)

步骤E-5：对解调频后的S_comp(n)信号进行低通滤波获得该分量的瞬时幅度向量IA(n)。滤波器采用线性相位FIR滤波器，窗函数选择为汉明窗，截止频率为CutFreq＝10Hz，滤波器长度为[1000×0.8]，其中，运算符号[x]表示取不大于x的最大整数运算。

步骤E-6：根据瞬时相位向量phase(n)和瞬时幅度向量IA(n)重构信号，得到该分量的时域形式Out_ss(n)。

Out_ss(n)＝IA(n)·exp(j·phase(n))(tstart≤n≤tend)

步骤F：修正平滑后的瞬时频率向量F_smooth(n)，得到最终的输出瞬时频率向量IF_Out(n)。

步骤F-1：计算标准频率向量：

步骤F-2：kk＝1。

步骤F-3：如果kk<1或kk>N跳转至步骤G，否则，dis5(kk,n)＝abs[F_smooth(kk)-F_std(n)-，求出dis5(kk,n)的最小值所对应的F_std(n)中的频率值f_kk。

步骤F-4：令IF_Out(kk)＝f_kk，kk＝kk+1，跳转至步骤F-3。

步骤G：将信号ss(n)减去已提取出的分量Out_ss(n)，得到新的新的ss(n)；number＝number+1；如果number≤NUM，跳转至步骤C；否则，分离程序结束。

图3至图5为多项式相位信号分量逐个提取出来后的时-频分布图，图6为所有分量的瞬时频率曲线表示图。从图中可以看出，分量被逐个从时-频分布图中提取出来，其中，第一个分量虽然在前一段时间内不存在，但仍旧被完整的提取出来；前一个分量提取出来后，会在它与其它分量的交点处出现断层，但从最终结果来看，剩余分量仍能被完整地提取出来。因此，可以得出结论，本方法能够很好地实现在时-频平面上对多分量多项式相位信号的分离。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。