一种多分量线性调频信号时频分析方法

文档序号:31053997发布日期:2022-08-06 10:23阅读:254来源:国知局
一种多分量线性调频信号时频分析方法

1.本发明涉及多分量非平稳信号时频估计技术领域,具体涉及一种多分量线性调频信号时频分析方法。


背景技术:

2.线性调频信号是一种非平稳信号,线性调频信号在雷达探测、转轴故障分析、水下通信系统中用途较为广泛。在水下通信系统中利用多个线性调频信号进行信道探测,由于水下信道的复杂性以及各种噪声的干扰,最终在接收端会收到一个受多普勒效应影响且带有噪声以及时延的多分量线性调频信号,所以对多分量线性调频信号进行时频分析研究具有重要意义。
3.目前已有的时频分析方法主要分为线性时频分析和双线性时频分析,传统的线性时频分析方法有短时傅里叶变换(stft)、小波变换(wt),双线性时频分析方法有wigner-ville分布(wvd)。目前一些基于stft改进的同步压缩算法和基于wvd改进的二次时频算法,处理时频域没有交叉的多分量线性调频信号可以得到能量聚集性较高、时频分辨率较高的时频分布,但对时频域存在交叉的多分量线性调频信号处理结果较差,在时频域交叉点处会发生时频脊线失真。针对此问题,可以用模式分解算法进行时频分析,比如基于自适应线性调频模式分解(acmd)算法结合信号残差获得自适应时频谱,acmd算法能较好的分解多分量线性调频信号,但处理时频域存在交叉的多分量线性调频信号时,残差的时频分布会出现时频脊线断裂问题,导致后续信号分量的时频脊线估计误差较大进而影响最终的自适应时频谱。


技术实现要素:

4.本发明的目的是为了解决现有技术中多分量线性调频信号时频分析方法存在的问题,提供一种多分量线性调频信号时频分析方法
5.本发明的目的可以通过采取如下技术方案达到:
6.一种多分量线性调频信号时频分析方法,所述时频分析方法包括以下步骤:
7.s1、对水声通信系统接收到的多分量线性调频信号以ts为采样周期采样得到离散的多分量线性调频信号s(n
t
),并定义非完全残差r(n
t
),利用非完全残差解决时频脊线断裂情况;
8.s2、对非完全残差r(n
t
)做短时傅里叶变换得到时频矩阵tf(n
t
,nf);
9.s3、采用改进的时频脊线估计算法估计tf(n
t
,nf)中具有最大能量的时频脊线向量i表示第i个分量,1≤i≤m,m是信号分量总个数;
10.s4、将步骤s3中得到的作为acmd算法的初始频率,利用acmd算法分解非完全残差r(n
t
)得到对应的信号分量瞬时幅度和更新的对更新的做线性拟合得到对应的调频斜率估计值起始频率估计值以及时频脊线向量最终估计值更新
非完全残差r(n
t
),并判断是否已存在,根据是否存在的判断结果确定是否需要更新阈值残差r
th
(n
t
);
11.s5、判断阈值残差r
th
(n
t
)是否达到设定的结束阈值(1-α)s(n
t
),如果阈值残差大于结束阈值则跳转到s2继续执行,否则整个时频分析方法结束得到各个分量信号的时频脊线向量瞬时幅度调频斜率和起始频率
12.进一步地,所述步骤s1过程如下:
13.对水声通信系统接收到的多分量线性调频信号以ts为采样周期采样得到离散的多分量线性调频信号s(n
t
),并定义非完全残差表达式,离散多分量线性调频信号s(n
t
)的表达式如下:
[0014][0015]
其中n
t
是s(n
t
)的长度,m是s(n
t
)的分量总个数,n(n
t
)为高斯白噪声,ai(n
t
)是第i个分量信号的幅度,fi是第i个分量信号的起始频率,ki是第i个分量信号的调频斜率;
[0016]
在基于acmd算法结合残差得到的自适应时频谱方法中,残差定义为接收信号与当前估计到的分量信号之间的差,记作:
[0017][0018]
其中表示分量信号i的估计值,为了处理基于acmd算法结合残差的自适应时频谱存在的时频脊线断裂问题,在计算残差时保留当前估计到的分量信号i的部分能量,定义非完全残差r(n
t
),表达式如下:
[0019][0020]
其中,0<α<1,α的取值由信噪比snr决定,即α=g(snr),噪声能量越小即信噪比snr越大,残差的时频脊线断裂越严重,所以在信噪比大时需要保留分量信号i更多的能量以解决残差时频脊线断裂问题,因此定义g(snr)在(-∞,+∞)上是单调递减函数,由一元线性函数、双边指数函数或argtan函数构成,非完全残差r(n
t
)的初始值是接收信号s(n
t
),即r(n
t
)=s(n
t
)。
[0021]
进一步地,所述步骤s2过程如下:
[0022]
对非完全残差r(n
t
)做短时傅里叶变换得到时频矩阵tf(n
t
,nf),表达式如下:
[0023][0024]
其中w(m)是窗函数,w(n
t-m)表示窗函数在时间上反转并有n
t
个样本偏移量,x(m)是短时傅里叶变换的输入信号,m是输入信号x(m)的采样点,取残差r(n
t
)作为输入信号,nf表示以tf为采样周期对频率采样得到的频率总点数,nf是取值从1到nf的频率采样点。
[0025]
进一步地,所述步骤s3过程如下:
[0026]
根据基于acmd算法结合残差的时频分析方法,利用与贪心算法类似的原理每次从非完全残差中分解出能带走最多能量的信号分量,执行acmd算法时需要提供初始频率值,所以需要估计非完全残差时频分布中具有最大能量的时频脊线,具体估计步骤如下:
[0027]
s301、在tf(n
t
,nf)中找到具有最大能量的时间频率点(m
t
,mf),表达式如下:
[0028][0029]
其中m
t
、mf分别是tf(n
t
,nf)中最大能量点对应的时间点和频率点,记录具有最大能量的时间频率点其中1≤m
t
≤n
t
,1≤mf≤nf;
[0030]
s302、以时间频率点(m
t
,mf)为起始点估计时间点m
t
右侧从时间点m
t
+1到时间点n
t
对应的频率点,初始化估计过程用下式表示:
[0031][0032]
其中δnf为设定的频率点增值,r
t
是时间点m
t
右侧的时间点,是时间点r
t
对应的频率点,r
t
从m
t
开始,每执行一次式(6),r
t
就需要加1,r
t
更新公式为r
t
=r
t
+1,直到r
t
=n
t
时结束估计过程,得到时间点m
t
右侧时频脊线向量对用线性拟合方法得到m
t
右侧时频脊线的斜率lk;
[0033]
为了使得m
t
左右侧估计的时频脊线属于同一分量信号,时间点m
t
左侧时频脊线的斜率应该和m
t
右侧估计的时频脊线斜率一致或者接近,所以下面按照m
t
右侧时频脊线斜率lk去估计时间点m
t
左侧的时频脊线;
[0034]
s303、以时间频率点(m
t
,mf)为起始点,按照斜率lk估计时间点m
t
左侧从时间点m
t-1到时间点1对应的频率点,初始化估计过程用下式表示:
[0035][0036]
其中δnf为设定的频率点增值表示向下取整,表示向上取整,l
t
是时间点m
t
左侧的时间点,是时间点l
t
对应的频率点,l
t
从m
t
开始,每执行一次式(7),l
t
减1,l
t
更新公式为l
t
=l
t-1,直到l
t
=1时结束估计过程,得到完整的时频脊线向量
[0037]
s304、对估计到的时频脊线向量用基于最小二乘法则的拟合方法进行一阶线性拟合更新。
[0038]
进一步地,所述步骤s4过程如下:
[0039]
s401、以步骤s3中得到的作为acmd算法的初始频率,对非完全残差r(n
t
)使用acmd算法分解得到分量信号对应的瞬时幅度更新的时频脊线向量并对时频脊线向量用基于最小二乘法则的拟合方法进行一阶线性拟合得到当前估计的分量信号i的时频脊线向量和调频斜率起始频率
[0040]
在使用非完全残差时,保留已估计分量信号的部分能量,但每个分量信号的能量大小不一致,如果非完全残差中已估计分量信号保留的部分能量要大于其他未被估计的分量信号能量,则在执行步骤s3时会得到重复的时频脊线,导致步骤s401中acmd分解非完全残差得到的信号分量重复,所以需要对s401得到的结果去重处理;
[0041]
s402、当前估计的时频脊线向量为判断是否已经存在:
[0042][0043]
其中j为正整数表示已经得到的估计分量信号,1≤j≤i-1,δif为自定义的频率增量向量,根据下式更新非完全残差r(n
t
),表达式如下:
[0044][0045]
如果式(8)成立则表示步骤s401得到的时频脊线向量重复直接跳转到步骤s2,如果式(8)不成立则按照以下式(10)更新阈值残差r
th
(n
t
):
[0046][0047]
阈值残差r
th
(n
t
)只有在当前估计到的分量信号i与已有的估计分量信号j不重复时,才执行更新。
[0048]
进一步地,所述步骤s5过程如下:
[0049]
无论时频脊线是否重复估计都需要更新非完全残差,所以不能根据非完全残差值判定算法是否结束,步骤s402定义了新的阈值残差用来判定算法是否结束;根据式(11)判断阈值残差是否达到设定的结束阈值(1-α)s(n
t
):
[0050]rth
(n
t
)≤(1-α)s(n
t
)
ꢀꢀ
(11)
[0051]
如果式(11)成立,得到基于非完全残差的自适应时频谱:
[0052][0053]
其中round()表示四舍五入取整数的函数,δ()表示dirac函数,是分量i的幅度估计值,是分量i的时频脊线估计值,如果式(11)不成立,则跳转到步骤s2继续执行。
[0054]
进一步地,所述步骤s5中当式(11)成立,利用已经得到的分量信号i对应的瞬时幅度估计值和时频脊线估计值以及理想时频分布原理得到基于非完全残差的自适应时频谱的同时,得到每个分量信号对应的调频斜率估计值和起始频率估计值
[0055]
本发明相对于现有技术具有如下的优点及效果:
[0056]
1、该方法适应性广。本发明不需要知道多分量线性调频信号的先验信息信号分量总个数,通过阈值残差和设定的结束阈值大小关系判定整个时频分析算法是否结束,同时对估计的时频脊线进行去重处理,避免得到重复的时频脊线。
[0057]
2、该方法实现步骤简单。本发明通过定义的非完全残差解决基于acmd算法结合残差的时频分析方法处理时频域有交叉的多分量线性调频信号时存在的时频脊线断裂问题,提高时频脊线估计以及时频谱准确度。
[0058]
3、该方法能获得每个线性调频信号分量较为精确的起始频率和调频斜率估计值。本发明利用改进的时频脊线估计算法避免时频分布中最大能量点右侧时频脊线向量和左侧时频脊线向量不属于同一信号分量,用基于最小二乘法则的线性拟合方法对时频脊线向量进行一阶线性拟合得到对应的起始频率和调频斜率。
附图说明
[0059]
此处所说明的附图用来提供对本发明的进一步理解,构成本技术的一部分,本发明的示意性实施例及其说明用于解释本发明,并不构成对本发明的不当限定。在附图中:
[0060]
图1是本发明实施1中信噪比10db下完全去除分量1之后的残差的stft变换结果;
[0061]
图2是本发明实施例1中信噪比10db下不完全去除分量1之后的非完全残差stft变换结果;
[0062]
图3是本发明实施例1中信噪比是10db下的第一个时频脊线估计结果图;
[0063]
图4是本发明实施例1中信噪比是10db下的基于非完全残差得到的自适应时频图;
[0064]
图5是本发明实施例2信噪比是5db下的多分量线性调频信号stft变换结果仿真图;
[0065]
图6是本发明实施例2信噪比是5db下的多分量线性调频信号基于非完全残差得到的自适应时频图;
[0066]
图7是本发明实施例3信噪比是5db下的多分量线性调频信号stft变换结果仿真图;
[0067]
图8是本发明实施例3信噪比是5db下的多分量线性调频信号基于非完全残差得到的自适应时频图;
[0068]
图9是本发明中公开的一种多分量线性调频信号时频分析方法的流程图。
具体实施方式
[0069]
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
[0070]
下面结合实施例对本发明作进一步的详细描述,但本发明的实施方式不限于此。
[0071]
实施例1
[0072]
如图9所示,本实施例公开了一种多分量线性调频信号时频分析方法,包含以下步骤:
[0073]
s1、对水声通信系统接收到的多分量线性调频信号以ts=0.002s为采样周期采样得到离散多分量线性调频信号s(n
t
),并定义非完全残差表达式,s(n
t
)表达式如下:
[0074][0075]
其中n
t
=3001,m=3,i=1,2,3,n(n
t
)是均值为0,方差为0.1189的高斯白噪声,信噪比snr=10,a1(n
t
)=exp(-0.03
·nt
·
t
t
),f1=25.098hz,k1=15.123hz/s,a2(n
t
)=exp(-0.03
·nt
·
t
t
),f2=40.8921hz,k2=8.23hz/s,a3(n
t
)=exp(-0.09
·nt
·
t
t
),f3=
60.023hz,k3=12.678hz/s。
[0076]
定义非完全残差,其表达式如下:
[0077]
其中表示分量信号i的估计值,1≤i≤m,0<α<1,α的取值由信噪比snr决定,即α=g(snr),g(snr)在(-∞,+∞)上是单调递减函数,由一元线性函数、双边指数函数或argtan函数构成,非完全残差r(n
t
)的初始值是s(n
t
),即r(n
t
)=s(n
t
)。本实施例中α定义为:
[0078][0079]
s2、对非完全残差r(n
t
)做短时傅里叶变换得到时频矩阵tf(n
t
,nf),表达式如下:
[0080]
其中w(m)是窗函数,w(n
t-m)表示窗函数在时间上反转并有n
t
个样本偏移量,x(m)是短时傅里叶变换的输入信号,m是输入信号x(m)的采样点,取残差r(n
t
)作为输入信号,nf表示以tf为采样周期对频率采样得到的频率总点数,nf是从1到nf的频率采样点,本实施例tf=0.488,nf=1024,tf(n
t
,nf)是大小为1024行3001列的矩阵。
[0081]
s3、采用改进的时频脊线估计算法估计tf(n
t
,nf)中具有最大能量的时频脊线向量1≤i≤3,改进的时频脊线估计算法具体步骤如下:
[0082]
s301、在时频矩阵tf(n
t
,nf)中找到具有最大能量的时间频率点(m
t
,mf),表达式如下:(m
t
,mf)=argmax|tf(n
t
,nf)|
ꢀꢀ
(5)
[0083]
其中m
t
,mf分别是tf(n
t
,nf)中最大能量点对应的时间点和频率点,记录具有最大能量的时间频率点其中1≤m
t
≤3001,1≤mf≤1024;
[0084]
s302、以时间频率点(m
t
,mf)为起始点估计时间点m
t
右侧从时间点m
t
+1到时间点n
t
对应的频率点,初始化估计过程用式(6)表示:
[0085]
其中δnf为设定的频率点增值,本实施例δnf=4,r
t
是时间点m
t
右侧的时间点,是时间点r
t
对应的频率点,r
t
从m
t
开始,每执行一次式(6),r
t
就需要加1,r
t
更新公式为r
t
=r
t
+1,直到r
t
=n
t
时结束估计过程,得到时间点m
t
右侧时频脊线向量对用线性拟合方法得到m
t
右侧时频脊线的斜率lk;
[0086]
s303、以时间频率点(m
t
,mf)为起始点,按照斜率lk估计时间点m
t
左侧从时间点m
t-1到时间点1对应的频率点,初始化估计过程表示如下:
[0087]
其中其中表示向下取整,表示向上取整,l
t
是时间点m
t
左侧的时间点,是时间点l
t
对应的频率点,l
t
从m
t
开始,每执行一次式(7),l
t
减1,l
t
更新公式为l
t
=l
t-1,直到l
t
=1时结束估计过程,得到完整的时频脊线向量
[0088]
s304、对时频脊线向量用基于最小二乘法则的拟合方法进行一阶线性拟合更新。
[0089]
s4、将步骤s3得到的作为acmd算法的初始频率,利用acmd算法分解非完全残差r(n
t
)得到对应的信号分量瞬时幅度更新的对用线性拟合得到并判断是否已经存在,更新非完全残差:
[0090]
s401、以步骤s3中得到的作为acmd算法的初始频率,对非完全残差r(n
t
)使用acmd算法分解得到分量信号以及对应的幅度更新的时频脊线向量并对时频脊线向量用基于最小二乘法则的拟合方法进行一阶线性拟合得到当前估计的分量信号i的时频脊线向量和调频斜率起始频率
[0091]
在使用非完全残差时,保留已估计分量信号的部分能量,但每个分量信号的能量大小不一致,如果非完全残差中已估计分量信号保留的部分能量要大于其他未被估计的分量信号能量,则在执行步骤s3时会得到重复的时频脊线,导致步骤s401中acmd分解非完全残差得到的信号分量重复,所以需要对s401得到的结果去重处理;
[0092]
s402、当前估计的时频脊线向量为判断是否已经存在:
[0093][0094]
其中j为正整数表示已经存在的估计分量信号,1≤j≤i-1,δif为自定义的频率增量向量,本实施例δif=[0.5,0.5,

,0.5]1×
3001
,根据下式更新非完全残差r(n
t
):
[0095]
如果式(8)成立则表示步骤s401得到的时频脊线向量重复跳转到步骤s2继续执行,如果式(8)不成立则按照以下式(10)更新阈值残差r
th
(n
t
):
[0096]
阈值残差r
th
(n
t
)只有在当前估计到的分量信号与已有的估计分量信号不重复时,
才执行更新。
[0097]
s5、根据下式判断阈值残差是否达到设定的结束阈值
[0098][0099]
如果式(11)成立,整个算法结束,得到基于非完全残差的自适应时频谱:
[0100]
其中round()表示对四舍五入取整数的函数,δ()表示dirac函数,是第i个分量的幅度估计值,是第i个分量的时频脊线向量估计值,如果式(11)不成立,则跳转到步骤s2继续执行。
[0101]
时频分析方法结束后得到如图4所示的时频图,每个线性调频信号分量的调频斜率以及起始频率的估计值如下表1所示:
[0102]
表1.实施例1中线性调频信号分量的调频斜率以及起始频率的估计值表
[0103][0104][0105]
实施例2
[0106]
如图9所示,基于实施例1,本实施例继续公开一种多分量线性调频信号时频分析方法,包括以下步骤:
[0107]
s1、在实施例1的基础上更改多分量线性调频信号各个分量信号的幅度、调频斜率以及起始频率,n(n
t
)是均值为0,方差为0.3712的高斯白噪声,信噪比snr=5,a1(n
t
)=exp(-0.03
·nt
·
t
t
),f1=90.258hz,k1=-10.5370hz/s,a2(n
t
)=exp(-0.03
·nt
·
t
t
),f2=100.621hz,k2=-9.4820hz/s,a3(n
t
)=exp(-0.09
·nt
·
t
t
),f3=140.023hz,k3=-20.678hz/s。
[0108]
s2、参照实施例1中对应步骤实施,此处不再阐述。
[0109]
s3、参照实施例1中对应步骤实施,此处不再阐述。
[0110]
s4、根据snr取值,步骤s402按照下式更新非完全残差r(n
t
):
[0111]
[0112]
其他步骤参考实施例1中对应步骤实施,此处不再阐述。
[0113]
s5、判断阈值残差是否达到设定的结束阈值
[0114][0115]
算法结束得到如图6所示的时频图,每个线性调频信号分量的调频斜率以及起始频率的估计值如下表2所示:
[0116]
表2.实施例2中线性调频信号分量的调频斜率以及起始频率的估计值表
[0117][0118][0119]
实施例3
[0120]
如图9所示,基于实施例1,本实施例继续公开一种多分量线性调频信号时频分析方法,包含以下步骤:
[0121]
s1、在实施例1的基础上更改多分量线性调频信号分量个数m=4,n(n
t
)是均值为0,方差为0.4800的高斯白噪声,信噪比snr=5,a1(n
t
)=exp(-0.03
·nt
·
t
t
),f1=35.540hz,k1=10.239hz/s,a2(n
t
)=exp(-0.03
·nt
·
t
t
),f2=50.201hz,k2=20.091hz/s,a3(n
t
)=exp(-0.09
·nt
·
t
t
),f3=62.761hz,k3=2.678hz/s,a4(n
t
)=exp(-0.06
·nt
·
t
t
),f4=80.122hz,k3=3.092hz/s。
[0122]
s2、参照实施例1中对应步骤实施,此处不再阐述。
[0123]
s3、参照实施例1中对应步骤实施,此处不再阐述。
[0124]
s4、根据snr取值,步骤s402按照下式更新非完全残差r(n
t
):
[0125][0126]
其他步骤参考实施例1中对应步骤实施,此处不再阐述。
[0127]
s5、判断阈值残差是否达到设定的结束阈值
[0128][0129]
算法结束得到如图8所示的时频图,每个线性调频信号分量的调频斜率
以及起始频率的估计值如下表3所示:
[0130]
表3.实施例3中线性调频信号分量的调频斜率以及起始频率的估计值表
[0131][0132]
在实施例1中,snr=10,利用本发明提出的一种多分量线性调频信号时频分析方法对时频域有交叉的三分量线性调频信号进行时频分析,得到的时频谱如图4所示,对比图3由stft变换得到的时频谱,本发明提出的方法提高了时频能量聚集性以及时频分辨率,在实施例1中的表1记录了三个分量信号对应的调频斜率以及起始频率真实值和估计值,每个分量信号的调频斜率估计相对误差在
±
4.1658
×
10-4

