CN114866159B - 一种多分量线性调频信号时频分析方法 - Google Patents

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

Info

Publication number
CN114866159B
CN114866159B CN202210355207.5A CN202210355207A CN114866159B CN 114866159 B CN114866159 B CN 114866159B CN 202210355207 A CN202210355207 A CN 202210355207A CN 114866159 B CN114866159 B CN 114866159B
Authority
CN
China
Prior art keywords
frequency
time
component
signal
point
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210355207.5A
Other languages
English (en)
Other versions
CN114866159A (zh
Inventor
宁更新
谢靓
冯义志
张军
陈芳炯
余华
季飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN202210355207.5A priority Critical patent/CN114866159B/zh
Publication of CN114866159A publication Critical patent/CN114866159A/zh
Application granted granted Critical
Publication of CN114866159B publication Critical patent/CN114866159B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04B—TRANSMISSION
    • H04B11/00—Transmission systems employing sonic, ultrasonic or infrasonic waves
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04B—TRANSMISSION
    • H04B13/00—Transmission systems characterised by the medium used for transmission, not provided for in groups H04B3/00 - H04B11/00
    • H04B13/02—Transmission systems in which the medium consists of the earth or a large mass of water thereon, e.g. earth telegraphy
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L27/00—Modulated-carrier systems
    • H04L27/10—Frequency-modulated carrier systems, i.e. using frequency-shift keying
    • H04L27/14—Demodulator circuits; Receiver circuits
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种多分量线性调频信号时频分析方法，过程如下：首先对多分量线性调频信号进行短时傅里叶变换得到时频分布，然后用改进的时频脊线估计算法估计时频分布中具有最大能量的时频脊线，用该时频脊线作为ACMD算法的初始瞬时频率信息，通过ACMD算法分解非完全残差得到信号分量估计值以及对应的瞬时幅度和时频脊线估计值，然后更新非完全残差继续循环迭代估计下一个信号分量及对应的瞬时幅度和时频脊线，直到非完全残差达到设定的阈值，结束整个时频分析算法，根据每个分量信号的瞬时幅度和时频脊线估计值得到改进的自适应时频谱以及每个分量信号对应的调频斜率和起始频率估计值。

Description

一种多分量线性调频信号时频分析方法
技术领域
本发明涉及多分量非平稳信号时频估计技术领域，具体涉及一种多分量线性调频信号时频分析方法。
背景技术
线性调频信号是一种非平稳信号，线性调频信号在雷达探测、转轴故障分析、水下通信系统中用途较为广泛。在水下通信系统中利用多个线性调频信号进行信道探测，由于水下信道的复杂性以及各种噪声的干扰，最终在接收端会收到一个受多普勒效应影响且带有噪声以及时延的多分量线性调频信号，所以对多分量线性调频信号进行时频分析研究具有重要意义。
目前已有的时频分析方法主要分为线性时频分析和双线性时频分析，传统的线性时频分析方法有短时傅里叶变换(STFT)、小波变换(WT)，双线性时频分析方法有Wigner-Ville分布(WVD)。目前一些基于STFT改进的同步压缩算法和基于WVD改进的二次时频算法，处理时频域没有交叉的多分量线性调频信号可以得到能量聚集性较高、时频分辨率较高的时频分布，但对时频域存在交叉的多分量线性调频信号处理结果较差，在时频域交叉点处会发生时频脊线失真。针对此问题，可以用模式分解算法进行时频分析，比如基于自适应线性调频模式分解(ACMD)算法结合信号残差获得自适应时频谱，ACMD算法能较好的分解多分量线性调频信号，但处理时频域存在交叉的多分量线性调频信号时，残差的时频分布会出现时频脊线断裂问题，导致后续信号分量的时频脊线估计误差较大进而影响最终的自适应时频谱。
发明内容
本发明的目的是为了解决现有技术中多分量线性调频信号时频分析方法存在的问题，提供一种多分量线性调频信号时频分析方法
本发明的目的可以通过采取如下技术方案达到：
一种多分量线性调频信号时频分析方法，所述时频分析方法包括以下步骤：
S1、对水声通信系统接收到的多分量线性调频信号以Ts为采样周期采样得到离散的多分量线性调频信号S(nt)，并定义非完全残差R(nt)，利用非完全残差解决时频脊线断裂情况；
S2、对非完全残差R(nt)做短时傅里叶变换得到时频矩阵TF(nt,nf)；
S3、采用改进的时频脊线估计算法估计TF(nt,nf)中具有最大能量的时频脊线向量
Figure BDA0003576948390000021
i表示第i个分量，1≤i≤M，M是信号分量总个数；
S4、将步骤S3中得到的
Figure BDA0003576948390000022
作为ACMD算法的初始频率，利用ACMD算法分解非完全残差R(nt)得到对应的信号分量
Figure BDA0003576948390000023
瞬时幅度
Figure BDA0003576948390000024
和更新的
Figure BDA0003576948390000025
对更新的
Figure BDA0003576948390000026
做线性拟合得到对应的调频斜率估计值
Figure BDA0003576948390000027
起始频率估计值
Figure BDA0003576948390000028
以及时频脊线向量最终估计值
Figure BDA0003576948390000029
更新非完全残差R(nt)，并判断
Figure BDA00035769483900000210
是否已存在，根据
Figure BDA00035769483900000211
是否存在的判断结果确定是否需要更新阈值残差Rth(nt)；
S5、判断阈值残差Rth(nt)是否达到设定的结束阈值(1-α)S(nt)，如果阈值残差大于结束阈值则跳转到S2继续执行，否则整个时频分析方法结束得到各个分量信号的时频脊线向量
Figure BDA00035769483900000212
瞬时幅度
Figure BDA00035769483900000213
调频斜率
Figure BDA00035769483900000214
和起始频率
Figure BDA00035769483900000215
进一步地，所述步骤S1过程如下：
对水声通信系统接收到的多分量线性调频信号以Ts为采样周期采样得到离散的多分量线性调频信号S(nt)，并定义非完全残差表达式，离散多分量线性调频信号S(nt)的表达式如下：
Figure BDA0003576948390000031
其中Nt是S(nt)的长度，M是S(nt)的分量总个数，n(nt)为高斯白噪声，Ai(nt)是第i个分量信号的幅度，fi是第i个分量信号的起始频率，ki是第i个分量信号的调频斜率；
在基于ACMD算法结合残差得到的自适应时频谱方法中，残差定义为接收信号与当前估计到的分量信号之间的差，记作：
Figure BDA0003576948390000032
其中
Figure BDA0003576948390000033
表示分量信号i的估计值，为了处理基于ACMD算法结合残差的自适应时频谱存在的时频脊线断裂问题，在计算残差时保留当前估计到的分量信号i的部分能量，定义非完全残差R(nt)，表达式如下：
Figure BDA0003576948390000034
其中，0＜α＜1，α的取值由信噪比SNR决定，即α＝g(SNR)，噪声能量越小即信噪比SNR越大，残差的时频脊线断裂越严重，所以在信噪比大时需要保留分量信号i更多的能量以解决残差时频脊线断裂问题，因此定义g(SNR)在(-∞,+∞)上是单调递减函数，由一元线性函数、双边指数函数或argtan函数构成，非完全残差R(nt)的初始值是接收信号S(nt)，即R(nt)＝S(nt)。
进一步地，所述步骤S2过程如下：
对非完全残差R(nt)做短时傅里叶变换得到时频矩阵TF(nt,nf)，表达式如下：
Figure BDA0003576948390000041
其中w(m)是窗函数，w(nt-m)表示窗函数在时间上反转并有nt个样本偏移量，x(m)是短时傅里叶变换的输入信号，m是输入信号x(m)的采样点，取残差R(nt)作为输入信号，Nf表示以Tf为采样周期对频率采样得到的频率总点数，nf是取值从1到Nf的频率采样点。
进一步地，所述步骤S3过程如下：
根据基于ACMD算法结合残差的时频分析方法，利用与贪心算法类似的原理每次从非完全残差中分解出能带走最多能量的信号分量，执行ACMD算法时需要提供初始频率值，所以需要估计非完全残差时频分布中具有最大能量的时频脊线，具体估计步骤如下：
S301、在TF(nt,nf)中找到具有最大能量的时间频率点(mt,mf)，表达式如下：
Figure BDA0003576948390000046
其中mt、mf分别是TF(nt,nf)中最大能量点对应的时间点和频率点，记录具有最大能量的时间频率点
Figure BDA0003576948390000042
其中1≤mt≤Nt，1≤mf≤Nf；
S302、以时间频率点(mt,mf)为起始点估计时间点mt右侧从时间点mt+1到时间点Nt对应的频率点，初始化
Figure BDA0003576948390000043
估计过程用下式表示：
Figure BDA0003576948390000044
其中Δnf为设定的频率点增值，rt是时间点mt右侧的时间点，
Figure BDA0003576948390000045
是时间点rt对应的频率点，rt从mt开始，每执行一次式(6)，rt就需要加1，rt更新公式为rt＝rt+1，直到rt＝Nt时结束估计过程，得到时间点mt右侧时频脊线向量
Figure BDA0003576948390000051
对
Figure BDA0003576948390000052
用线性拟合方法得到mt右侧时频脊线的斜率Lk；
为了使得mt左右侧估计的时频脊线属于同一分量信号，时间点mt左侧时频脊线的斜率应该和mt右侧估计的时频脊线斜率一致或者接近，所以下面按照mt右侧时频脊线斜率Lk去估计时间点mt左侧的时频脊线；
S303、以时间频率点(mt,mf)为起始点，按照斜率Lk估计时间点mt左侧从时间点mt-1到时间点1对应的频率点，初始化
Figure BDA0003576948390000053
估计过程用下式表示：
Figure BDA0003576948390000054
其中
Figure BDA0003576948390000055
Δnf为设定的频率点增值
Figure BDA0003576948390000056
表示向下取整，
Figure BDA0003576948390000057
表示向上取整，lt是时间点mt左侧的时间点，
Figure BDA0003576948390000058
是时间点lt对应的频率点，lt从mt开始，每执行一次式(7)，lt减1，lt更新公式为lt＝lt-1，直到lt＝1时结束估计过程，得到完整的时频脊线向量
Figure BDA0003576948390000059
S304、对估计到的时频脊线向量
Figure BDA00035769483900000510
用基于最小二乘法则的拟合方法进行一阶线性拟合更新。
进一步地，所述步骤S4过程如下：
S401、以步骤S3中得到的
Figure BDA00035769483900000511
作为ACMD算法的初始频率，对非完全残差R(nt)使用ACMD算法分解得到分量信号
Figure BDA00035769483900000512
对应的瞬时幅度
Figure BDA00035769483900000513
更新的时频脊线向量
Figure BDA00035769483900000514
并对时频脊线向量
Figure BDA00035769483900000515
用基于最小二乘法则的拟合方法进行一阶线性拟合得到当前估计的分量信号i的时频脊线向量
Figure BDA0003576948390000061
和调频斜率
Figure BDA0003576948390000062
起始频率
Figure BDA0003576948390000063
在使用非完全残差时，保留已估计分量信号的部分能量，但每个分量信号的能量大小不一致，如果非完全残差中已估计分量信号保留的部分能量要大于其他未被估计的分量信号能量，则在执行步骤S3时会得到重复的时频脊线，导致步骤S401中ACMD分解非完全残差得到的信号分量重复，所以需要对S401得到的结果去重处理；
S402、当前估计的时频脊线向量为
Figure BDA0003576948390000064
判断
Figure BDA0003576948390000065
是否已经存在：
Figure BDA0003576948390000066
其中j为正整数表示已经得到的估计分量信号，1≤j≤i-1，ΔIF为自定义的频率增量向量，根据下式更新非完全残差R(nt)，表达式如下：
Figure BDA0003576948390000067
如果式(8)成立则表示步骤S401得到的时频脊线向量
Figure BDA0003576948390000068
重复直接跳转到步骤S2，如果式(8)不成立则按照以下式(10)更新阈值残差Rth(nt)：
Figure BDA0003576948390000069
阈值残差Rth(nt)只有在当前估计到的分量信号i与已有的估计分量信号j不重复时，才执行更新。
进一步地，所述步骤S5过程如下：
无论时频脊线是否重复估计都需要更新非完全残差，所以不能根据非完全残差值判定算法是否结束，步骤S402定义了新的阈值残差用来判定算法是否结束；根据式(11)判断阈值残差是否达到设定的结束阈值(1-α)S(nt)：
Rth(nt)≤(1-α)S(nt)  (11)
如果式(11)成立，得到基于非完全残差的自适应时频谱：
Figure BDA0003576948390000071
其中round()表示四舍五入取整数的函数，δ()表示Dirac函数，
Figure BDA0003576948390000072
是分量i的幅度估计值，
Figure BDA0003576948390000073
是分量i的时频脊线估计值，如果式(11)不成立，则跳转到步骤S2继续执行。
进一步地，所述步骤S5中当式(11)成立，利用已经得到的分量信号i对应的瞬时幅度估计值
Figure BDA0003576948390000074
和时频脊线估计值
Figure BDA0003576948390000075
以及理想时频分布原理得到基于非完全残差的自适应时频谱的同时，得到每个分量信号对应的调频斜率估计值
Figure BDA0003576948390000076
和起始频率估计值
Figure BDA0003576948390000077
本发明相对于现有技术具有如下的优点及效果：
1、该方法适应性广。本发明不需要知道多分量线性调频信号的先验信息信号分量总个数，通过阈值残差和设定的结束阈值大小关系判定整个时频分析算法是否结束，同时对估计的时频脊线进行去重处理，避免得到重复的时频脊线。
2、该方法实现步骤简单。本发明通过定义的非完全残差解决基于ACMD算法结合残差的时频分析方法处理时频域有交叉的多分量线性调频信号时存在的时频脊线断裂问题，提高时频脊线估计以及时频谱准确度。
3、该方法能获得每个线性调频信号分量较为精确的起始频率和调频斜率估计值。本发明利用改进的时频脊线估计算法避免时频分布中最大能量点右侧时频脊线向量和左侧时频脊线向量不属于同一信号分量，用基于最小二乘法则的线性拟合方法对时频脊线向量进行一阶线性拟合得到对应的起始频率和调频斜率。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是本发明实施1中信噪比10dB下完全去除分量1之后的残差的STFT变换结果；
图2是本发明实施例1中信噪比10dB下不完全去除分量1之后的非完全残差STFT变换结果；
图3是本发明实施例1中信噪比是10dB下的第一个时频脊线估计结果图；
图4是本发明实施例1中信噪比是10dB下的基于非完全残差得到的自适应时频图；
图5是本发明实施例2信噪比是5dB下的多分量线性调频信号STFT变换结果仿真图；
图6是本发明实施例2信噪比是5dB下的多分量线性调频信号基于非完全残差得到的自适应时频图；
图7是本发明实施例3信噪比是5dB下的多分量线性调频信号STFT变换结果仿真图；
图8是本发明实施例3信噪比是5dB下的多分量线性调频信号基于非完全残差得到的自适应时频图；
图9是本发明中公开的一种多分量线性调频信号时频分析方法的流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合实施例对本发明作进一步的详细描述，但本发明的实施方式不限于此。
实施例1
如图9所示，本实施例公开了一种多分量线性调频信号时频分析方法，包含以下步骤：
S1、对水声通信系统接收到的多分量线性调频信号以Ts＝0.002s为采样周期采样得到离散多分量线性调频信号S(nt)，并定义非完全残差表达式，S(nt)表达式如下：
Figure BDA0003576948390000091
其中Nt＝3001，M＝3，i＝1,2,3，n(nt)是均值为0，方差为0.1189的高斯白噪声，信噪比SNR＝10，A1(nt)＝exp(-0.03·nt·Tt)，f1＝25.098Hz，k1＝15.123Hz/s，A2(nt)＝exp(-0.03·nt·Tt)，f2＝40.8921Hz，k2＝8.23Hz/s，A3(nt)＝exp(-0.09·nt·Tt)，f3＝60.023Hz，k3＝12.678Hz/s。
定义非完全残差，其表达式如下：
Figure BDA0003576948390000092
其中
Figure BDA0003576948390000093
表示分量信号i的估计值，1≤i≤M，0＜α＜1，α的取值由信噪比SNR决定，即α＝g(SNR)，g(SNR)在(-∞,+∞)上是单调递减函数，由一元线性函数、双边指数函数或argtan函数构成，非完全残差R(nt)的初始值是S(nt)，即R(nt)＝S(nt)。本实施例中α定义为：
Figure BDA0003576948390000094
S2、对非完全残差R(nt)做短时傅里叶变换得到时频矩阵TF(nt,nf)，表达式如下：
Figure BDA0003576948390000095
其中w(m)是窗函数，w(nt-m)表示窗函数在时间上反转并有nt个样本偏移量，x(m)是短时傅里叶变换的输入信号，m是输入信号x(m)的采样点，取残差R(nt)作为输入信号，Nf表示以Tf为采样周期对频率采样得到的频率总点数，nf是从1到Nf的频率采样点，本实施例Tf＝0.488，Nf＝1024，TF(nt,nf)是大小为1024行3001列的矩阵。
S3、采用改进的时频脊线估计算法估计TF(nt,nf)中具有最大能量的时频脊线向量
Figure BDA0003576948390000101
1≤i≤3，改进的时频脊线估计算法具体步骤如下：
S301、在时频矩阵TF(nt,nf)中找到具有最大能量的时间频率点(mt,mf)，表达式如下：(mt,mf)＝argmax|TF(nt,nf)|  (5)
其中mt，mf分别是TF(nt,nf)中最大能量点对应的时间点和频率点，记录具有最大能量的时间频率点
Figure BDA0003576948390000102
其中1≤mt≤3001，1≤mf≤1024；
S302、以时间频率点(mt,mf)为起始点估计时间点mt右侧从时间点mt+1到时间点Nt对应的频率点，初始化
Figure BDA0003576948390000103
估计过程用式(6)表示：
Figure BDA0003576948390000104
其中Δnf为设定的频率点增值，本实施例Δnf＝4，rt是时间点mt右侧的时间点，
Figure BDA0003576948390000105
是时间点rt对应的频率点，rt从mt开始，每执行一次式(6)，rt就需要加1，rt更新公式为rt＝rt+1，直到rt＝Nt时结束估计过程，得到时间点mt右侧时频脊线向量
Figure BDA0003576948390000106
对
Figure BDA0003576948390000107
用线性拟合方法得到mt右侧时频脊线的斜率Lk；
S303、以时间频率点(mt,mf)为起始点，按照斜率Lk估计时间点mt左侧从时间点mt-1到时间点1对应的频率点，初始化
Figure BDA0003576948390000108
估计过程表示如下：
Figure BDA0003576948390000111
其中
Figure BDA0003576948390000112
Figure BDA0003576948390000113
表示向下取整，
Figure BDA0003576948390000114
表示向上取整，lt是时间点mt左侧的时间点，
Figure BDA0003576948390000115
是时间点lt对应的频率点，lt从mt开始，每执行一次式(7)，lt减1，lt更新公式为lt＝lt-1，直到lt＝1时结束估计过程，得到完整的时频脊线向量
Figure BDA0003576948390000116
S304、对时频脊线向量
Figure BDA0003576948390000117
用基于最小二乘法则的拟合方法进行一阶线性拟合更新。
S4、将步骤S3得到的
Figure BDA0003576948390000118
作为ACMD算法的初始频率，利用ACMD算法分解非完全残差R(nt)得到对应的信号分量
Figure BDA0003576948390000119
瞬时幅度
Figure BDA00035769483900001110
更新的
Figure BDA00035769483900001111
对
Figure BDA00035769483900001112
用线性拟合得到
Figure BDA00035769483900001113
并判断
Figure BDA00035769483900001114
是否已经存在，更新非完全残差：
S401、以步骤S3中得到的
Figure BDA00035769483900001115
作为ACMD算法的初始频率，对非完全残差R(nt)使用ACMD算法分解得到分量信号
Figure BDA00035769483900001116
以及对应的幅度
Figure BDA00035769483900001117
更新的时频脊线向量
Figure BDA00035769483900001118
并对时频脊线向量
Figure BDA00035769483900001119
用基于最小二乘法则的拟合方法进行一阶线性拟合得到当前估计的分量信号i的时频脊线向量
Figure BDA00035769483900001120
和调频斜率
Figure BDA00035769483900001121
起始频率
Figure BDA00035769483900001122
在使用非完全残差时，保留已估计分量信号的部分能量，但每个分量信号的能量大小不一致，如果非完全残差中已估计分量信号保留的部分能量要大于其他未被估计的分量信号能量，则在执行步骤S3时会得到重复的时频脊线，导致步骤S401中ACMD分解非完全残差得到的信号分量重复，所以需要对S401得到的结果去重处理；
S402、当前估计的时频脊线向量为
Figure BDA00035769483900001123
判断
Figure BDA00035769483900001124
是否已经存在：
Figure BDA0003576948390000121
其中j为正整数表示已经存在的估计分量信号，1≤j≤i-1，ΔIF为自定义的频率增量向量，本实施例ΔIF＝[0.5,0.5,…,0.5]1×3001，根据下式更新非完全残差R(nt)：
Figure BDA0003576948390000122
如果式(8)成立则表示步骤S401得到的时频脊线向量
Figure BDA0003576948390000123
重复跳转到步骤S2继续执行，如果式(8)不成立则按照以下式(10)更新阈值残差Rth(nt)：
Figure BDA0003576948390000124
阈值残差Rth(nt)只有在当前估计到的分量信号与已有的估计分量信号不重复时，才执行更新。
S5、根据下式判断阈值残差是否达到设定的结束阈值
Figure BDA0003576948390000125
如果式(11)成立，整个算法结束，得到基于非完全残差的自适应时频谱：
Figure BDA0003576948390000126
其中round()表示对
Figure BDA0003576948390000127
四舍五入取整数的函数，δ()表示Dirac函数，
Figure BDA0003576948390000128
是第i个分量的幅度估计值，
Figure BDA0003576948390000129
是第i个分量的时频脊线向量估计值，如果式(11)不成立，则跳转到步骤S2继续执行。
时频分析方法结束后得到如图4所示的时频图，每个线性调频信号分量的调频斜率
Figure BDA00035769483900001210
以及起始频率
Figure BDA00035769483900001211
的估计值如下表1所示：
表1.实施例1中线性调频信号分量的调频斜率以及起始频率的估计值表
Figure BDA00035769483900001212
Figure BDA0003576948390000131
实施例2
如图9所示，基于实施例1，本实施例继续公开一种多分量线性调频信号时频分析方法，包括以下步骤：
S1、在实施例1的基础上更改多分量线性调频信号各个分量信号的幅度、调频斜率以及起始频率，n(nt)是均值为0，方差为0.3712的高斯白噪声，信噪比SNR＝5，A1(nt)＝exp(-0.03·nt·Tt)，f1＝90.258Hz，k1＝-10.5370Hz/s，A2(nt)＝exp(-0.03·nt·Tt)，f2＝100.621Hz，k2＝-9.4820Hz/s，A3(nt)＝exp(-0.09·nt·Tt)，f3＝140.023Hz，k3＝-20.678Hz/s。
S2、参照实施例1中对应步骤实施，此处不再阐述。
S3、参照实施例1中对应步骤实施，此处不再阐述。
S4、根据SNR取值，步骤S402按照下式更新非完全残差R(nt)：
Figure BDA0003576948390000132
其他步骤参考实施例1中对应步骤实施，此处不再阐述。
S5、判断阈值残差是否达到设定的结束阈值
Figure BDA0003576948390000133
Figure BDA0003576948390000134
算法结束得到如图6所示的时频图，每个线性调频信号分量的调频斜率
Figure BDA0003576948390000135
以及起始频率
Figure BDA0003576948390000136
的估计值如下表2所示：
表2.实施例2中线性调频信号分量的调频斜率以及起始频率的估计值表
Figure BDA0003576948390000137
Figure BDA0003576948390000141
实施例3
如图9所示，基于实施例1，本实施例继续公开一种多分量线性调频信号时频分析方法，包含以下步骤：
S1、在实施例1的基础上更改多分量线性调频信号分量个数M＝4，n(nt)是均值为0，方差为0.4800的高斯白噪声，信噪比SNR＝5，A1(nt)＝exp(-0.03·nt·Tt)，f1＝35.540Hz，k1＝10.239Hz/s，A2(nt)＝exp(-0.03·nt·Tt)，f2＝50.201Hz，k2＝20.091Hz/s，A3(nt)＝exp(-0.09·nt·Tt)，f3＝62.761Hz，k3＝2.678Hz/s，A4(nt)＝exp(-0.06·nt·Tt)，f4＝80.122Hz，k3＝3.092Hz/s。
S2、参照实施例1中对应步骤实施，此处不再阐述。
S3、参照实施例1中对应步骤实施，此处不再阐述。
S4、根据SNR取值，步骤S402按照下式更新非完全残差R(nt)：
Figure BDA0003576948390000142
其他步骤参考实施例1中对应步骤实施，此处不再阐述。
S5、判断阈值残差是否达到设定的结束阈值
Figure BDA0003576948390000143
Figure BDA0003576948390000144
算法结束得到如图8所示的时频图，每个线性调频信号分量的调频斜率
Figure BDA0003576948390000145
以及起始频率
Figure BDA0003576948390000146
的估计值如下表3所示：
表3.实施例3中线性调频信号分量的调频斜率以及起始频率的估计值表
Figure BDA0003576948390000151
在实施例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环境下对多分量线性调频信号的时频分析，同时对不同分量个数的多分量线性调频信号进行了时频分析仿真验证，仿真结果表明本发明提出的一种多分量线性调频信号时频分析方法的有效性，能准确的估计每个分量信号的时频脊线以及调频斜率和起始频率。
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (6)

1.一种多分量线性调频信号时频分析方法，其特征在于，所述时频分析方法包括以下步骤：
S1、对水声通信系统接收到的多分量线性调频信号以Ts为采样周期采样得到离散的多分量线性调频信号S(nt)，并定义非完全残差R(nt)解决基于ACMD算法结合残差的自适应时频分析方法中存在的时频脊线断裂问题；
S2、对非完全残差R(nt)做短时傅里叶变换得到时频矩阵TF(nt,nf)；
S3、采用改进的时频脊线估计算法估计TF(nt,nf)中具有最大能量的时频脊线向量
Figure FDA0004036424340000011
i表示第i个分量，1≤i≤M；
S4、将步骤S3中得到的
Figure FDA0004036424340000012
作为ACMD算法的初始频率，利用ACMD算法分解非完全残差R(nt)得到对应的信号分量
Figure FDA0004036424340000013
瞬时幅度
Figure FDA0004036424340000014
和更新的
Figure FDA0004036424340000015
对更新的
Figure FDA0004036424340000016
做线性拟合得到对应的调频斜率估计值
Figure FDA0004036424340000017
起始频率估计值
Figure FDA0004036424340000018
以及时频脊线向量最终估计值
Figure FDA0004036424340000019
更新非完全残差R(nt)，并判断
Figure FDA00040364243400000110
是否已存在，根据
Figure FDA00040364243400000111
是否存在的判断结果确定是否需要更新阈值残差Rth(nt)；
S5、判断阈值残差Rth(nt)是否达到设定的结束阈值(1-α)S(nt)，如果阈值残差大于结束阈值则跳转到S2继续执行，否则整个时频分析方法结束得到各个分量信号的时频脊线向量
Figure FDA00040364243400000112
瞬时幅度
Figure FDA00040364243400000113
调频斜率
Figure FDA00040364243400000114
和起始频率
Figure FDA00040364243400000115
其中，所述步骤S1过程如下：
对水声通信系统接收到的多分量线性调频信号以Ts为采样周期采样得到离散的多分量线性调频信号S(nt)，并定义非完全残差表达式，多分量线性调频信号S(nt)的表达式如下：
Figure FDA0004036424340000021
其中Nt是S(nt)的长度，M是S(nt)的分量总个数，n(nt)为高斯白噪声，Ai(nt)是第i个分量信号的幅度，fi是第i个分量信号的起始频率，ki是第i个分量信号的调频斜率；
接收信号与当前估计到的分量信号之间的差记作：
Figure FDA0004036424340000022
其中
Figure FDA0004036424340000023
表示第i个分量信号的估计值，为了处理基于ACMD算法结合残差的自适应时频谱存在的时频脊线断裂问题，在计算残差时保留当前估计到的第i个分量信号的部分能量，定义非完全残差R(nt)，表达式如下：
Figure FDA0004036424340000024
其中，0＜α＜1，α的取值由信噪比SNR决定，即α＝g(SNR)，g(SNR)在(-∞,+∞)上是单调递减函数，由一元线性函数、双边指数函数或argtan函数构成，非完全残差R(nt)的初始值是接收信号S(nt)，即R(nt)＝S(nt)。
2.根据权利要求1所述的一种多分量线性调频信号时频分析方法，其特征在于，所述步骤S2过程如下：
对非完全残差R(nt)做短时傅里叶变换得到时频矩阵TF(nt,nf)，表达式如下：
Figure FDA0004036424340000025
其中w(m)是窗函数，w(nt-m)表示窗函数在时间上反转并有nt个样本偏移量，x(m)是短时傅里叶变换的输入信号，m是输入信号x(m)的采样点，取残差R(nt)作为输入信号，Nf表示以Tf为采样周期对频率采样得到的频率总点数，nf是取值从1到Nf的频率采样点。
3.根据权利要求2所述的一种多分量线性调频信号时频分析方法，其特征在于，所述步骤S3过程如下：
S301、在TF(nt,nf)中找到具有最大能量的时间频率点(mt,mf)，表达式如下：
Figure FDA0004036424340000031
其中mt是TF(nt,nf)中最大能量点对应的时间点，mf是TF(nt,nf)中最大能量点对应的频率点，记录具有最大能量的时间频率点
Figure FDA0004036424340000032
其中1≤mt≤Nt，1≤mf≤Nf；
S302、以时间频率点(mt,mf)为起始点估计时间点mt右侧从时间点mt+1到时间点Nt对应的频率点，初始化
Figure FDA0004036424340000033
估计过程用下式表示：
Figure FDA0004036424340000034
其中Δnf为设定的频率点增值，rt是时间点mt右侧的时间点，
Figure FDA0004036424340000035
是时间点rt对应的频率点，rt从mt开始，每执行一次式(6)，rt就需要加1，rt更新公式为rt＝rt+1，直到rt＝Nt时结束估计过程，得到时间点mt右侧时频脊线向量
Figure FDA0004036424340000036
对
Figure FDA0004036424340000037
用线性拟合方法得到mt右侧时频脊线的斜率Lk；
S303、以时间频率点(mt,mf)为起始点，按照斜率Lk估计时间点mt左侧从时间点mt-1到时间点1对应的频率点，初始化
Figure FDA0004036424340000038
估计过程用下式表示：
Figure FDA0004036424340000039
其中
Figure FDA0004036424340000041
Δnf为设定的频率点增值，
Figure FDA0004036424340000042
表示向下取整，
Figure FDA0004036424340000043
表示向上取整，lt是时间点mt左侧的时间点，
Figure FDA0004036424340000044
是时间点lt对应的频率点，lt从mt开始，每执行一次式(7)，lt减1，lt更新公式为lt＝lt-1，直到lt＝1时结束估计过程，得到完整的时频脊线向量
Figure FDA0004036424340000045
S304、对步骤S303得到的时频脊线向量
Figure FDA0004036424340000046
用基于最小二乘法则的拟合方法进行一阶线性拟合更新。
4.根据权利要求3所述的一种多分量线性调频信号时频分析方法，其特征在于，所述步骤S4过程如下：
S401、以步骤S3中得到的
Figure FDA0004036424340000047
作为ACMD算法的初始频率，对非完全残差R(nt)使用ACMD算法分解得到分量信号
Figure FDA0004036424340000048
对应的瞬时幅度
Figure FDA0004036424340000049
更新的时频脊线向量
Figure FDA00040364243400000410
并对时频脊线向量
Figure FDA00040364243400000411
用基于最小二乘法则的拟合方法进行一阶线性拟合得到当前估计到的信号分量i的时频脊线向量
Figure FDA00040364243400000412
和调频斜率
Figure FDA00040364243400000413
起始频率
Figure FDA00040364243400000414
S402、当前估计的时频脊线向量为
Figure FDA00040364243400000415
判断
Figure FDA00040364243400000416
是否已经存在：
Figure FDA00040364243400000417
其中j为正整数表示已经得到的估计分量信号，1≤j≤i-1，ΔIF为自定义的频率增量向量，根据下式更新非完全残差R(nt)，表达式如下：
Figure FDA00040364243400000418
如果式(8)成立则表示步骤S401得到的时频脊线向量
Figure FDA00040364243400000419
重复直接跳转到步骤S2，如果式(8)不成立则按照式(10)更新阈值残差Rth(nt)：
Figure FDA00040364243400000420
阈值残差Rth(nt)只有在当前估计到的分量信号i与已有的估计分量信号j不重复时，才执行更新。
5.根据权利要求4所述的一种多分量线性调频信号时频分析方法，其特征在于，所述步骤S5过程如下：
根据式(11)判断阈值残差是否达到设定的结束阈值(1-α)S(nt)：
Rth(nt)≤(1-α)S(nt) (11)
如果式(11)成立，得到基于非完全残差的自适应时频谱：
Figure FDA0004036424340000051
其中round()表示四舍五入取整数的函数，δ()表示Dirac函数，
Figure FDA0004036424340000052
是第i个分量信号的幅度估计值，
Figure FDA0004036424340000053
是第i个分量信号的时频脊线估计值，如果式(11)不成立，则跳转到步骤S2继续执行。
6.根据权利要求5所述的一种多分量线性调频信号时频分析方法，其特征在于，所述步骤S5中当式(11)成立，得到基于非完全残差的自适应时频谱的同时，得到每个分量信号对应的调频斜率估计值
Figure FDA0004036424340000054
和起始频率估计值
Figure FDA0004036424340000055
CN202210355207.5A 2022-04-01 2022-04-01 一种多分量线性调频信号时频分析方法 Active CN114866159B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210355207.5A CN114866159B (zh) 2022-04-01 2022-04-01 一种多分量线性调频信号时频分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210355207.5A CN114866159B (zh) 2022-04-01 2022-04-01 一种多分量线性调频信号时频分析方法

Publications (2)

Publication Number Publication Date
CN114866159A CN114866159A (zh) 2022-08-05
CN114866159B true CN114866159B (zh) 2023-04-21

Family

ID=82629467

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210355207.5A Active CN114866159B (zh) 2022-04-01 2022-04-01 一种多分量线性调频信号时频分析方法

Country Status (1)

Country Link
CN (1) CN114866159B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113472390A (zh) * 2021-07-07 2021-10-01 哈尔滨工程大学 一种基于深度学习的跳频信号参数估计方法
CN113887360A (zh) * 2021-09-23 2022-01-04 同济大学 一种基于迭代扩展频散模态分解的频散波提取方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107248869B (zh) * 2017-05-24 2019-05-14 电子科技大学 一种基于吕分布的多分量线性调频信号去噪方法
KR102033831B1 (ko) * 2017-12-05 2019-11-08 한양대학교 에리카산학협력단 Eq 기법을 이용한 수동 시역전 수중통신 성능 예측 방법
CN109510787B (zh) * 2018-10-15 2021-08-17 中国人民解放军战略支援部队信息工程大学 水声信道下线性调频信号参数估计方法及装置
CN111665469B (zh) * 2020-06-11 2022-08-23 浙江大学 一种基于空间时频分布的水下多径信号参数估计方法
CN111935038B (zh) * 2020-08-03 2022-08-19 中国人民解放军国防科技大学 基于分数阶傅里叶变换的线性调频干扰消除方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113472390A (zh) * 2021-07-07 2021-10-01 哈尔滨工程大学 一种基于深度学习的跳频信号参数估计方法
CN113887360A (zh) * 2021-09-23 2022-01-04 同济大学 一种基于迭代扩展频散模态分解的频散波提取方法

Also Published As

Publication number Publication date
CN114866159A (zh) 2022-08-05

Similar Documents

Publication Publication Date Title
US7835456B2 (en) Spur mitigation techniques
US8130852B2 (en) Method for estimating channel in radio communication system and device therefor
CN107666451B (zh) 用于lte系统的信道估计方法
CN111935046B (zh) 一种低复杂度的频移键控信号符号率估计方法
CN114785379A (zh) 一种水声janus信号参数估计方法及系统
US7499497B2 (en) Impulse suppression apparatus applied in OFDM system and related method thereof
CN114866159B (zh) 一种多分量线性调频信号时频分析方法
Liu et al. Ziv–Zakai time-delay estimation bounds for frequency-hopping waveforms under frequency-selective fading
CN109167744B (zh) 一种相位噪声联合估计方法
WO2020183544A1 (ja) 受信装置、無線通信システムおよび干渉電力推定方法
CN111490954A (zh) 信道脉冲响应的重要时延抽头选择方法及系统
US20020196877A1 (en) Method and apparatus for characterizing the distortion produced by electronic equipment
CN112883787B (zh) 一种基于频谱匹配的短样本低频正弦信号参数估计方法
CN114884777A (zh) 一种基于变换域的信道估计方法
Tamim et al. Hilbert transform of FFT pruned cross correlation function for optimization in time delay estimation
Li et al. Interrupted-sampling repeater jamming (ISRJ) suppression based on cyclostationarity
CN109302360B (zh) 信道估计方法及装置、计算机可读存储介质、终端
US11336317B2 (en) Radio communication system, interference suppression method, control circuit, and program storage medium
Roienko et al. An Overview of the Adaptive Robust DFT and it’s Applications
CN113721201B (zh) 一种线性调频信号调频率的估计方法
RU2812822C1 (ru) Способ выделения полезной составляющей из входного сигнала, содержащего полезную составляющую и шум
CN115037388B (zh) 一种基于改进梯度下降法的lfm信号瞬时频率提取方法
US20240171301A1 (en) Machine learning device and jamming signal generation apparatus
CN110943952B (zh) 调幅信号的检测方法及装置
Fritz et al. Analyzing the Impact of HF-Specific Signal Degradation on Automatic Speech Recognition

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant