CN107402326A - 一种改进s变换的有限窗长时频分析方法 - Google Patents

一种改进s变换的有限窗长时频分析方法 Download PDF

Info

Publication number
CN107402326A
CN107402326A CN201710593238.3A CN201710593238A CN107402326A CN 107402326 A CN107402326 A CN 107402326A CN 201710593238 A CN201710593238 A CN 201710593238A CN 107402326 A CN107402326 A CN 107402326A
Authority
CN
China
Prior art keywords
mrow
window
frequency
mfrac
time
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.)
Granted
Application number
CN201710593238.3A
Other languages
English (en)
Other versions
CN107402326B (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.)
Nanjing University of Science and Technology
Original Assignee
Nanjing University of Science and Technology
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 Nanjing University of Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN201710593238.3A priority Critical patent/CN107402326B/zh
Publication of CN107402326A publication Critical patent/CN107402326A/zh
Application granted granted Critical
Publication of CN107402326B publication Critical patent/CN107402326B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Complex Calculations (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种改进S变换的有限窗长时频分析方法，在频率与窗长之间建立基于反正切函数的映射，使得窗长在一定范围内随频率自适应变化，将改进后的窗函数代入S变换，实现高分辨率时频分析，包括以下步骤：信号采样，得到离散序列后进行快速傅里叶变换获得信号频谱，再根据信号频谱特点以及分辨率要求确定控制因子，将得到的控制函数代入窗函数中，窗函数经过快速傅里叶变换后与扩维后的信号频谱相乘，再经傅里叶逆变换得到时频谱。本方法既实现了窗长随信号频率变化而变化，又限定了窗长的变化范围，保证了时频谱各个部分都有很高的分辨率。

Description

一种改进S变换的有限窗长时频分析方法
技术领域
本发明属于一种雷达信号处理技术，特别是一种改进S变换的有限窗长高精度时频分析方法。
背景技术
非平稳信号是雷达信号处理中最常见的信号，而时频分析是分析该类信号的重要工具。为了准确分析信号的局部特性，时频分析将一维时域信号映射到二维时频平面，从而获得信号的时频分布。目前，常用的时频分析方法主要包括短时傅里叶变换(STFT)、希尔伯特-黄变换(HHT)、S变换(ST)、小波变换(WT)等。
Dennis Gabor于1946年提出了短时傅里叶变换，其基本思想是通过加窗实现信号的分段傅里叶变换，从而得到信号的时变特性。但STFT所用窗函数固定，与时间和频率无关，是一种单一分辨率的分析方法。1998年，黄锷等人提出的希尔伯特-黄变换，将经验模态分解与Hilbert谱分析相结合，将信号分解为若干固有模态函数(IMF)后映射到时频域中，但其理论与算法还未完善，存在模态混叠、端点效应等问题。而小波变换的思想来源于伸缩与平移方法，是一种窗口面积固定但形状可改变的时频局部化分析方法，能根据高低频信号特点自适应调整时频窗，有着“数学显微镜”美称。但小波基设计难度较大，还有容许性条件的约束，同时存在时频分辨率不足、尺度频率转换复杂等缺陷。
为了弥补短时傅里叶变换和小波变换的不足，Stockwell提出了S变换，引入了可变高斯窗函数，且时窗宽度与频率导数成反比。该方法得到的时频谱在低频部分频率分辨率高，在高频部分频率分辨率低，即分辨率可变。但是这种反比关系使得窗函数在局部出现窗长过宽和过窄的问题，导致低频处时间定位失效，高频处频率定位失效。
专利申请号为CN201611158226.X，发明名称为“一种基于非线性模式分解和自适应最优核的时频分析方法”的中国专利，首先运用非线性模式分解算法，将多分量非平稳信号分解为一组有物理意义的非线性模式分量，再利用自适应最优核的时频分析方法使核函数自适应地随着信号的变化而变化。该方法可以改善时频聚焦能力，局限性在于自适应最优核得到比较困难。
专利申请号为CN201610946585.5，发明名称为“一种基于反褶积广义S变换的地震频谱成像方法”的中国专利，通过将原始信号与高斯窗各自的魏格纳分布进行二维褶积得到时频谱。该方法可以压制Wigner-Ville分布的交叉项的产生，同时使广义S变换谱获得了较高的时频聚集性，但局限性在于低频和高频处的分辨率不足问题得不到解决。
由上可知，现有的时频分析方法还存在不足，需进一步改进来实现高精度的时频分析。
发明内容
本发明所解决的技术问题在于提供一种改进S变换的有限窗长时频分析方法。
实现本发明目的的技术解决方案为：一种改进S变换的有限窗长时频分析方法，包括以下步骤：
步骤1、对信号进行采样，采样点数为N＝t/T,其中t为信号时长，T为采样周期，得到信号的离散序列h[kT](k＝1,2,…,N)；
步骤2、对信号的离散序列h[kT]进行快速傅里叶变换，得到离散信号频谱其中(n＝0,1,2,…,N-1)；
步骤3、确定控制因子a、b、c的值，得到窗长控制函数后确定窗函数；
步骤4、对窗函数进行快速傅里叶变化得到窗函数频谱；
步骤5、将扩维后的信号频谱跟窗函数频谱相乘，再对其进行傅里叶逆变换；
步骤6、重复步骤4、步骤5，直至所有频率点全部计算完成，得到高分辨率时频谱。
本发明与现有技术相比，其显著优点为：1)窗函数可以根据频率自适应做出调整，在低频处有较高的频率分辨率，在高频处有较高的时间分辨率；2)通过反正切函数限制了窗长的变化范围，解决了在低频处由于时窗宽度过宽引起的时间定位不准确问题，以及在高频处由于频窗宽度过宽引起的频率定位不准确问题；3)对于不同类型的信号和不同的分辨率要求，通过调整控制因子取值，都能实现高分辨率时频分析，具有很强的灵活性。
下面结合附图对本发明作进一步详细描述。
附图说明
图1是本发明一种改进S变换的有限窗长时频分析方法流程图。
图2是本发明一种改进S变换的有限窗长时频分析方法的算法流程图。
图3是本发明实施例1的时频分析结果图。
图4是本发明实施例2的时频分析结果图。
图5是本发明实施例3的时频分析结果图。
具体实施方式
结合图1和图2所示，本发明的一种改进S变换的有限窗长时频分析方法，包括以下步骤：
步骤1、对信号进行采样，采样点数为N＝t/T,其中t为信号时长，T为采样周期，得到信号的离散序列h[kT](k＝1,2,…,N)；
步骤2、对信号的离散序列h[kT]进行快速傅里叶变换，得到离散信号频谱其中(n＝0,1,2,…,N-1)；
步骤3、确定控制因子a、b、c的值，得到窗长控制函数后确定窗函数，具体步骤为：
步骤3-1、确定时窗取值范围[Δtmin,Δtmax]，其中Δtmin为最小时窗长度，Δtmax为最大时窗长度，通过下列不等式确定a和c的取值范围：
步骤3-2、在取值范围内确定a和c的值，并取其中fs为采样频率，并将a、b、c的值代入窗长控制函数中：
步骤3-3、将窗长控制函数带入到窗函数中，得到改进后的窗函数表达式：
步骤4、对窗函数进行快速傅里叶变化得到窗函数频谱，对窗函数进行快速傅里叶变化得到窗函数频谱时，是以为起始频率点，具体公式为：
其中，n从0开始取值。
步骤5、将扩维后的信号频谱跟窗函数频谱相乘，再对其进行傅里叶逆变换，具体步骤为：
步骤5-1、将信号频谱扩维得到扩维后的信号频谱其中m＝0,1,2,…N-1；
步骤5-2、将扩维后的信号频谱与窗函数频谱G(m,n)相乘；
步骤5-3、对步骤5-2的结果进行傅里叶逆变换，得到
步骤6、重复步骤4、步骤5，直至所有频率点全部计算完成，得到高分辨率时频谱。判断所有频率点是否全部计算完成,具体是判断n是否满足大于等于N，若满足，输出时频谱；若不满足，n自加1后重复步骤4、步骤5；全部完成后判断时频谱的分辨率是否满足设定的分辨率要求，若满足，输出时频分析结果；若不满足，调整控制因子的值，重复步骤3、步骤4、步骤5，直至满足设定的分辨率要求，输出高分辨率时频谱。
本发明一方面实现了窗函数大小随着分析信号的频率进行自适应改变；另一方面，通过反正切函数控制时窗宽度在一定范围内随着频率的增大而减小，实现在时频谱的各个区域都有较高的分辨率。
下面以三种信号为例来说明本发明的时频分析方法。
实施例1
仿真信号为两个正弦信号的叠加，信号频率分别为100Hz和400Hz，解析式为：
h(t)＝sin(200πt)+sin(800πt) t∈[0,1]
信号采样频率fs＝1024Hz,图3为采用反正切函数联合S变换的有限窗长时频方法得到的时频谱。该信号只有两个固定的频率分量，对于频率不随时间变化的信号，只需考虑频率分辨率，因此通过调整控制因子，取a＝1,b＝50,c＝30将频窗宽度控制在较小范围内，保证时频谱有较好的频率分辨率性能。由图3可知，该方法可以实现很好的频率分辨率。
实施例2
仿真信号为调频斜率为k＝400的线性调频信号，解析式为：
信号采样频率fs＝1024Hz，图4为采用反正切函数联合S变换的有限窗长时频方法对线性调频信号的分析结果对于线性调频信号，频率变化较大，所以调整控制因子，取a＝5,b＝50,c＝50可以得到如图的分析结果，解决了原S变换高频处信号发散、能量聚集性差的问题，具有很好的时频性能。
实施例3
仿真信号为频率呈正弦变化的非线性调频信号，解析式为：
h(t)＝ej2π[6cos(10πt)+260t] t∈[0,1]
信号采样频率fs＝1024Hz，图5为非线性调频信号h(t)经过本发明方法处理后得到的时频图，因为该信号频率随时间呈正弦变化，所以对时间分辨率的要求相对较高，最小值控制因子c应较大，调整控制因子，从a＝6,b＝50,c＝100时候的时频结果中可以清楚地看到频率随时间的变化轨迹，具有很好的时频分析性能。

