短时傅里叶变换与WVD相结合的自适应时频分析方法

短时傅里叶变换与wvd相结合的自适应时频分析方法
技术领域
1.本发明属于信号处理技术领域，尤其涉及一种短时傅里叶变换与wvd相 结合的自适应时频分析方法。


背景技术：

2.非合作脉冲信号时频分析是在没有任何先验知识的情况下，通过信号处 理算法实现信号时频分布的计算，观察信号的时间和频率特征。非合作水声 脉冲信号通常是单频脉冲信号、线性调频信号以及由这两种信号构成的组合 脉冲信号，对几种组合脉冲信号进行高质量时频分析，同时观察出水声脉冲 信号的准确的时频关系，可为实现截获敌方的发射主动脉冲信号提供技术支 持，这在海洋军事战争中对判断敌方目标意图以及获取军事情报有着极为重 要的意义。
3.目前国内外学者提出了许多脉冲信号时频分析方法，主要有：(1)在线 性时频分析中，最为经典的短时傅里叶变换(short time fourier transform， stft)，它是根据gabor提出的在傅里叶变换上加窗而得到局部傅里叶变换 的思想形成的，stft算法简单，容易实现，但其时频分辨率无法同时实现； (2)继承了stft的局部化思想，小波变换(wavelet transform，wt)应运 而生，它克服了stft窗口大小不随频率变化的缺点，弥补了时频分辨率之间 的矛盾，但冗余度和计算量较大；(3)在二次型时频分析中，wigner-ville 分布(wvd)作为最核心部分，通过其非线性特性，打破了stft时频分辨率 的矛盾，具有很高的时频分辨率，也具有良好的时频聚集性和时频边缘特性， 但也正因为其非线性特性，在应用于多分量信号时会出现交叉项问题。


技术实现要素：

