CN107402326B - 一种改进s变换的有限窗长时频分析方法 - Google Patents
一种改进s变换的有限窗长时频分析方法 Download PDFInfo
- Publication number
- CN107402326B CN107402326B CN201710593238.3A CN201710593238A CN107402326B CN 107402326 B CN107402326 B CN 107402326B CN 201710593238 A CN201710593238 A CN 201710593238A CN 107402326 B CN107402326 B CN 107402326B
- Authority
- CN
- China
- Prior art keywords
- window
- frequency
- time
- signal
- spectrum
- 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.)
- Expired - Fee Related
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 34
- 238000001228 spectrum Methods 0.000 claims abstract description 44
- 238000005070 sampling Methods 0.000 claims abstract description 10
- 230000003190 augmentative effect Effects 0.000 claims abstract description 5
- 238000000034 method Methods 0.000 abstract description 11
- 230000003044 adaptive effect Effects 0.000 abstract description 3
- 230000008859 change Effects 0.000 abstract description 2
- 238000013507 mapping Methods 0.000 abstract 1
- 238000000354 decomposition reaction Methods 0.000 description 4
- 230000009466 transformation Effects 0.000 description 4
- 238000009826 distribution Methods 0.000 description 3
- 230000007812 deficiency Effects 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 238000000701 chemical imaging Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
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变换的有限窗长高精度时频分析方法。
背景技术
非平稳信号是雷达信号处理中最常见的信号,而时频分析是分析该类信号的重要工具。为了准确分析信号的局部特性,时频分析将一维时域信号映射到二维时频平面,从而获得信号的时频分布。目前,常用的时频分析方法主要包括短时傅里叶变换(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 (4)
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的值,得到窗长控制函数后确定窗函数,具体步骤为:
步骤3-1、确定时窗取值范围[Δtmin,Δtmax],其中Δtmin为最小时窗长度,Δtmax为最大时窗长度,通过下列不等式确定a和c的取值范围:
步骤3-2、在取值范围内确定a和c的值,并取其中fs为采样频率,并将a、b、c的值代入窗长控制函数中:
步骤3-3、将窗长控制函数带入到窗函数中,得到改进后的窗函数表达式:
步骤4、对窗函数进行快速傅里叶变化得到窗函数频谱;
步骤5、将扩维后的信号频谱跟窗函数频谱相乘,再对其进行傅里叶逆变换;
步骤6、重复步骤4、步骤5,直至所有频率点全部计算完成,得到高分辨率时频谱。
2.根据权利要求1所述的改进S变换的有限窗长时频分析方法,其特征在于,步骤4中对窗函数进行快速傅里叶变化得到窗函数频谱时,是以为起始频率点,具体公式为:
其中,n从0开始取值。
3.根据权利要求1所述的改进S变换的有限窗长时频分析方法,其特征在于,步骤5中将扩维后的信号频谱跟窗函数频谱相乘,再对其进行傅里叶逆变换,具体步骤为:
步骤5-1、将信号频谱扩维得到扩维后的信号频谱其中m=0,1,2,…N-1;
步骤5-2、将扩维后的信号频谱与窗函数频谱G(m,n)相乘;
步骤5-3、对步骤5-2的结果进行傅里叶逆变换,得到
4.根据权利要求1所述的改进S变换的有限窗长时频分析方法,其特征在于,步骤6中判断所有频率点是否全部计算完成,具体是判断n是否满足大于等于N,若满足,输出时频谱;若不满足,n自加1后重复步骤4、步骤5;全部完成后判断时频谱的分辨率是否满足设定的分辨率要求,若满足,输出时频分析结果;若不满足,调整控制因子的值,重复步骤3、步骤4、步骤5,直至满足设定的分辨率要求,输出高分辨率时频谱。
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 CN107402326A (zh) | 2017-11-28 |
CN107402326B true CN107402326B (zh) | 2019-08-23 |
Family
ID=60401008
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710593238.3A Expired - Fee Related CN107402326B (zh) | 2017-07-20 | 2017-07-20 | 一种改进s变换的有限窗长时频分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107402326B (zh) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108920418B (zh) * | 2018-05-08 | 2021-06-01 | 电子科技大学 | 一种基于偏度的自适应变窗长短时时频变换方法 |
CN109270573B (zh) * | 2018-09-14 | 2020-01-31 | 同济大学 | 一种快速保频保幅s变换方法 |
CN109343020B (zh) * | 2018-11-16 | 2023-04-07 | 南京理工大学 | 一种基于改进窗函数的s变换时频分析方法 |
CN111198357A (zh) * | 2019-12-19 | 2020-05-26 | 南京理工大学 | 一种基于可调窗函数的s变换时频分析方法 |
CN113092850B (zh) * | 2021-04-16 | 2022-03-15 | 湖南师范大学 | 一种简化s变换的时频谱分析方法及系统 |
Citations (5)
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 |
-
2017
- 2017-07-20 CN CN201710593238.3A patent/CN107402326B/zh not_active Expired - Fee Related
Patent Citations (5)
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)
Title |
---|
基于方向性S变换的多分量FM信号瞬时频率估计;朱明哲 等;《系统工程与电子技术》;20131031;第35卷(第1期);第29-34页 |
Also Published As
Publication number | Publication date |
---|---|
CN107402326A (zh) | 2017-11-28 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107402326B (zh) | 一种改进s变换的有限窗长时频分析方法 | |
CN108009347A (zh) | 基于同步压缩联合改进广义s变换的时频分析方法 | |
CN109343020A (zh) | 一种基于改进窗函数的s变换时频分析方法 | |
CN104749432B (zh) | 基于聚焦s变换的多分量非平稳信号瞬时频率估计方法 | |
CN110046323B (zh) | 同步压缩变换与重构的快速计算方法 | |
CN111198357A (zh) | 一种基于可调窗函数的s变换时频分析方法 | |
CN104182762A (zh) | 一种基于pir探测器的人与非人识别方法 | |
Zhang et al. | Micro-motion frequency estimation of radar targets with complicated translations | |
CN107248869B (zh) | 一种基于吕分布的多分量线性调频信号去噪方法 | |
Wang et al. | Approach for high‐resolution inverse synthetic aperture radar imaging of ship target with complex motion | |
Jiang et al. | Instantaneous frequency estimation of nonlinear frequency-modulated signals under strong noise environment | |
Jang et al. | Exploiting early time response using the fractional Fourier transform for analyzing transient radar returns | |
Abratkiewicz | Double‐adaptive chirplet transform for radar signature extraction | |
Li et al. | FFBV algorithm design for soft synchronization FMCW radar | |
CN104914413A (zh) | 一种随机序列线形调频信号加窗脉冲压缩方法 | |
Wang et al. | A sparse reconstruction algorithm based on constrained inhomogeneous grid optimization | |
Stanković et al. | A real-time time-frequency based instantaneous frequency estimator | |
Ba et al. | Maximum likelihood time delay estimation based on Monte Carlo importance sampling in multipath environment | |
Guo | Wavelet packet transform‐based time of arrival estimation method for orthogonal frequency division multiplexing ultra‐wideband signal | |
Liu et al. | An efficient algorithm for frequency estimation of sinusoid signal based on improved quinn | |
He et al. | Blind estimation for PN code of LFM-PRBC signal based on DPT and spectrum shifting | |
Bai et al. | Robust speech recognition model using multi-source federal learning after distillation and deep edge intelligence | |
CN111522000A (zh) | 一种基于OFDM-chirp波形的目标检测方法 | |
Yang et al. | A parameter estimation algorithm for frequency-hopping signals with a stable noise | |
Tian et al. | Parameter estimation of manoeuvring targets based on segment integration and Lv's transform |
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 | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190823 |
|
CF01 | Termination of patent right due to non-payment of annual fee |