Claims (5)

1.一种改进S变换的有限窗长时频分析方法，其特征在于，包括以下步骤：
步骤1、对信号进行采样，采样点数为N＝t/T,其中t为信号时长，T为采样周期，得到信号的离散序列h[kT](k＝1,2,…,N)；
步骤2、对信号的离散序列h[kT]进行快速傅里叶变换，得到离散信号频谱其中(n＝0,1,2,…,N-1)；
步骤3、确定控制因子a、b、c的值，得到窗长控制函数后确定窗函数；
步骤4、对窗函数进行快速傅里叶变化得到窗函数频谱；
步骤5、将扩维后的信号频谱跟窗函数频谱相乘，再对其进行傅里叶逆变换；
步骤6、重复步骤4、步骤5，直至所有频率点全部计算完成，得到高分辨率时频谱。
2.根据权利要求1所述的改进S变换的有限窗长时频分析方法，其特征在于，步骤3中确定控制因子a、b、c的值，得到窗长控制函数后确定窗函数，具体步骤为：
步骤3-1、确定时窗取值范围[Δtmin,Δtmax]，其中Δtmin为最小时窗长度，Δtmax为最大时窗长度，通过下列不等式确定a和c的取值范围：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mi>&amp;pi;</mi> <mn>2</mn> </mfrac> <mi>a</mi> <mo>+</mo> <mi>c</mi> <mo>&gt;</mo> <mfrac> <mn>1</mn> <mrow> <msqrt> <mn>2</mn> </msqrt> <msub> <mi>t</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mi>&amp;pi;</mi> <mn>2</mn> </mfrac> <mi>a</mi> <mo>+</mo> <mi>c</mi> <mo>&lt;</mo> <mfrac> <mn>1</mn> <mrow> <msqrt> <mn>2</mn> </msqrt> <msub> <mi>t</mi> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
步骤3-2、在取值范围内确定a和c的值，并取其中fs为采样频率，并将a、b、c的值代入窗长控制函数中：
<mrow> <mi>w</mi> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>a</mi> <mo>&amp;lsqb;</mo> <mi>a</mi> <mi>r</mi> <mi>c</mi> <mi>t</mi> <mi>a</mi> <mi>n</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>f</mi> <mo>-</mo> <mfrac> <mrow> <mi>f</mi> <mi>m</mi> </mrow> <mn>2</mn> </mfrac> </mrow> <mi>b</mi> </mfrac> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>+</mo> <mi>c</mi> </mrow>
步骤3-3、将窗长控制函数带入到窗函数中，得到改进后的窗函数表达式：
<mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>w</mi> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <msqrt> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> </mrow> </msqrt> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mi>t</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>w</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
3.根据权利要求1所述的改进S变换的有限窗长时频分析方法，其特征在于，步骤4中对窗函数进行快速傅里叶变化得到窗函数频谱时，是以为起始频率点，具体公式为：
<mrow> <mi>G</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mn>2</mn> <msup> <mi>&amp;pi;</mi> <mn>2</mn> </msup> <msup> <mi>m</mi> <mn>2</mn> </msup> </mrow> <mrow> <msup> <mi>w</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mi>N</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中，n从0开始取值。
4.根据权利要求1所述的改进S变换的有限窗长时频分析方法，其特征在于，步骤5中将扩维后的信号频谱跟窗函数频谱相乘，再对其进行傅里叶逆变换，具体步骤为：
步骤5-1、将信号频谱扩维得到扩维后的信号频谱其中m＝0,1,2,…N-1；
步骤5-2、将扩维后的信号频谱与窗函数频谱G(m,n)相乘；
步骤5-3、对步骤5-2的结果进行傅里叶逆变换，得到
5.根据权利要求1所述的改进S变换的有限窗长时频分析方法，其特征在于，步骤6中判断所有频率点是否全部计算完成，具体是判断n是否满足大于等于N，若满足，输出时频谱；若不满足，n自加1后重复步骤4、步骤5；全部完成后判断时频谱的分辨率是否满足设定的分辨率要求，若满足，输出时频分析结果；若不满足，调整控制因子的值，重复步骤3、步骤4、步骤5，直至满足设定的分辨率要求，输出高分辨率时频谱。
CN201710593238.3A 2017-07-20 2017-07-20 一种改进s变换的有限窗长时频分析方法 Active CN107402326B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710593238.3A CN107402326B (zh) 2017-07-20 2017-07-20 一种改进s变换的有限窗长时频分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710593238.3A CN107402326B (zh) 2017-07-20 2017-07-20 一种改进s变换的有限窗长时频分析方法

Publications (2)

Publication Number Publication Date
CN107402326A true CN107402326A (zh) 2017-11-28
CN107402326B CN107402326B (zh) 2019-08-23

Family

ID=60401008

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710593238.3A Active CN107402326B (zh) 2017-07-20 2017-07-20 一种改进s变换的有限窗长时频分析方法

Country Status (1)

Country Link
CN (1) CN107402326B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108920418A (zh) * 2018-05-08 2018-11-30 电子科技大学 一种基于偏度的自适应变窗长短时时频变换技术
CN109270573A (zh) * 2018-09-14 2019-01-25 同济大学 一种快速保频保幅s变换方法
CN109343020A (zh) * 2018-11-16 2019-02-15 南京理工大学 一种基于改进窗函数的s变换时频分析方法
CN111198357A (zh) * 2019-12-19 2020-05-26 南京理工大学 一种基于可调窗函数的s变换时频分析方法
CN113092850A (zh) * 2021-04-16 2021-07-09 湖南师范大学 一种简化s变换的时频谱分析方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101975910A (zh) * 2010-09-07 2011-02-16 昆明理工大学 一种特高压直流输电线路故障智能分类与测距方法
CN104749432A (zh) * 2015-03-12 2015-07-01 西安电子科技大学 基于聚焦s变换的多分量非平稳信号瞬时频率估计方法
CN106405654A (zh) * 2016-10-26 2017-02-15 成都理工大学 一种基于反褶积广义s变换的地震频谱成像方法
CN106645947A (zh) * 2016-12-14 2017-05-10 南京航空航天大学 一种基于非线性模式分解和自适应最优核的时频分析方法
EP3168638A1 (en) * 2014-07-10 2017-05-17 Sfft Company Limited Peak frequency detection device, method, and program

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101975910A (zh) * 2010-09-07 2011-02-16 昆明理工大学 一种特高压直流输电线路故障智能分类与测距方法
EP3168638A1 (en) * 2014-07-10 2017-05-17 Sfft Company Limited Peak frequency detection device, method, and program
CN104749432A (zh) * 2015-03-12 2015-07-01 西安电子科技大学 基于聚焦s变换的多分量非平稳信号瞬时频率估计方法
CN106405654A (zh) * 2016-10-26 2017-02-15 成都理工大学 一种基于反褶积广义s变换的地震频谱成像方法
CN106645947A (zh) * 2016-12-14 2017-05-10 南京航空航天大学 一种基于非线性模式分解和自适应最优核的时频分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
朱明哲 等: "基于方向性S变换的多分量FM信号瞬时频率估计", 《系统工程与电子技术》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108920418A (zh) * 2018-05-08 2018-11-30 电子科技大学 一种基于偏度的自适应变窗长短时时频变换技术
CN108920418B (zh) * 2018-05-08 2021-06-01 电子科技大学 一种基于偏度的自适应变窗长短时时频变换方法
CN109270573A (zh) * 2018-09-14 2019-01-25 同济大学 一种快速保频保幅s变换方法
CN109270573B (zh) * 2018-09-14 2020-01-31 同济大学 一种快速保频保幅s变换方法
CN109343020A (zh) * 2018-11-16 2019-02-15 南京理工大学 一种基于改进窗函数的s变换时频分析方法
CN111198357A (zh) * 2019-12-19 2020-05-26 南京理工大学 一种基于可调窗函数的s变换时频分析方法
CN113092850A (zh) * 2021-04-16 2021-07-09 湖南师范大学 一种简化s变换的时频谱分析方法及系统

