CN117520782A - 一种窗长自适应的短时傅里叶变换方法 - Google Patents

一种窗长自适应的短时傅里叶变换方法 Download PDF

Info

Publication number
CN117520782A
CN117520782A CN202311545976.2A CN202311545976A CN117520782A CN 117520782 A CN117520782 A CN 117520782A CN 202311545976 A CN202311545976 A CN 202311545976A CN 117520782 A CN117520782 A CN 117520782A
Authority
CN
China
Prior art keywords
window length
fourier transform
time
time point
quality evaluation
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.)
Pending
Application number
CN202311545976.2A
Other languages
English (en)
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.)
Southeast University
Original Assignee
Southeast University
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 Southeast University filed Critical Southeast University
Priority to CN202311545976.2A priority Critical patent/CN117520782A/zh
Publication of CN117520782A publication Critical patent/CN117520782A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2131—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on a transform domain processing, e.g. wavelet transform

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computational Mathematics (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Discrete Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种窗长自适应的短时傅里叶变换方法，包括以下步骤：第一步：获取待分析信号，对信号进行镜像延拓；第二步：改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换得到频谱；第三步：根据谱质量评价指标选取各时间点的初始窗长；第四步：在初始窗长的基础上进一步细化窗长间隔，改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换得到频谱；第五步：根据谱质量评价指标选取各时间点的最终窗长；第六步：将各时间点的最终窗长对应的傅里叶变换结果进行插值并拼接，得到时频图。

Description

一种窗长自适应的短时傅里叶变换方法
技术领域
本发明属于信号处理领域，尤其涉及一种窗长自适应的短时傅里叶变换方法。
背景技术
在信号处理领域，短时傅里叶变换可以同时反映信号在时域和频域的特性，是对非平稳信号进行时频分析的经典方法。由于测不准原理的约束，短时傅里叶变换结果的时间分辨率和频率分辨率之间存在矛盾。短时傅里叶变换窗长的选择决定了时频图的时间分辨率和频率分辨率：当窗长较大时，频率分辨率高而时间分辨率低；当窗长较小时，时间分辨率高而频率分辨率低。因此，短时傅里叶变换窗长的确定在信号处理领域有着重要研究价值与意义。
在传统的短时傅里叶变换中，往往根据经验选择窗长，很难得到高时频分辨率的结果。并且一般以固定窗长对整个信号进行时频分析，对于时变信号，固定的时频分辨率难以适应信号在不同时间段的特性，分析效果较差。
发明内容
发明目的：针对以上现有技术存在的问题，本发明提出一种窗长自适应的短时傅里叶变换方法。该方法基于最小谱质量评价指标准则，根据信号在不同时间段的特性自适应确定傅里叶变换的窗长。本发明突破了窗长固定的短时傅里叶变换的时频分辨率之间的约束，具有良好的时频分析性能。
技术方案：为实现上述目的，本发明提出一种窗长自适应的短时傅里叶变换方法，该方法包括以下步骤：
步骤1，获取待分析信号，对信号进行镜像延拓；
步骤2，改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换得到频谱；
步骤3，根据谱质量评价指标选取各时间点的初始窗长；
步骤4，在初始窗长的基础上进一步细化窗长间隔，改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换得到频谱；
步骤5，根据谱质量评价指标选取各时间点的最终窗长；
步骤6，将各时间点的最终窗长对应的傅里叶变换结果进行插值并拼接，以此得到时频图。
进一步的，在步骤1中，获取待分析信号x(n)，n＝1，2，…，N，n为信号采样时刻序号，N为信号采样点数，按照如下公式对信号进行镜像延拓：
其中，s(i)为镜像延拓后的信号，i＝1，2，…，3N为镜像延拓后的信号采样时刻序号。
进一步的，在步骤2中，以2的整数次幂为傅里叶变换的窗长，即L＝2m+2，其中，m＝1，2，…，M，为窗长序号，M为满足的最大整数，以K为时间间隔，分别以信号s(i)的第N+jK个时间点为中心做傅里叶变换得到频谱，其中，K为小于N的正整数，j＝1，2，…，J，为时间点序号，J为不大于N/K的最大整数。
进一步的，在步骤3中，根据谱质量评价指标选取各时间点的初始窗长，包括以下步骤：
(3.1)谱质量评价指标的计算方法如下：
其中，P(j，m)、AS(j，m)、AM(j，m)、V(j，m)和f(j，m)分别为信号s(i)以第N+jK个时间点为中心、窗长为2m+2时的傅里叶变换的频谱的谱质量评价指标、旁瓣幅度、主瓣幅度、主瓣瞬时带宽和瞬时频率，并规定V(j，m)为频谱归一化幅度的处的带宽；
(3.2)按照最小谱质量评价指标的准则选取各时间点的初始窗长为：
Linitial(j)＝2^[2+minitial(j)]
其中，
上式表示谱质量评价指标P(j，m)取最小值时的窗长序号m。
进一步的，在步骤4中，以初始窗长为基础进一步细化窗长间隔，改变窗长进行傅里叶变换，包括以下步骤：
(4.1)将各时间点的初始窗长向下向上各延伸一个档作为窗长区间，即
(4.2)将窗长区间等间隔取R个点，得到第j个时间点的窗长序列Lj＝{Lj(1)，Lj(2)，…，Lj(R)}，其中，
其中，r＝1，2，…，R为窗长序号，R为大于2的正整数，
(4.3)以信号s(i)的第N+jK个时间点为中心做傅里叶变换，窗长遍历序列Lj，将傅里叶变换结果存储为Fj(r，f)，其表示窗长为Lj(r)时的频谱，其中，f为频率。
进一步的，在步骤5中，根据谱质量评价指标选取各时间点的最终窗长，包括以下步骤：
(5.1)谱质量评价指标的计算方法如下：
其中，P(j，r)、AS(j，r)、AM(j，r)、V(j，r)和f(j，r)分别为信号s(i)以第N+jK个时间点为中心、窗长为Lj(r)时的傅里叶变换的频谱的谱质量评价指标、旁瓣幅度、主瓣幅度、主瓣瞬时带宽和瞬时频率，并规定V(j，r)为频谱归一化幅度的处的带宽；
(5.2)按照最小谱质量评价指标的准则选取各时间点的最终窗长为：
其中，
上式表示谱质量评价指标P(j，r)取最小值时的窗长序号r。
进一步的，在步骤6中，将各时间点的最终窗长对应的傅里叶变换结果进行插值并拼接得到时频图，包括以下步骤：
(6.1)提取各时间点的最终窗长对应的傅里叶变换结果为：
Sj(f)＝Fj(rfinal(j)，f)
(6.2)将上述结果进行线性插值，插值为2M+1个点，得到Sj(k)，k＝1，2，…，2M+1；
(6.3)将插值后的结果拼接为：
以此得到时频图。
有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：
在传统的短时傅里叶变换中，一般以固定窗长对信号进行时频处理，固定的时频分辨率难以适应信号在不同时间段的特性。本方法以各时间点为中心，通过粗略和细化窗长间隔的结合，两次基于最小谱质量评价指标准则选取合适的傅里叶变换窗长，减小计算量并提升算法性能，能够适应信号的时变特性，具有良好的时频分析效果。
附图说明
图1为本发明的流程框图；
图2为实例中生成的信号瞬时频率图；
图3为实例中生成的初始窗长图；
图4为实例中生成的最终窗长图；
图5为实例中生成的信号时频图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
如图1所示，本发明的一种窗长自适应的短时傅里叶变换方法，包括以下步骤：
步骤1：获取待分析信号x(n)，n＝1，2，…，N，n为信号采样时刻序号，N为信号采样点数。按照
对信号进行镜像延拓，其中，s(i)为镜像延拓后的信号，i＝1，2，…，3N为镜像延拓后的信号采样时刻序号。
步骤2：以2的整数次幂为傅里叶变换的窗长，即L＝2m+2，其中，m＝1，2，…，M为窗长序号，M为满足的最大整数。以K为时间间隔，分别以信号s(i)的第N+jK个时间点为中心做傅里叶变换，得到频谱，其中，K为小于N的正整数，j＝1，2，…，J为时间点序号，J为不大于N/K的最大整数。
步骤3：根据谱质量评价指标选取各时间点的初始窗长，包括以下步骤：
首先，计算谱质量评价指标，具体公式如下：
其中，P(j，m)、AS(j，m)、AM(j，m)、V(j，m)和f(j，m)分别为信号s(i)以第N+jK个时间点为中心、窗长为2m+2时的傅里叶变换的频谱的谱质量评价指标、旁瓣幅度、主瓣幅度、主瓣瞬时带宽和瞬时频率，并规定V(j，m)为频谱归一化幅度的处的带宽。
其次，按照最小谱质量评价指标的准则选取各时间点的初始窗长为
Linitial(j)＝2^[2+minitial(j)]
其中，
表示谱质量评价指标P(j，m)取最小值时的窗长序号m。
步骤4：在初始窗长的基础上进一步细化窗长间隔，改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换，包括以下步骤：
首先，将各时间点的初始窗长向下向上各延伸一个档作为窗长区间，即
其次，将窗长区间等间隔取R个点，得到第j个时间点的窗长序列Lj＝{Lj(1)，Lj(2)，…，Lj(R)}，其中，
其中，r＝1，2，…，R为窗长序号，R为大于2的正整数，
最后，以信号s(i)的第N+jK个时间点为中心做傅里叶变换，窗长遍历序列Lj，将傅里叶变换结果存储为Fj(r，f)，表示窗长为Lj(r)时的频谱，其中，f为频率。
步骤5：根据谱质量评价指标选取各时间点的最终窗长，包括以下步骤：
首先，计算谱质量评价指标，具体公式如下：
其中，P(j，r)、AS(j，r)、AM(j，r)、V(j，r)和f(j，r)分别为信号s(i)以第N+jK个时间点为中心、窗长为Lj(r)时的傅里叶变换的频谱的谱质量评价指标、旁瓣幅度、主瓣幅度、主瓣瞬时带宽和瞬时频率，并规定V(j，r)为频谱归一化幅度的处的带宽。
其次，按照最小谱质量评价指标的准则选取各时间点的最终窗长为
其中，
表示谱质量评价指标P(j，r)取最小值时的窗长序号r。
步骤6：将各时间点的最终窗长对应的傅里叶变换结果进行插值并拼接，得到时频图，包括以下步骤：
首先，提取各时间点的最终窗长对应的傅里叶变换结果为
Sj(f)＝Fj(rfinal(j)，f)
其次，将上述结果进行线性插值，插值为2M+1个点，得到Sj(k)，k＝1，2，…，2M+1。
最后，将插值后的结果拼接为
并得到时频图。下面例举一实施例。
实施例
现有一个信号，其由4个调频率不同的线性调频信号组成，信号的时域表达式为：
其瞬时频率如图2所示。
依据步骤1，采样得到x(n)，n＝1，2，…，N，信号采样率为Fs＝1000Hz，信号采样点数为N＝4000。对信号进行镜像延拓，得到s(i)，i＝1，2，…，12000。
依据步骤2，以2的整数次幂为傅里叶变换的窗长，即L＝2m+2，m＝1，2，…，8。以10为时间间隔，分别以信号s(i)的第4000+10j个点为中心做傅里叶变换，j＝1，2，…，400。
依据步骤3，计算谱质量评价指标，按照最小谱质量评价指标的准则选取各时间点的初始窗长Linitial(j)，如图3所示。
依据步骤4，将各时间点的初始窗长向下向上各延伸一个档作为窗长区间，将窗长区间等间隔分为15份，得到窗长序列为Lj。以信号s(i)的第N+jK个时间点为中心做傅里叶变换，窗长遍历序列Lj，将傅里叶变换结果保存。
依据步骤5，计算谱质量评价指标，按照最小谱质量评价指标的准则选取各时间点的最终窗长Lfinal(j)，如图4所示。
依据步骤6，提取各时间点的最终窗长对应的傅里叶变换结果，线性插值为512个点，然后进行拼接，并得到时频图如图5所示。
以上对本发明实施例所提供的一种窗长自适应的短时傅里叶变换方法，进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (7)

1.一种窗长自适应的短时傅里叶变换方法，其特征在于，该方法包括以下步骤：
步骤1，获取待分析信号，对信号进行镜像延拓；
步骤2，改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换得到频谱；
步骤3，根据谱质量评价指标选取各时间点的初始窗长；
步骤4，在初始窗长的基础上进一步细化窗长间隔，改变傅里叶变换的窗长，以各时间点为中心对信号做傅里叶变换得到频谱；
步骤5，根据谱质量评价指标选取各时间点的最终窗长；
步骤6，将各时间点的最终窗长对应的傅里叶变换结果进行插值并拼接以此得到时频图。
2.根据权利要求1所述的一种窗长自适应的短时傅里叶变换方法，其特征在于，在步骤1中，获取待分析信号x(n)，n＝1，2，…，N，n为信号采样时刻序号，N为信号采样点数，按照如下公式对信号进行镜像延拓：
其中，s(i)为镜像延拓后的信号，i＝1，2，…，3N为镜像延拓后的信号采样时刻序号。
3.根据权利要求2所述的一种窗长自适应的短时傅里叶变换方法，其特征在于，在步骤2中，以2的整数次幂为傅里叶变换的窗长，即L＝2m+2，其中，m＝1，2，…，M，为窗长序号，M为满足的最大整数，以K为时间间隔，分别以信号s(i)的第N+jK个时间点为中心做傅里叶变换得到频谱，其中，K为小于N的正整数，j＝1，2，…，J，为时间点序号，J为不大于N/K的最大整数。
4.根据权利要求3所述的一种窗长自适应的短时傅里叶变换方法，其特征在于，在步骤3中，根据谱质量评价指标选取各时间点的初始窗长，包括以下步骤：
(3.1)谱质量评价指标的计算方法如下：
其中，P(j，m)、AS(j，m)、AM(j，m)、V(j，m)和f(j，m)分别为信号s(i)以第N+jK个时间点为中心、窗长为2m+2时的傅里叶变换的频谱的谱质量评价指标、旁瓣幅度、主瓣幅度、主瓣瞬时带宽和瞬时频率，并规定V(j，m)为频谱归一化幅度的处的带宽；
(3.2)按照最小谱质量评价指标的准则选取各时间点的初始窗长为：
Linitial(j)＝2^[2+minitial(j)]
其中，
上式表示谱质量评价指标P(j，m)取最小值时的窗长序号m。
5.根据权利要求4所述的一种窗长自适应的短时傅里叶变换方法，其特征在于，在步骤4中，以初始窗长为基础进一步细化窗长间隔，改变窗长进行傅里叶变换，包括以下步骤：
(4.1)将各时间点的初始窗长向下向上各延伸一个档作为窗长区间，即
(4.2)将窗长区间等间隔取R个点，得到第j个时间点的窗长序列Lj＝(Lj(1)，Lj(2)，…，Lj(R)}，其中，
其中，r＝1，2，…，R为窗长序号，R为大于2的正整数，
(4.3)以信号s(i)的第N+jK个时间点为中心做傅里叶变换，窗长遍历序列Lj，将傅里叶变换结果存储为Fj(r，f)，其表示窗长为Lj(r)时的频谱，其中，f为频率。
6.根据权利要求5所述的一种窗长自适应的短时傅里叶变换方法，其特征在于，在步骤5中，根据谱质量评价指标选取各时间点的最终窗长，包括以下步骤：
(5.1)谱质量评价指标的计算方法如下：
其中，P(j，r)、AS(j，r)、AM(j，r)、V(j，r)和f(j，r)分别为信号s(i)以第N+jK个时间点为中心、窗长为Lj(r)时的傅里叶变换的频谱的谱质量评价指标、旁瓣幅度、主瓣幅度、主瓣瞬时带宽和瞬时频率，并规定V(j，r)为频谱归一化幅度的处的带宽；
(5.2)按照最小谱质量评价指标的准则选取各时间点的最终窗长为：
其中，
上式表示谱质量评价指标P(j，r)取最小值时的窗长序号r。
7.根据权利要求6所述的一种窗长自适应的短时傅里叶变换方法，其特征在于，在步骤6中，将各时间点的最终窗长对应的傅里叶变换结果进行插值并拼接得到时频图，包括以下步骤：
(6.1)提取各时间点的最终窗长对应的傅里叶变换结果为：
Sj(f)＝Fj(rfinal(j)，f)
(6.2)将上述结果进行线性插值，插值为2M+1个点，得到Sj(k)，k＝1，2，…，2M+1；
(6.3)将插值后的结果拼接为：
以此得到时频图。
CN202311545976.2A 2023-11-20 2023-11-20 一种窗长自适应的短时傅里叶变换方法 Pending CN117520782A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311545976.2A CN117520782A (zh) 2023-11-20 2023-11-20 一种窗长自适应的短时傅里叶变换方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311545976.2A CN117520782A (zh) 2023-11-20 2023-11-20 一种窗长自适应的短时傅里叶变换方法

Publications (1)

Publication Number Publication Date
CN117520782A true CN117520782A (zh) 2024-02-06

Family

ID=89743441

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311545976.2A Pending CN117520782A (zh) 2023-11-20 2023-11-20 一种窗长自适应的短时傅里叶变换方法

Country Status (1)

Country Link
CN (1) CN117520782A (zh)

Similar Documents

Publication Publication Date Title
EP2413313B1 (en) Method and device for audio signal classification
CN102323518B (zh) 一种基于谱峭度的局部放电信号识别方法
US20110286605A1 (en) Noise suppressor
RU2011105976A (ru) Устройство и способы для обработки аудио сигнала, с целью повышения разборчивости речи, используя функцию выделения нужных характеристик
CN102144258B (zh) 促进确定信号边界频率的方法和装置
CN103424183B (zh) 机械振动信号检测异常干扰消除方法
CN108267657B (zh) 一种基于s变换的电能质量扰动检测方法及系统
CN116626408B (zh) 基于机器学习的电源纹波噪声检测方法
CN113325277A (zh) 一种局部放电处理方法
CN114785379A (zh) 一种水声janus信号参数估计方法及系统
CN114492538A (zh) 一种城市中压配电电缆局部放电信号去噪方法
CN108680782B (zh) 基于极值点对称模式分解的电压闪变参数检测方法
CN114268531B (zh) 一种单音干扰检测和消除方法
CN108761202B (zh) 极点对称模态分解和希尔伯特变换相结合的谐波检测方法
Boashash et al. Time-frequency signal analysis and instantaneous frequency estimation: methodology, relationships and implementations
CN117520782A (zh) 一种窗长自适应的短时傅里叶变换方法
CN110112757B (zh) 基于sure小波消噪和改进hht的低频振荡分析方法
CN110287853B (zh) 一种基于小波分解的暂态信号去噪方法
Takagi et al. A method for pitch extraction of speech signals using autocorrelation functions through multiple window lengths
JPH05505282A (ja) デジタル・フィルターとその設計方法
CN113899976B (zh) 一种复合电能质量扰动可视化方法
Meng et al. An empirical envelope estimation algorithm
CN111241902A (zh) 一种高精度多重同步压缩广义s变换时频分析方法
CN110412522B (zh) 一种nlfm波形设计方法
CN109239458B (zh) 一种高噪声背景下电能质量扰动信号降噪方法

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