4.本发明目的在于提供一种短时傅里叶变换与wvd相结合的自适应时频分 析方法,以解决脉冲信号常规时频分析方法时频分辨率无法同时提高，二次型 时频分析方法冗余度和计算量较大，且在应用于多分量信号时会出现交叉项 的技术问题。
5.为解决上述技术问题，本发明的具体技术方案如下：
6.一种短时傅里叶变换与wvd相结合的自适应时频分析方法，包括如下步 骤：
7.步骤1、获取待分析脉冲信号采样数据序列x(n),n＝0,1,
…
,n-1，所述的n 为待分析脉冲信号脉宽长度所对应的采样点个数；
8.步骤2、对自适应分析所需时间和频率参数进行初始化；
9.步骤3、依据所述采样数据序列x(n)计算脉冲信号功率谱p(k)；
10.步骤4、依据脉冲信号功率谱p(k)估计信号的频率下限f
l
和频率上限fh， 并对短时傅里叶变换的分析窗长l和窗长的步进值r进行初始化；
11.步骤5、依据自适应生成的f
l
，fh，l和r计算采样数据序列的归一化短 时傅里叶变换yo(λ,q)；
12.步骤6、计算采样数据序列的wvdwo(n,l)，并选取其在频带f
l
～fh内的归 一化wvdw
(n,l)；
13.步骤7、依据归一化wvdw(n,l)的维数，对归一化短时傅里叶变换yo(λ,q) 进行插值得到映射短时傅里叶变换y(n,l)；
14.步骤8、将映射短时傅里叶变换结果y(n,l)与归一化wvd矩阵w(n,l)经过 处理后进行矩阵点乘，得到归一化自适应时频分析结果z(n,l)。
15.进一步的，所述步骤1中获取待分析的脉冲信号采样数据序列x(n)包括 以下步骤：
16.从传感器接收n个采样点的实时采集数据作为待分析的采样数据序列 x(n)。
17.进一步的，所述步骤1中获取待分析的脉冲信号采样数据序列x(n)包括 以下步骤：
18.从存储器中提取包含整个脉冲信号的n个采样点数据作为待分析的采 样数据序列x(n)，n为大于等于128的整数。
19.进一步的，所述步骤2中对自适应分析所需时间和频率参数进行初始化 包括以下步骤：
20.设置功率变化率谱阈值系数ξ，频率分辨率控制系数α和窗长的步进值 控制系数η，其中ξ、α和η为正实数，0.02≤ξ≤0.04，0.01≤α≤0.03，0.05≤η≤0.2。
21.进一步的，所述步骤3中，依据采样数据序列x(n)计算脉冲信号功率谱 p(k)包括以下步骤：
22.对所述采样数据序列x(n)做快速傅里叶变换，计算得到采样数据序列的 离散傅里叶变换x(d)和脉冲信号功率谱p(k)，包括如下步骤：
23.步骤3.1、计算x(n)的离散傅里叶变换为
[0024][0025]
其中d为x(d)的离散频率索引，j表示虚数单位，即该式通过 快速傅里叶变换实现；
[0026]
步骤3.2、依据x(d)计算x(n)的脉冲信号功率谱
[0027][0028]
其中k为p(k)的离散频率索引，| |代表取模值运算，起始离散频率索引 为k
l
，终止离散频率索引为kh，此时有k
l
＝1，kh＝n/2，离散频率索引表示为 k＝k
l
,k
l
1,
…
,kh。
[0029]
进一步的，所述步骤4中依据脉冲信号功率谱p(k)估计信号的频率下限 f
l
和频率上限fh，并对短时傅里叶变换的分析窗长l和窗长的步进值r进行初 始化包括以下步骤：
[0030]
根据脉冲信号功率谱p(k)计算信号功率变化率谱p(k)，以功率变化率谱 阈值p0为门界搜索p(k)的两个突变点对应的离散频率索引k'
l
和k'h，k'
l
作为估 计的离散频率索引下限赋值给k
l
，k'h作为估计的离散频率索引上限赋值给kh， 估计出f
l
与fh，并根据f
l
与fh初始化合适的l与r，包括如下步骤：
[0031]
步骤4.1、计算脉冲信号功率谱p(k)的变化率谱
[0032][0033]
其中δfn＝fs/n为离散傅里叶变换所对应的频率分辨率，fs为采样数据序 列的采样频率。设置功率变化率谱阈值p0＝ξmax(p(k))，max()为取序列或 矩阵中的最大值；
[0034]
步骤4.2、以p0为门界搜索p(k)的正序第一个突变点的离散频率索引k'
l
， 搜索过程包括如下步骤：
[0035]
步骤4.2.1、初始化起始搜索离散频率索引k＝k
l
；
[0036]
步骤4.2.2、判断p(k)＞p0或k＝kh是否成立，若成立则跳到步骤4.2.4， 否则进入步骤4.2.3；
[0037]
步骤4.2.3、令k＝k 1，并返回步骤4.2.2；
[0038]
步骤4.2.4、令k'
l
＝k，得到p(k)正序第一个突变点对应的离散频率索引 k'
l
；
[0039]
步骤4.3、以p0为门界搜索p(k)的倒序第一个突变点的离散频率索引k'h， 搜索过程包括如下步骤：
[0040]
步骤4.3.1、初始化起始搜索离散频率索引k＝kh；
[0041]
步骤4.3.2、判断p(k)＞p0或k＝k
l
是否成立，若成立则跳到步骤4.3.4， 否则进入步骤4.3.3；
[0042]
步骤4.3.3令k＝k-1，并返回步骤4.3.2；
[0043]
步骤4.3.4令k'h＝k，得到p(k)倒序第一个突变点对应的离散频率索引k'h；
[0044]
步骤4.4、将两个突变点对应的离散频率索引k'
l
和k'h赋值给离散频率索引 下限和上限并赋值给k
l
＝k'
l
和kh＝k'h；
[0045]
步骤4.5、采样数据序列的采样频率为fs，采样数据序列的频率分辨率 δfn＝fs/n，δfn为正实数，根据得到的新的k
l
和kh，计算出精细信号频率下 限f'
l
＝δfn(k
l-1)和精细信号频率上限f'h＝δfn(k
h-1)，并对f'
l
和f'h进行修正得 到估计的信号频率f
l
＝f'
l-(f'
h-f'
l
)/10和频率上限fh＝f'h (f'
h-f'
l
)/10，f
l
和fh为 正实数；
[0046]
步骤4.6、计算出合适的采样数据序列的频率分辨率为
[0047]
δf
t
＝αmin((f
l
fh)/2,|f
h-f
l
|)
[0048]
其中min(,)表示求两个数的较小数，α为频率分辨率控制系数，其取值范 围为0.01≤α≤0.03，根据得到的δf
t
计算出自适应的短时傅里叶变换分析窗长 l＝round(fs/δf
t
)，窗长的步进值r＝round(ηl)，其中round()代表四舍五入取 整，l为不小于20的正整数，r为正整数。
[0049]
进一步的，所述在步骤5中依据自适应生成的f
l
，fh，l和r计算采样数 据序列的归一化短时傅里叶变换yo(λ,k)，包括如下步骤：
[0050]
步骤5.1、依据自适应生成的f
l
，fh，l和r，计算x(n)在频率从f
l
至fh的短时傅里叶变换为
[0051][0052]
其中lc为短时傅里叶变换的离散频率索引长度，lc为不小于2的正整数， h(μ
λ
)为
长度为l的时频分析窗，μ
λ
为分析窗截取的采样数据序列时间索引， λ为分析窗序号，q是短时傅里叶变换的离散频率索引，起始离散频率索引 为q
l
，终止离散频率索引为qh，具体计算过程包括如下步骤：
[0053]
步骤5.1.1、单个分析窗内的采样数据序列的频率分辨率δf＝fs/l，δf为 正实数，起始离散频率索引为q
l
＝round(f
l
/δf) 1，终止离散频率索引为 qh＝round(fh/δf) 1，q
l
和qh为正整数，短时傅里叶变换的离散频率索引 q＝q
l
,q
l
1,...,qh，短时傅里叶变换的离散频率索引长度lc＝q
h-q
l
1，lc为不 小于2的正整数，分析窗个数为
[0054]
m＝fix((n-l)/r)-1
[0055]
其中fix()为向零方向取整，分析窗序号λ＝1,2,...,m，m为不小于2的 正整数；
[0056]
步骤5.1.2、对λ逐步赋值从1到m，每次对λ的赋值结束后对分析窗截 取的采样数据序列起始索引μ
λl
和终止索引μ
λh
进行赋值，μ
λl
＝(λ-1)r 1， μ
λh
＝(λ-1)r lc，μ
λl
和μ
λh
为正整数，μ
λ
＝μ
λl
,
…
,μ
λh
，μ
λ
的长度为lc，接着对 x(μ
λ
)序列做离散傅里叶变换，得到该序列的离散傅里叶变换为
[0057][0058]
q＝q
l
,q
l
1,...,qhλ＝1,2,...,m
[0059]
计算x(μ
λ
)的离散傅里叶变换x(λ,q)可通过快速傅里叶变换实现，
[0060]
进而得到一个m行lc列的短时傅里叶变换矩阵
[0061][0062]
步骤5.2、对短时傅里叶变换矩阵ys(λ,q)做归一化处理后，得到归一化 的短时傅里叶变换yo(λ,q)，即yo(λ,q)＝ys(λ,q)/max(ys(λ,q))。
[0063]
进一步的，所述步骤6中计算采样数据序列的wvdwo(n,l)，并选取其在 频带f
l
～fh内的归一化wvdw(n,l)具体包括以下步骤：
[0064]
步骤6.1、对采样数据序列x(n)做希尔伯特变换，得到复信号数据序列 y(n)；
[0065]
步骤6.2、计算采样数据序列的wvdwo(n,l)，计算过程包括如下步骤：
[0066]
步骤6.2.1、设半延时索引τm＝0,1,...,k-1，k为半延时索引序列长度， 初始化k＝round(n/8)，k为不小于16的正整数，wvd的频率分辨率 δl＝fs/(2k)，δl为正实数，wvd的离散频率索引为l；
[0067]
步骤6.2.2、初始化采样数据序列索引n＝0；
[0068]
步骤6.2.3、令正向半延时索引为τ1＝n,n 1,...,min(n k-1,n)，反向半延 时索引为τ2＝n,n-1,...,max(n-k 1,0)，计算两个索引序列长度较短的一个为公 共半延时索引长度为
[0069]kl
＝min[(min(n k-1,n)-n 1),(n-max(n-k 1,0) 1)]
[0070]
正向半延时索引τ1＝n,n 1,...,n k
l-1，反向半延时索引 τ2＝n,n-1,...,n-k
l
1，对两个索引序列取公共索引p＝0,1,...,k
l-1，其中τ1(0)＝n， τ1(1)＝n 1，τ2(0)＝n，τ2(1)＝n-1，以此类推；
[0071]
步骤6.2.4、计算采样数据序列索引n对应的瞬时自相关函数为
[0072]kz
(τ1,τ2)＝y(τ1(p))y
*
(τ2(p)),p＝0,1,...,k
l-1
[0073]
其中*表示取共轭；
[0074]
步骤6.2.5、若k
l
＜k，将kz(τ1,τ2)补零至长度为k，即将两个索引序列 公共索引补值为p＝0,1,...,k-1，对p＞k
l-1，kz(τ1,τ2)＝0；若k
l
＝k，则进入 步骤6.2.6；
[0075]
步骤6.2.6、对kz(τ1,τ2)做离散傅里叶变换，得到瞬时自相关函数的离散 傅里叶变换为
[0076][0077]
计算kz(τ1,τ2)的离散傅里叶变换w(l)可通过快速傅里叶变换实现，
[0078]
进而得到该n取值下的关于离散频率索引l的wvd为
[0079][0080]
步骤6.2.7、判断n＝n-1是否成立，若成立跳到步骤6.2.9，否则进入 步骤6.2.8；
[0081]
步骤6.2.8、令n＝n 1，并返回步骤6.2.3；
[0082]
步骤6.2.9、wvd的离散频率索引上限l1＝round(f
l
/δl) 1，离散频率索引 下限l2＝round(fh/δl) 1，l1和l2为正整数，令离散频率索引为l＝l1,l1 1,...,l2， 该索引序列长度kc＝l
2-l1 1，kc为不小于2的正整数，得到一个n行kc列的 在自适应频带内的wvd矩阵wo(n,l)；
[0083]
步骤6.3、计算wvd矩阵的最大值wm＝max(wo(n,l))，进而计算出归一化 wvdw(n,l)＝wo(n,l)/wm。
[0084]
进一步的，所述步骤7中对归一化短时傅里叶变换yo(λ,q)进行插值得到 映射短时傅里叶变换y(n,l)具体包括以下步骤：对m行lc列的yo(λ,q)的每一 行做一维线性插值，得到m行k列的列插值归一化短时傅里叶变换矩阵 y
t
(λ,l)，再对y
t
(λ,l)的每一列做一维线性插值，最终得到与w(n,l)同等大小的 n行kc列的映射短时傅里叶变换y(n,l)。
[0085]
进一步的，所述步骤8中将映射短时傅里叶变换y(n,l)与归一化wvd w(n,l)经过处理后进行矩阵点乘得到归一化的自适应时频分析结果z(n,l)具 体包括以下步骤：对y(n,l)和w(n,l)两个矩阵的所有元素做加1取对数处理， 得到对数处理后的归一化短时傅里叶变换矩阵z1(n,l)和对数处理后的归一化 wvd矩阵z2(n,l)，将两个对数处理后的矩阵各个元素对应相乘得到zo(n,l)， 最终得到归一化自适应时频分析结果z(n,l)，包括如下步骤：
[0086]
步骤8.1、对y(n,)l和w(n,l)两个矩阵进行算数处理，得到 z1(n,l)＝lg(y(n,l) 1)，z2(n,l)＝lg(w(n,l) 1)，其中lg()表示以10为底取对数， 两个等式均对矩阵的元素做加1和取对数处理；
[0087]
步骤8.2、根据对数处理后的矩阵z1(n,l)和z2(n,l)得到自适应时频分布 为
[0088]zo
(n,l)＝z1(n,l).*z2(n,l)
[0089]
其中.*表示矩阵各个元素对应相乘；
[0090]
步骤8.3、计算得到归一化自适应时频分析结果 z(n,l)＝zo(n,l)/max(zo(n,l))。
[0091]
本发明的短时傅里叶变换与wvd相结合的自适应时频分析方法，具有以 下优点：
[0092]
1、本发明的时频分析方法在求取信号的短时傅里叶变换时，利用了信号 的功率分布谱特性，除去了噪声信号的频带，自适应选择了待截获信号的频 带进行短时傅里叶变换，提高了短时傅里叶变换的时频分辨率，在步骤3和 步骤4中，通过计算采样数据序列的功率谱，通过功率变化率谱阈值来滤除 噪声信号频带，自适应地得到待截获信号的频带；
[0093]
2、本发明的时频分析方法在求取信号的短时傅里叶变换和wvd时，利用 得到的待截获信号的频带，对分析窗长和窗长步进值进行自适应生成，使得 短时傅里叶变换结果能够清晰地分辨出单频脉冲信号的频率，在步骤4中， 通过得到的待截获信号的频带，计算出合适的短时傅里叶变换频率分辨率， 进而自适应地得到分析窗长和窗长步进值，在步骤5中基于得到的分析窗长、 窗长步进值和待截获信号的频带，能够计算出较为理想的信号的短时傅里叶 变换结果；
[0094]
3、本发明的时频分析方法在求取自适应时频分析结果时，利用短时傅里 叶变换和wvd两个方法的特性，可在较低信噪比下，得到弱交叉项、高时频 分辨率的时频分布，在步骤6中，计算得到在待截获信号频带内的wvd矩阵， 在步骤7中，将短时傅里叶变换结果插值成与wvd矩阵同等大小，在步骤8 中将两个矩阵线性点乘并做运算处理，最终得到理想的时频分析结果，通过 两个矩阵点乘解决了短时傅里叶变换时频分辨率无法同时提高和wvd应用于 多分量信号时会出现交叉项的技术问题，同时涉及的运算为线性运算，避免 了二次型时频分析冗余度和计算量较大的问题。
附图说明
[0095]
图1为本发明的短时傅里叶变换与wvd相结合的自适应时频分析方法流 程示意图；
[0096]
图2为本发明的实例1仿真脉冲信号归一化短时傅里叶变换图结构示意 图；
[0097]
图3为实施例1仿真脉冲信号归一化wvd图；
[0098]
图4为实施例1仿真脉冲信号归一化自适应时频分析结果图；
[0099]
图5为实施例2仿真脉冲信号归一化短时傅里叶变换图；
[0100]
图6为实施例2仿真脉冲信号归一化wvd图；
[0101]
图7为实施例2仿真脉冲信号归一化自适应时频分析结果图。
具体实施方式
[0102]
为了更好地了解本发明的目的、结构及功能，下面结合附图，对本发明 一种短时傅里叶变换与wvd相结合的自适应时频分析方法做进一步详细的描 述。
[0103]
如图1所示，本发明一种短时傅里叶变换与wvd相结合的自适应时频分 析方法，包括如下步骤：
[0104]
步骤1、获取待分析脉冲信号采样数据序列x(n),n＝0,1,
…
,n-1，所述的n 为待分析脉冲信号脉宽长度所对应的采样点个数；
[0105]
获取待分析的脉冲信号采样数据序列包括以下步骤：从传感器接收n个 采样点的实时采集数据作为待分析的数据序列x(n),n＝0,1,
…
,n-1；或从存储 器中提取包含整个脉冲信号的n个采样点数据作为待分析的数据序列 x(n),n＝0,1,
…
,n-1，所述的n为检测到的脉冲信号脉宽长度所对应的采样点 个数，n为大于等于128的整数。
[0106]
步骤2、对自适应分析所需时间和频率参数进行初始化，包括以下步骤：
[0107]
设置功率变化率谱阈值系数ξ，频率分辨率控制系数α和窗长的步进值 控制系数η，其中ξ、α和η为正实数，0.02≤ξ≤0.04，0.01≤α≤0.03，0.05≤η≤0.2。
[0108]
为了兼顾本发明的运算量和时频分布的准确度，优选的ξ取值为0.03，α 取值为0.02，η取值为0.1。
[0109]
步骤3、依据所述采样数据序列x(n)计算脉冲信号功率谱p(k)，包括以 下步骤：
[0110]
对所述采样数据序列x(n)做快速傅里叶变换，计算得到采样数据序列的 离散傅里叶变换x(d)和脉冲信号功率谱p(k)，包括如下步骤：
[0111]
步骤3.1、计算x(n)的离散傅里叶变换为
[0112][0113]
其中d为x(d)的离散频率索引，j表示虚数单位，即该式可以 直接计算实现，但运算量较大。本发明计算x(n)的离散傅里叶变换x(d)通过 快速傅里叶变换实现。
[0114]
步骤3.2、依据x(d)计算x(n)的脉冲信号功率谱
[0115][0116]
其中k为p(k)的离散频率索引，| |代表取模值运算，起始离散频率索引 为k
l
，终止离散频率索引为kh，此时有k
l
＝1，kh＝n/2，离散频率索引表示为 k＝k
l
,k
l
1,
…
,kh。
[0117]
步骤4、依据脉冲信号功率谱p(k)估计信号的频率下限f
l
和频率上限fh， 并对短时傅里叶变换的分析窗长l和窗长的步进值r进行初始化，包括以下 步骤：
[0118]
根据脉冲信号功率谱p(k)计算信号功率变化率谱p(k)，以功率变化率谱 阈值p0为门界搜索p(k)的两个突变点对应的离散频率索引k'
l
和k'h，k'
l
作为估 计的离散频率索引下限赋值给k
l
，k'h作为估计的离散频率索引上限赋值给kh， 估计出f
l
与fh，并根据f
l
与fh初始化合适的l与r，包括如下步骤：
[0119]
步骤4.1、计算脉冲信号功率谱p(k)的变化率谱
[0120][0121]
其中δfn＝fs/n为离散傅里叶变换所对应的频率分辨率，fs为采样数据序 列的采样频率。设置功率变化率谱阈值p0＝ξmax(p(k))，max()为取序列或 矩阵中的最大值；
[0122]
步骤4.2、以p0为门界搜索p(k)的正序第一个突变点的离散频率索引k'
l
， 搜索过程包括如下步骤：
[0123]
步骤4.2.1、初始化起始搜索离散频率索引k＝k
l
；
[0124]
步骤4.2.2、判断p(k)＞p0或k＝kh是否成立，若成立则跳到步骤4.2.4， 否则进入步骤4.2.3；
[0125]
步骤4.2.3、令k＝k 1，并返回步骤4.2.2；
[0126]
步骤4.2.4、令k'
l
＝k，得到p(k)正序第一个突变点对应的离散频率索引 k'
l
；
[0127]
步骤4.3、以p0为门界搜索p(k)的倒序第一个突变点的离散频率索引k'h， 搜索过程包括如下步骤：
[0128]
步骤4.3.1、初始化起始搜索离散频率索引k＝kh；
[0129]
步骤4.3.2、判断p(k)＞p0或k＝k
l
是否成立，若成立则跳到步骤4.3.4， 否则进入步骤4.3.3；
[0130]
步骤4.3.3令k＝k-1，并返回步骤4.3.2；
[0131]
步骤4.3.4令k'h＝k，得到p(k)倒序第一个突变点对应的离散频率索引k'h；
[0132]
步骤4.4、将两个突变点对应的离散频率索引k'
l
和k'h赋值给离散频率索引 下限和上限并赋值给k
l
＝k'
l
和kh＝k'h；
[0133]
步骤4.5、采样数据序列的采样频率为fs，采样数据序列的频率分辨率 δfn＝fs/n，δfn为正实数，根据得到的新的k
l
和kh，计算出精细信号频率下 限f'
l
＝δfn(k
l-1)和精细信号频率上限f'h＝δfn(k
h-1)，并对f'
l
和f'h进行修正得 到估计的信号频率f
l
＝f'
l-(f'
h-f'
l
)/10和频率上限fh＝f'h (f'
h-f'
l
)/10，f
l
和fh为 正实数；
[0134]
步骤4.6、计算出合适的采样数据序列的频率分辨率为
[0135]
δf
t
＝αmin((f
l
fh)/2,|f
h-f
l
|)
[0136]
其中min(,)表示求两个数的较小数，α为频率分辨率控制系数，其取值范 围为0.01≤α≤0.03，根据得到的δf
t
计算出自适应的短时傅里叶变换分析窗长 l＝round(fs/δf
t
)，窗长的步进值r＝round(ηl)，其中round()代表四舍五入取 整，l为不小于20的正整数，r为正整数。
[0137]
步骤5、依据自适应生成的f
l
，fh，l和r计算采样数据序列的归一化短 时傅里叶变换yo(λ,q)，包括如下步骤：
[0138]
步骤5.1、依据自适应生成的f
l
，fh，l和r，计算x(n)在频率从f
l
至fh的短时傅里叶变换为
[0139][0140]
其中lc为短时傅里叶变换的离散频率索引长度，lc为不小于2的正整数，h(μ
λ
)为长度为l的时频分析窗，h(μ
λ
)可选择矩形分析窗，μ
λ
为分析窗截取 的采样数据序列时间索引，λ为分析窗序号，q是短时傅里叶变换的离散频 率索引，起始离散频率索引为q
l
，终止离散频率索引为qh，具体计算过程包 括如下步骤：
[0141]
步骤5.1.1、单个分析窗内的采样数据序列的频率分辨率δf＝fs/l，δf为 正实数，起始离散频率索引为q
l
＝round(f
l
/δf) 1，终止离散频率索引为 qh＝round(fh/δf) 1，q
l
和qh为正整数，短时傅里叶变换的离散频率索引 q＝q
l
,q
l
1,...,qh，短时傅里叶变换的离散频率索引长度lc＝q
h-q
l
1，lc为不 小于2的正整数，分析窗个数为
[0142]
m＝fix((n-l)/r)-1
[0143]
其中fix()为向零方向取整，分析窗序号λ＝1,2,...,m，m为不小于2的 正整数；
[0144]
步骤5.1.2、对λ逐步赋值从1到m，每次对λ的赋值结束后对分析窗截 取的采样数据序列起始索引μ
λl
和终止索引μ
λh
进行赋值，μ
λl
＝(λ-1)r 1， μ
λh
＝(λ-1)r lc，μ
λl
和μ
λh
为正整数，μ
λ
＝μ
λl
,...,μ
λh
，μ
λ
的长度为lc，接着对 x(μ
λ
)序列做离散傅里叶变换，得到该序列的离散傅里叶变换为
[0145]
＝round(fh/δl) 1，l1和l2为正整数，令离散频率索引为l＝l1,l1 1,...,l2， 该索引序列长度kc＝l
2-l1 1，kc为不小于2的正整数，得到一个n行kc列的 在自适应频带内的wvd矩阵wo(n,l)；
[0171]
步骤6.3、计算wvd矩阵的最大值wm＝max(wo(n,l))，进而计算出归一化 wvdw(n,l)＝wo(n,l)/wm。
[0172]
步骤7、依据归一化wvdw(n,l)的维数，对归一化短时傅里叶变换yo(λ,q) 进行插值得到映射短时傅里叶变换y(n,l)，具体包括以下步骤：对m行lc列 的yo(λ,q)的每一行做一维线性插值，得到m行k列的列插值归一化短时傅里 叶变换矩阵y
t
(λ,l)，再对y
t
(λ,l)的每一列做一维线性插值，最终得到与w(n,l) 同等大小的n行kc列的映射短时傅里叶变换y(n,l)。
[0173]
步骤8、将映射短时傅里叶变换结果y(n,l)与归一化wvd矩阵w(n,l)经过 处理后进行矩阵点乘，得到归一化自适应时频分析结果z(n,l)，具体包括以 下步骤：对y(n,l)和w(n,l)两个矩阵的所有元素做加1取对数处理，得到对数 处理后的归一化短时傅里叶变换矩阵z1(n,l)和对数处理后的归一化wvd矩阵 z2(n,l)，将两个对数处理后的矩阵各个元素对应相乘得到zo(n,l)，最终得到 归一化自适应时频分析结果z(n,l)，包括如下步骤：
[0174]
步骤8.1、对y(n,)l和w(n,l)两个矩阵进行算数处理，得到 z1(n,l)＝lg(y(n,l) 1)，z2(n,l)＝lg(w(n,l) 1)，其中lg()表示以10为底取对数， 两个等式均对矩阵的元素做加1和取对数处理；
[0175]
步骤8.2、根据对数处理后的矩阵z1(n,l)和z2(n,l)得到自适应时频分布 为
[0176]zo
(n,l)＝z1(n,l).*z2(n,l)
[0177]
其中.*表示矩阵各个元素对应相乘；
[0178]
步骤8.3、计算得到归一化自适应时频分析结果 z(n,l)＝zo(n,l)/max(zo(n,l))。
[0179]
本发明的实施例中，仿真接收脉冲信号模型为：
[0180][0181]
b1＜t1＜b1 t1,b2＜t2＜b2 t2,0＜t＜t
[0182]
其中a为信号幅度，为初始相位，t为脉冲宽度，b1为单频脉冲信号的 起始时刻，b2为线性调频脉冲信号的起始时刻，t1为单频脉冲信号的脉冲宽 度，t2为线性调频脉冲信号的脉冲宽度，f0为单频脉冲信号的中心频率，f1为 线性调频脉冲信号的起始频率，k0＝(f
2-f1)/t2为线性脉冲调频信号调频率，f2为线性调频脉冲信号的终止频率。w(t)为均值为0，方差为σ2的大小由信噪 比snr决定：snr＝10lg(a2/2σ2)。
[0183]
以采样频率fs对上述脉冲信号进行离散采样可得到脉冲信号采样数据 序列：
[0184][0185][0186]
其中 n
t
＝round(fst)，round()代表四舍五入运算。
[0187]
实施例1：
[0188]
仿真信号参数分别设置为：信号幅度a＝1，初始相位脉冲宽度t＝4 s，单频脉冲信号起始时刻b1＝1.7s，线性调频脉冲信号起始时刻b2＝0.2s，单 频脉冲信号脉冲宽度t1＝2s，线性调频脉冲信号脉冲宽度t2＝3s，单频脉冲信 号中心频率f0＝475hz，线性调频脉冲信号起始频率f1＝400hz，线性调频脉冲 信号终止频率f2＝450hz，线性脉冲调频信号调频率k0＝16.67hz/s，即该仿真 信号为单频脉冲信号和线性调频脉冲信号的组合信号，采样频率fs＝4000hz， 观测数据序列点数n＝16000，信噪比snr＝0db。
[0189]
在步骤2中，设置功率变化率谱阈值系数ξ＝0.03，频率分辨率控制系数 α＝0.02，窗长的步进值控制系数η＝0.1。
[0190]
依据步骤3，计算采样数据序列x(n)的功率谱p(k)。
[0191]
依据步骤4，根据脉冲信号功率谱p(k)计算信号的频率下限f
l
＝379.84hz 和频率上限fh＝511.56hz，并对短时傅里叶变换的分析窗长l＝1518和窗长的 步进值r＝152。
[0192]
依据步骤5，根据自适应生成的f
l
，fh，l和r计算数据序列的归一化短 时傅里叶变换yo(λ,q)，如图2所示，可以看出：基于自适应生成的待截获信 号的频带、分析窗长和窗长步进值计算得到短时傅里叶变换结果较为理想， 时域分辨率与频域分辨率达到平衡，但都不高。
[0193]
依据步骤6，计算数据序列在频带f
l
～fh内的wvdwo(n,l)，并计算归一化 wvdw(n,l)，如图3所示，可以看出wvd有明显的交叉项存在。
[0194]
依据步骤7，根据归一化wvdw(n,l)的维数，对归一化短时傅里叶变换 yo(λ,q)进行插值得到映射短时傅里叶变换y(n,l)。
[0195]
依据步骤8，将映射短时傅里叶变换y(n,l)与归一化wvdw(n,l)经过处理 后进行矩阵点乘得到归一化的自适应时频分析结果z(n,l)，如图4所示，得 到了高质量的信号时频分布，且结果呈现出的单频脉冲信号与线性调频脉冲 信号所设置的时频特征一致，将图2与图4对比，可以看出：最终得到的时 频分析结果与短时傅里叶变换结果相比，其时频分辨率得到了明显地提高， 解决了短时傅里叶变换时频分辨率无法同时提高的技术问题；将图3与图4 对比，可以看出：最终得到的时频分析结果与wvd结果相比，其交叉项得到 了显著的抑制，解决了wvd应用于多分量信号时会出现交叉项的技术问题。
[0196]
实施例2：
[0197]
仿真信号参数分别设置为：信号幅度a＝2，初始相位脉冲宽度 t＝2s，单频脉冲信号起始时刻b1＝0.2s，线性调频脉冲信号起始时刻b2＝0.5s， 单频脉冲信号脉冲宽度t1＝1s，线性调频脉冲信号脉冲宽度t2＝1s，单频脉冲 信号中心频率f0＝750hz，线性调频脉冲信号起始频率f1＝850hz，线性调频脉 冲信号终止频率f2＝800hz，线性脉冲调频信号调频率k0＝50hz/s，即该仿真 信号为单频脉冲信号和线性调频脉冲信号的组合信号，采样频率fs＝4000hz， 观测数据序列点数n＝8000，信噪比snr＝-3db。
[0198]
在步骤2中，设置功率变化率谱阈值系数ξ＝0.03，频率分辨率控制系数 α＝0.02，窗长的步进值控制系数η＝0.1。
[0199]
依据步骤3，计算采样数据序列x(n)的功率谱p(k)。
[0200]
依据步骤4，根据脉冲信号功率谱p(k)计算信号的频率下限f
l
＝707.83hz 和频率上限fh＝884.70hz，并对短时傅里叶变换的分析窗长l＝1131和窗长的 步进值r＝113。
[0201]
依据步骤5，根据自适应生成的f
l
，fh，l和r计算数据序列的归一化短 时傅里叶变
换yo(λ,q)，如图5所示，可以看出：基于自适应生成的待截获信 号的频带、分析窗长和窗长步进值计算得到短时傅里叶变换结果较为理想， 时域分辨率与频域分辨率达到平衡，但都不高。
[0202]
依据步骤6，计算数据序列在频带f
l
～fh内的wvdwo(n,l)，并计算归一化 wvdw(n,l)，如图6所示，可以看出wvd有明显的交叉项存在。
[0203]
依据步骤7，根据归一化wvdw(n,l)的维数，对归一化短时傅里叶变换 yo(λ,q)进行插值得到映射短时傅里叶变换y(n,l)。
[0204]
依据步骤8，将映射短时傅里叶变换y(n,l)与归一化wvdw(n,l)经过处理 后进行矩阵点乘得到归一化的自适应时频分析结果z(n,l)，如图7所示，得 到了高质量的信号时频分布，且结果呈现出的单频脉冲信号与线性调频脉冲 信号所设置的时频特征一致，将图5与图7对比，可以看出：最终得到的时 频分析结果与短时傅里叶变换结果相比，其时频分辨率得到了明显地提高， 解决了短时傅里叶变换时频分辨率无法同时提高；将图5与图7对比，可以 看出：最终得到的时频分析结果与wvd结果相比，其交叉项得到了显著的抑 制，解决了wvd应用于多分量信号时会出现交叉项的技术问题。
[0205]
可以理解，本发明是通过一些实施例进行描述的，本领域技术人员知悉 的，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行 各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例 进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此， 本发明不受此处所公开的具体实施例的限制，所有落入本技术的权利要求范 围内的实施例都属于本发明所保护的范围内。