Also Published As

Publication number Publication date
CN107402326B (zh) 2019-08-23

Similar Documents

Publication Publication Date Title
CN107402326B (zh) 一种改进s变换的有限窗长时频分析方法
CN108009347A (zh) 基于同步压缩联合改进广义s变换的时频分析方法
Djurović et al. Frequency-based window width optimization for S-transform
CN109343020A (zh) 一种基于改进窗函数的s变换时频分析方法
CN103353550A (zh) 一种测量电力系统信号频率及谐波参数的方法
CN101701985B (zh) 定频变点电网谐波检测方法及其测量仪
CN106771567A (zh) 一种基于多分辨率短时傅里叶变换的动态谐波电能计量方法
CN104459398A (zh) 一种采用二维形态学降噪的电能质量复合扰动识别方法
CN105675956A (zh) 一种基于加窗插值短时傅里叶变换的电压闪变检测方法
CN103399204A (zh) 一种基于Rife-Vincent（II）窗插值FFT的谐波与间谐波检测方法
CN104714075A (zh) 一种电网电压闪变包络参数提取方法
CN105445541A (zh) 一种任意频率下自适应功率计算方法
CN105373708A (zh) 一种基于参数优化的改进广义s变换的时频分析方法
CN102724155B (zh) 基于分数傅里叶变换的高频域能量集中度同步方法
CN104880697B (zh) 基于稀疏约束的线性调频信号参数估计方法
Hague Nonlinear frequency modulation using fourier sine series
CN111198357A (zh) 一种基于可调窗函数的s变换时频分析方法
CN108761202A (zh) 极点对称模态分解和希尔伯特变换相结合的谐波检测方法
CN104914413A (zh) 一种随机序列线形调频信号加窗脉冲压缩方法
CN104008292A (zh) 宽带天线超宽带电磁脉冲响应快速准确预测方法
Razzaq et al. Instantaneous Frequency Estimation for Frequency-Modulated Signals under Gaussian and Symmetric α-Stable Noise
Xu et al. Harmonic estimation base on center frequency shift algorithm
Guo Wavelet packet transform‐based time of arrival estimation method for orthogonal frequency division multiplexing ultra‐wideband signal
CN104849551A (zh) 一种谐相角分析方法
CN111522000A (zh) 一种基于OFDM-chirp波形的目标检测方法

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