CN104007318A - 获取信号时频函数的方法 - Google Patents

Abstract

本发明提供了一种获取信号时频函数的方法。该方法将分数阶傅里叶变换和WVD变换结合起来，其中，分数阶傅里叶变换能够克服传统STFT的不足，分数阶傅里叶变换尤其对线性调频信号具有很好的时频聚集性，因此在对非平稳信号进行时频分析时具有特别的优势，而且分数阶傅里叶变换不存在交叉项干扰。基于上述理由，本发明具有计算精度高、抗噪性能强等优点，可以用于地震、声纳、振动检测等许多领域。

Description

获取信号时频函数的方法

技术领域

本发明涉及电子行业信号处理技术领域，尤其涉及一种获取信号时频函数的方法。

背景技术

现代信号处理中的信号一般具有时间参数和频率参数。对于非平稳信号即信号的频率可能会随时间的改变而变化，此时，时间和频率参数转换为时间频率函数，简称时频函数。而获取信号时频函数的方法，本领域内称之为时频分析。对雷达信号计算出信号的时间频率函数，可以用来检查信号的线性度、分离信号的成分、计算出信号的多普勒频率。非平稳信号的时频分析还被广泛应用于地震、声纳、振动检测等许多领域。因此对于信号来说，时间和频率的函数非常重要，有必要对信号做时频分析。

信号处理领域使用最频繁的傅里叶变换可以计算出信号的频率，但是它是一种整体性的从时间域到频率域的变换，不具有时间分辨率，基于此缺陷，Gabor于1946年引入了短时傅里叶变换(STFT)，其基本思想是利用窗函数来截取一段信号，对截取的信号进行傅里叶变换，由于所截取的信号时间较短，因此可以作为平稳信号处理，通过窗函数的不断移位，便可以计算出各个时刻的信号频率，STFT没有交叉项的干扰。分数阶傅里叶变换尤其对线性调频信号具有很好的时频聚集性，因此在对非平稳信号进行时频分析时具有特别的优势，而且分数阶傅里叶变换不存在交叉项干扰，分数阶傅里叶变换在时频分析上有其独到的优势。

然而，在利用现有的STFT进行信号时频分析时，虽然没有交叉项的干扰，但是其自项不是十分集中，其傅里叶变换的频谱带宽较宽，STFT不能精确计算出其对应时刻的频率，因此会产生较大的计算误差。此外，还有一些时频分析方法其算法精度有限，不能满足一些高精度的要求。

发明内容

(一)要解决的技术问题

鉴于上述技术问题，本发明提供了一种高精度的获取信号时频函数的方法。

(二)技术方案

本发明获取信号时频函数的方法包括：步骤A，对待处理信号s(t)加窗长为T_w的窗函数w(t)，得到分段信号s_i(t)，其中，i＝1、2、3、……；步骤B，针对每个分段信号，计算该分段信号s_i(t)的调频斜率K′_i；步骤C，针对每个分段信号，依据预设调频斜率k，构建线性调频信号模型s_ci(t)，获取其时频函数f_ci(t)；步骤D，针对每个分段信号，利用分段信号s_i(t)与线性调频信号模型s_ci(t)混频得到调频信号x_i(t)；步骤E：针对每个分段信号，对于混频得到的调频信号x_i(t)进行离散的WVD变换，得到该调频信号x_i(t)的离散的二维频谱W_xi(t，f)和跳变点处频率F_q，其中，q＝1、2、……、Q，Q为跳变点的个数；步骤F：针对每个分段信号，利用窗长T_w内的Q个的跳变点处的频率进行线性插值计算并将插值计算结果除以2，得到分段信号的时频函数f′_i(t)，t∈[0，T_w]；步骤G：针对每个分段信号，将时频函数f′_i(t)减去线性调频信号模型s_ci(t)的时频函数f_ci(t)，得到该分段信号的时频函数F′_i(t)；以及步骤H：将各段的分段信号的时频函数按分段顺序组合起来得到整段待处理信号s(t)的时频函数。

(三)有益效果

从上述技术方案可以看出，本发明获取信号时频函数具有以下有益效果：

(1)本发明在原有采样N点时频率分辨率为其中Fs为采样率，N为采样点数，经过本发明后，时频分析的精度可以提高几十倍，具体根据***的需求和信号的质量来决定。

(2)分数阶傅里叶变换能够克服传统STFT的不足，分数阶傅里叶变换尤其对线性调频信号具有很好的时频聚集性，因此在对非平稳信号进行时频分析时具有特别的优势，而且分数阶傅里叶变换不存在交叉项干扰，分数阶傅里叶变换在时频分析上有其独到的优势。

附图说明

图1为根据本发明实施例获取信号时频函数方法的流程图；

图2为分数阶傅里叶变坐标旋转换示意图；

图3为调频斜率较大时的时频图；

图4为调频斜率较小时的时频图；

图5为调频斜率合适时的时频图；

图6为栅栏效应示意图；

图7为栅栏效应整数倍频率分辨率频谱图；

图8为栅栏效应频谱最大值和次大值相当时频谱图；

图9为MATLAB仿真本方法时频图和误差图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向。因此，使用的方向用语是用来说明并非用来限制本发明的保护范围。

本发明获取信号时频函数的方法是基于分数阶傅里叶变换和WVD变换提出的一种新型的时频分析方法。该方法通过仿真实验和理论验证，两者的结果均证明了本发明的有效性。

在本发明的一个示例性实施例中，提供了一种获取信号时频函数的方法。图1为根据本发明实施例获取信号时频函数方法的流程图。如图1所示，本实施例获取信号时频函数的方法包括：

步骤A，对待处理信号s(t)加窗长为T_w的窗函数w(t)，得到分段信号s_i(t)，其中，i＝1、2、3、……；

待处理的信号可以表示为：

$s\left(t\right)=A\·e^{j2\mathrm{\π}(f_{0}t+\frac{k_2{t}^2}{2}+\frac{k_3{t}^3}{3}...+\frac{k_n{t}^n}{n})}---\left(1\right)$

其中，t为时间，且0＜t＜T，T为信号时长，A为信号幅度，f₀为信号初始频率，k₂、k₃、…k_n分别为信号频率的二次项、三次项、n次项系数，e为自然对数底数，j为虚数单位。

需要说明的是，该待处理信号s(t)的采样率为Fs，采样点数为N，该采样率Fs和采样点数N在后续步骤中会采用。一般情况下，窗长T_w可以根据***需求和计算精度，一般取值为信号时长的1/50到1/100之间。

对信号s(t)加窗长为T_w的窗函数，对信号s(t)加第i个窗时，窗向右移动因此，加第i个窗函数时间内的信号为：

$s_i\left(t\right)=[A\·e^{j2\mathrm{\π}(f_{0}t+\frac{k_2{t}^2}{2}+\frac{k_3{t}^3}{3}...)}]\·w(t-(2i-1)\·\frac{T_w}{2})---\left(2\right)$

其中，w(t)为窗函数，可选的窗函数有矩形窗、汉宁窗、汉明窗等。

步骤B，针对每个分段信号，计算分段信号s_i(t)的调频斜率K′_i；

把每段信号看作时间从0到T_w，可以将s_i(t)看做一近似线性调频信号，由于时窗函数的窗长T_w较短，对应s_i(t)中的幂级数大于3的高次项求导后所对应的频率变化较小，对计算的影响不大，可以忽略，只剩下一次项和二次项，因此可以将信号s_i(t)独立的表示为一线性调频信号：

$s_i\left(t\right)=[A\·e^{j2\mathrm{\π}(f_i\·t+\frac{K_i\·t^2}{2})}]\·w(t-\frac{T_w}{2})---\left(3\right)$

其中，fi为该分段信号的起始频率，K_i为该分段信号的调频斜率。

图2为分数阶傅里叶变换坐标旋转的示意图。利用分数阶傅里叶变换计算加窗分段信号s_i(t)各段的调频斜率K′_i的步骤如下：

子步骤B1：对待处理信号s_i(t)进行分数阶傅里叶变换，得到S_p(u)：

$S_p\left(u\right)=F^p\left[s_i\left(t\right)\right]={\∫}_{-\∞}^{+\∞}{s}_i\left(t\right)K_p(t,u)\mathrm{dt}---\left(4\right)$

式中，p为分数阶傅里叶变换的阶，t为时间，u为分数阶域的自变量，变换核为K_p(t，u)，且

$K_p(t,u)=\left\{\begin{array}{c}\sqrt{\frac{1-j\cot \mathrm{\α}}{2\mathrm{\π}}}\·e^{(j\frac{t^2+u^2}{2}\cot \mathrm{\α}-\mathrm{jtu}\csc \mathrm{\α})}...\mathrm{\α}\≠\mathrm{n\π}\\ \mathrm{\δ}(t-u)...\mathrm{\α}=2\mathrm{n\π}\\ \mathrm{\δ}(t+u)...\mathrm{\α}=(2n\±1)\mathrm{\π}\end{array}\right.---\left(5\right)$

其中，δ为单位冲击脉冲函数，j为虚数单位，n为整数，α为旋转角度且cot()和csc()为三角函数。

子步骤B2：计算S_p(u)的绝对值|S_p(u)|，并找到|S_p(u)|中最大值所对应的旋转角度α_i；

本实施例中，寻找|S_p(u)|中最大值所对应的旋转角度α_i的方法为本领域技术人员所熟知，此处不再详细描述。

子步骤B3，根据计算出调频斜率K′_i，其中，Fs为待处理信号s(t)的采样率，N为待处理信号s(t)的采样点数。

需要说明的是，利用分数阶傅里叶变换计算加窗分段信号s_i(t)各段的调频斜率K′_i，由于计算存在误差和精度限制，所计算的调频斜率K′_i和信号理论的调频斜率K_i不完全相等，存在一定的误差，可以满足近似相等，不影响后续计算。

步骤C，针对每个分段信号，依据预设调频斜率k，构建线性调频信号模型s_ci(t)，获取其时频函数f_ci(t)；

根据***和算法精度要求选择合适的预设调频斜率k可以选择实际***的需求来选择，当k取较大值时，在窗长T_w内阶梯跳跃次数较多，便于后续插值处理但是会导致算法的精度下降，如图3所示；当k取较小值时，在时间T_w内阶梯跳跃次数较少，可以有效提高算法的精度但是会导致跳变点较少而不利于后续的插值处理，如图4所示。一般以在窗长T_w内有十个跳变点左右为宜，因此k值一般取如图5所示。

线性调频信号模型的调频斜率为：k-K′_i，其中k为所选择的合适的预设调频斜率，K′_i是利用分数阶傅里叶变换所计算的分段信号s_i(t)的调频斜率。因此，构建线性调频信号模型的表达式为：

$s_{\mathrm{ci}}\left(t\right)=e^{j2\mathrm{\π}}\frac{(k-K_i^{\′})}{2}{t}^2---\left(6\right)$

其时频函数为：f_ci(t)＝(k-K′_i)·t，t∈[0，T_w]，其中，j为虚数单位，t为时间，e为自然对数底数。

步骤D，针对每个分段信号，利用分段信号s_i(t)与线性调频信号模型s_ci(t)混频得到调频信号x_i(t)；

依照下式，利用线性调频信号模型s_ci(t)与分段的信号s_i(t)混频得到调频信号：

$x_i\left(t\right)=s_{\mathrm{ci}}\left(t\right)\·s_i\left(t\right)=A\·e^{j2\mathrm{\π}[f_i\·t+\frac{(k-K_i+K_i^{\′})}{2}\·t^2]}---\left(7\right)$

其时频函数为：f_xi(t)＝f_i+(k-K_i+K′_i)·t，t∈[0，T_w]，其中，j为虚数单位，t为时间，e为自然对数底数，K_i是分段信号s_i(t)的理论调频斜率，因K′_i-K_i≈0，因此混频后的信号：

$x_i\left(t\right)=s_{\mathrm{ci}}\left(t\right)\·s_i\left(t\right)\≈A\·e^{j2\mathrm{\π}(f_i\·t+\frac{k}{2}\·t^2)}---\left(8\right)$

其时频函数约等于：f′_xi(t)＝f_i+k·t，t∈[0，T_w]。

步骤E：针对每个分段信号，对于混频得到的调频信号x_i(t)进行离散的WVD变换，得到该信号x_i(t)的离散的二维频谱W_xi(t，f)和跳变点处频率F_q，其中，q＝1、2、……、Q，Q为跳变点的个数；

信号x_i(t)的WVD定义为：

$W_{\mathrm{xi}}(t,f)={\∫}_{-\∞}^{+\∞}{x}_i(t+\frac{\mathrm{\τ}}{2})\·{x_i}^{*}(t-\frac{\mathrm{\τ}}{2})\·e^{j2\mathrm{\πf\τ}}\mathrm{d\τ}---\left(9\right)$

其中x_i(t)为信号的解析信号，令上式中积分部分等于r_x，x(t，τ)：

$r_{x,x}(t,\mathrm{\τ})=x_i(t+\frac{\mathrm{\τ}}{2})\·{x_i}^{*}(t-\frac{\mathrm{\τ}}{2})---\left(10\right)$

上式被称作信号的瞬时自相关，因此WVD变换可以看作为瞬时自相关的傅里叶变换。瞬时自相关会导致信号瞬时频率增加一倍，因此在计算信号的时频函数时所计算的信号频率应该除以2。由于对信号的离散化会产生栅栏效应，栅栏效应是由于离散采样造成的频谱离散化。

栅栏效应是指当信号频率不为频率分辨率的整数倍时，中频信号的离散频谱为连续频率的抽样，并且离散频谱最大值与连续频谱最大值不相等，因此离散频谱最大值所对应的频率F_m不等于信号的频率值F，有一定的误差，如图6所示；只有当信号频率为频率分辨率的整数倍时，离散频谱最大值与连续频谱最大值重合，此时离散频谱最大值所对应的频率F_m才是信号的频率值F，如图7所示。

当离散频谱最大值与次大值幅度差不多时，其信号的频率F近似等于离散频谱最大值所对应的频率F_m和离散频谱次大值所对应的频率F_c二者和的二分之一：如图8所示，并且由于离散频谱最大值和离散频谱次大值相差一条谱线，因此其频率差的绝对值等于频率分辨率：其中Fs为离散采样率，N为离散采样点数。

在x_i(t)的离散WVD变换后，由于x_i(t)为线性调频信号，并且其预设调频斜率k可以保证在窗长T_w内其频率的变化等于十个频率分辨率左右，如图5所示，因此可以保证其离散频谱发生渐变，使其离散频谱最大值与次大值出现十次左右的跳变，在这些跳变点处离散频谱最大值与离散频谱次大值幅度差不多，在这些跳变点处频率F_q近似等于二维频谱W_xi(t，f)中相邻离散频谱最大值所对应的频率F_mq和相邻离散频谱次大值所对应的频率F_cq二者和的二分之一：

$F_q=\frac{F_{\mathrm{cq}}+F_{\mathrm{mq}}}{2}---\left(11\right)$

其中，q＝1、2、……、Q，Q为跳变点的个数。Q的大小与预设的调皮斜率k有关，一般情况下，8≤Q≤12。

步骤F：针对每个分段信号，利用窗长T_w内的Q个的跳变点处的频率进行线性插值计算，得到窗长T_w内的时频函数f″_i(t)，t∈[0，T_w]，由于WVD变换引入瞬时自相关使信号频率增加一倍，因此所计算的频率函数f″_i(t)需要除以2得到分段信号的时频函数f′_i(t)。

步骤G：针对每个分段信号，将时频函数f′_i(t)减去线性调频信号模型s_ci(t)的时频函数f_ci(t)，得到该分段信号的时频函数F′_i(t)：

F′_i(t)＝f′_i(t)-f_ci(t)，t∈[0，T_w]    (12)

步骤H：将各段的分段信号的时频函数按分段顺序组合起来得到整段待处理信号s(t)的时频函数。

通过之前的步骤可以计算出分段信号各段的时频函数，每段时频函数窗长为T_w，将这些时频函数按分段顺序组合起来得到整段信号的时频函数。

为了验证本实施例获取信号时频函数方法的有效性，采用MATLAB软件对算法进行仿真实现。设信号的频率为一正余弦变换的信号，其中心频率为600Hz，最大频率变化为500Hz，该信号的表达式为：

$s\left(t\right)=\cos \left(2\mathrm{\π}(600t-500\·\frac{T}{2\mathrm{\π}}\·\cos \left(2\mathrm{\π}\frac{t}{T}\right))\right)+A_n\·N\left(t\right)---\left(13\right)$

其中，T为信号总的时间长度，仿真中T＝5s，N(t)为高斯白噪声，A_n为白噪声系数，调整A_n的大小可以得到不同信噪比的信号，因此对信号s(t)的相位求导得到其频率的表达式为：

$F\left(t\right)=600+500\sin \left(2\mathrm{\π}\frac{t}{T}\right)---\left(14\right)$

以采样率Fs＝20480Hz对信号进行离散化处理，构建时间长度为的短时窗函数，利用短时窗函数将信号分为50段。在信噪比为5dB时，本文方法对信号s(t)的频率估计时频图和误差图如图9所示。从误差分布图中可以看出本方法时频分析精度高，在信噪比为5dB时的时频法分析误差在0.5Hz以内。如果用STFT做时频分析的精度为频率分辨率因此在低信噪比下比原来的时频分析算法提高了20倍。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明获取信号时频函数的方法有了清楚的认识。

此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。

综上所述，本发明将分数阶傅里叶变换和WVD变换结合起来，提供一种获取信号时频函数的方法，具有计算精度高、抗噪性能强等优点，可以用于地震、声纳、振动检测等许多领域。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种获取信号时频函数的方法，其特征在于，包括：

步骤A，对待处理信号s(t)加窗长为T_w的窗函数w(t)，得到分段信号s_i(t)，其中，i＝1、2、3、……；

步骤B，针对每个分段信号，计算该分段信号s_i(t)的调频斜率K′_i；

步骤C，针对每个分段信号，依据预设调频斜率k，构建线性调频信号模型s_ci(t)，获取其时频函数f_ci(t)；

步骤D，针对每个分段信号，利用分段信号s_i(t)与线性调频信号模型s_ci(t)混频得到调频信号x_i(t)；

步骤E：针对每个分段信号，对于混频得到的调频信号x_i(t)进行离散的WVD变换，得到该调频信号x_i(t)的离散的二维频谱W_xi(t，f)和跳变点处频率F_q，其中，q＝1、2、……、Q，Q为跳变点的个数；

步骤F：针对每个分段信号，利用窗长T_w内的Q个的跳变点处的频率进行线性插值计算并将插值计算结果除以2，得到分段信号的时频函数f′_i(t)，t∈[0，T_w]；

步骤G：针对每个分段信号，将时频函数f′_i(t)减去线性调频信号模型s_ci(t)的时频函数f_ci(t)，得到该分段信号的时频函数F′_i(t)；以及

步骤H：将各段的分段信号的时频函数按分段顺序组合起来得到整段待处理信号s(t)的时频函数。

2.根据权利要求1所述的方法，其特征在于，所述步骤B中计算分段信号s_i(t)的调频斜率K′_i具体包括：

子步骤B1：对待处理信号s_i(t)进行分数阶傅里叶变换，得到S_p(u)：

子步骤B2：计算S_p(u)的绝对值|S_p(u)|，并找到|S_p(u)|中最大值所对应的旋转角度α_i；以及

子步骤B3，根据计算出调频斜率K′_i，其中，Fs为待处理信号s(t)的采样率，N为待处理信号s(t)的采样点数。

3.根据权利要求2所述的方法，其特征在于，所述子步骤B1中：

$S_p\left(u\right)=F^p\left[s_i\left(t\right)\right]={\∫}_{-\∞}^{+\∞}{s}_i\left(t\right)K_p(t,u)\mathrm{dt}$

式中，p为分数阶傅里叶变换的阶，u为分数阶域的自变量，变换核为K_p(t，u)，且

$K_p(t,u)=\left\{\begin{array}{c}\sqrt{\frac{1-j\cot \mathrm{\α}}{2\mathrm{\π}}}\·e^{(j\frac{t^2+u^2}{2}\cot \mathrm{\α}-\mathrm{jtu}\csc \mathrm{\α})}...\mathrm{\α}\≠\mathrm{n\π}\\ \mathrm{\δ}(t-u)...\mathrm{\α}=2\mathrm{n\π}\\ \mathrm{\δ}(t+u)...\mathrm{\α}=(2n\±1)\mathrm{\π}\end{array}\right.$

其中，δ为单位冲击脉冲函数，n为整数，α为旋转角度且

4.根据权利要求1所述的方法，其特征在于，所述步骤C中：

$s_{\mathrm{ci}}\left(t\right)=e^{j2\mathrm{\π}}\frac{(k-K_i^{\′})}{2}{t}^2$   f_ci(t)＝(k-K′_i)·t。

5.根据权利要求1所述的方法，所述步骤C中，所述根据***和计算精度预设调频斜率k，k的取值一般可以取

6.根据权利要求1所述的方法，所述步骤D中，混频得到调频信号x_i(t)为：

$x_i\left(t\right)=s_{\mathrm{ci}}\left(t\right)\·s_i\left(t\right)\≈A\·e^{j2\mathrm{\π}(f_i\·t+\frac{k}{2}\·t^2)}$

其中，该调频信号x_i(t)的时频函数等于：f_xi(t)＝f_i+k·t，t∈[0，T_w]。

7.根据权利要求1所述的方法，所述步骤E中：

$F_q=\frac{F_{\mathrm{cq}}+F_{\mathrm{mq}}}{2}$

其中，F_mq和F_cq分别二维频谱W_xi(t，f)中第q个跳变点相邻离散频谱最大值和次大值所对应的频率。

8.根据权利要求1所述的方法，所述步骤E中一般情况下：8≤Q≤12。

9.根据权利要求1至8中任一项所述的方法，其特征在于，所述步骤A中窗函数为矩形窗、汉宁窗或汉明窗，所述窗长T_w介于信号时长T的1/50～1/100之间。

10.根据权利要求1至8中任一项所述的方法，其特征在于，所述步骤A中，待处理的信号为：

$s\left(t\right)=A\·e^{j2\mathrm{\π}(f_{0}t+\frac{k_2{t}^2}{2}+\frac{k_3{t}^3}{3}...+\frac{k_n{t}^n}{n})}$

其中，t为时间，且0＜t＜T，T为信号时长，A为信号幅度，f₀为信号初始频率，k₂、k₃、…k_n分别为信号频率的二次项、三次项、n次项系数。