±
9.7205
×
10-4
之间,起始频率估计相对误差在
±
5.6978
×
10-4

±
8.3672
×
10-4
之间,因此利用本发明提出的一种多分量线性调频信号时频分析方法得到的每个分量信号的调频斜率以及起始频率可信度较高。在实施例2中,snr=5,利用本发明提出的一种多分量线性调频信号时频分析方法对时频域有交叉的三分量线性调频信号进行时频分析,相对实施例1该三个分量在时频域相隔较近且有两个交点,得到的时频谱如图6所示,对比图5由stft变换得到的时频谱,本发明方法明显提高了时频分辨率以及能量聚集性,在实施例2的表2中记录了三个分量对应的调频斜率以及起始频率估计值,每个分量的调频斜率估计相对误差在
±
1.7710
×
10-4

±
12.00
×
10-4
之间,起始频率估计相对误差在
±
1.8783
×
10-4

±
4.3064
×
10-4
之间,因此对每个分量的调频斜率以及起始频率估计值可信度较高。在实施例3中,snr=5,利用本发明提出的一种多分量线性调频信号时频分析方法对时频域有交叉的四分量线性调频信号进行时频分析,相对实施例1和实施例2,该实施例的多分量线性调频信号由四个线性调频信号分量组成,得到的时频谱如图8所示,对比其图7stft变换得到的时频谱,很好的提高了时频能量聚集性以及时频分辨率,在实施例3的表3中记录了四个分量对应的调频斜率以及起始频率估计值,每个分量的调频斜率估计相对误差在
±2×
10-4

±9×
10-4
,起始频率估计相对误差在
±1×
10-4

±4×
10-4
之间,因此对每个分量的调频斜率以及起始频率估计值可信度较高。综上所述,实施例1、实施例2、实施例3分别仿真验证了本发明提出的方法在信噪比srn=10以及snr=5环境下对多分量线性调频信号的时频分析,同时对不同分量个数的多分量线性调频信号进行了时频分析仿真验证,仿真结果表明本发明提出的一种多分量线性调频信号时频分析方法的有效性,能准确的估计每个分量信号的时频脊线以及调频斜率和起始频率。
[0133]
上述实施例为本发明较佳的实施方式,但本发明的实施方式并不受上述实施例的
限制,其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化,均应为等效的置换方式,都包含在本发明的保护范围之内。
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1