CN106134472B - 基于微分的am-fm信号瞬时频率计算方法 - Google Patents

Abstract

基于微分的AM-FM信号瞬时频率计算方法，步骤为：(1)输入AM-FM信号x(t)＝a(t)cos[θ(t)]，并进行去噪及离散采样处理；(2)求取信号x(t)＝a(t)cos[θ(t)]的局部极大、极小值及发生时刻，记为(t_j，a(t_j))；(3)根据得到的各不同时刻的a(t_j)值，采用样条插值函数拟合a(t_j)的绝对值，得到输入信号的包络函数a(t)；(4)求取输入信号x(t)＝a(t)cos[θ(t)]的一阶导数x′(t)；(5)求取x′(t)的局部极大、极小值及发生时刻，记为(t_i，x′(t_i))；(6)根据公式计算t_i时刻输入信号的瞬时频率f(t_i)；(7)对x(t)＝a(t)cos[θ(t)]进行移相处理，分别得到x₁(t)，x₂(t)，…x_n-1(t)，其中n≥1；(8)对得到的x₁(t)，x₂(t)，…x_n-1(t)，重复步骤(2)～(6)，得到(t_i，f(t_i))点序列；(9)根据得到的(t_i，f(t_i))点序列，采用样条插值函数拟合f(t_i)，得到x(t)＝a(t)cos[θ(t)]的瞬时频率函数f(t)。

Description

基于微分的AM-FM信号瞬时频率计算方法

技术领域

本发明涉及一种瞬时频率计算方法，特别是一种基于微分的AM-FM信号瞬时频率计算方法。

背景技术

非平稳信号是统计量(相关函数、功率谱等)随时间变化的信号，早期对非平稳信号几乎没有什么正规的分析方法，发展最早、最成熟的信号分析方法是傅立叶变换。傅立叶变换理论中表征信号交变的基本量是与时间无关的频率，基本时域信号是平稳的简谐波信号，这些概念是全局性的，用其分析平稳信号很有效，但用它们分析非平稳信号容易产生虚假信号和假频等矛盾现象，并且用于分析非平稳信号时缺乏物理意义。因此时频分析方法是分析非平稳信号的有效手段。对非平稳信号比较直观的分析方法是使用具有局域性的基本量和基本函数。瞬时频率是容易想到的具有局域性的基本量，也是很早就已提出的概念。目前与本发明方法有关的求取瞬时频率的方法主要有相位差分法、过零点法以及HHT法，过零点法实际上等价于前向相位差分方法。差分法首先要对实信号x(t)进行一次Hilbert变换，构成一个解析信号z(t)＝x(t)+jH[x(t)]，求得z(t)的相位函数然后再应用差分公式对φ(t)进行计算。该方法的主要不足之处在于构造解析信号计算量大，瞬时频率结果误差大。HHT法是美国的N.E.Huang于1998年在美国专利号：US5983162的专利中公开的一种非线性、非平稳信号处理方法。HHT方法首先通过用经验模态分解方法(Empirical Mode Decomposition，EMD)将信号分解成一系列内蕴模态函数(Intrinsic Mode Function，IMF)，通过对每一个IMF构造其解析信号，再通过相位的一阶导数求出瞬时频率，从而达到获取瞬时频率、包络的目的。该方法同样是对IMF构造解析信号求取瞬时频率，不足之处也在于构造解析信号计算量大，瞬时频率结果误差大。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种具有自适应性、计算量小、易于通过软、硬件实现的基于微分的AM-FM信号瞬时频率计算方法。

本发明的技术解决方案是：基于微分的AM-FM信号瞬时频率计算方法，其特征在于步骤如下：

(1)输入AM-FM信号x(t)＝a(t)cos[θ(t)]，其中a(t)为输入信号的幅值，θ(t)为输入信号的相位；

(2)求取输入信号x(t)＝a(t)cos[θ(t)]的局部极大值和极小值及发生时刻，记为(t_j，a(t_j))；

(3)根据步骤(2)得到的各不同时刻的a(t_j)值，采用样条插值函数拟合所述a(t_j)的绝对值，得到输入信号的包络函数a(t)；

(4)求取输入信号x(t)＝a(t)cos[θ(t)]的一阶导数x′(t)；

(5)求取x′(t)的局部极大值和极小值及发生时刻，记为(t_i，x′(t_i))；

(6)根据公式计算t_i时刻输入信号的瞬时频率f(t_i)；

(7)对输入信号x(t)＝a(t)cos[θ(t)]进行移相处理，分别得到x₁(t)，x₂(t)，…x_n-1(t)，其中n≥1；

(8)对于步骤(7)得到的x₁(t)，x₂(t)，…x_n-1(t)，重复步骤(2)～(6)，得到(t_i，f(t_i))点序列；

(9)根据步骤(8)得到的(t_i，f(t_i))点序列，采用样条插值函数拟合所述的f(t_i)，得到输入信号x(t)＝a(t)cos[θ(t)]的瞬时频率函数f(t)。

所述的输入信号x(t)＝a(t)cos[θ(t)]应满足的要求为：

(1)a(t)＞0；

(2)θ(t)≥0，但不要求单调；

(3)a(t)的Fourier分析所得最高频率必须小于cosθ(t)的Fourier分析所得最低频率

(4)x′(t)·(cosθ(t))′≥0。

所述步骤(7)中对输入信号x(t)＝a(t)cos[θ(t)]进行移相处理的公式为：

$x_i\left(t\right)=x\left(t\right)cos\left(\frac{i\·\mathrm{\π}}{n}\right)+H\[x\left(t\right)\]sin\left(\frac{i\·\mathrm{\π}}{n}\right),i=1,2,...n-1,$

其中H[·]为Hilbert变换，n≥1。

本发明与现有技术相比的优点在于：

1)本发明的瞬时频率计算方法立足于信号本身进行处理，充分利用信号的自身信息，具有比较高的自适应性；

2)与HHT方法计算瞬时频率相比，本发明方法不必进行Hilbert变换构造解析信号，进而利用解析信号的相位微分求取瞬时频率，计算量小；

3)本发明利用AM-FM信号移相处理得到不同相位信号的过零点时刻的瞬时频率，提高了拟合准确度；

4)本发明方法易于通过软、硬件实现。

附图说明

图1为简谐振动s(t)＝r cosθ(t)＝r cosωt的位移示意图；

图2为一般振动s(t)＝a(t)cosθ(t)在x-y平面投影示意图；

图3为正交两路AM-FM信号示意图；

图4为n路AM-FM信号示意图；

图5为信号x(t)＝cos(πt²)不同相位零点-极值点关系示意图；

图6为本发明基于微分的AM-FM信号瞬时频率计算流程图；

图7为正弦信号x(t)＝sin(2π·t)及其瞬时频率示意图；

图8为调频信号x(t)＝sin(10πt²)及其瞬时频率示意图；

图9为指数调频信号x(t)＝cos[2π·exp(3t)]及其瞬时频率示意图；

图10为复合调频信号及其瞬时频率示意图；

图11为正弦信号x(t)＝sin(2π·t)通过本发明、HHT所求瞬时频率及真实频率的对比图；

图12为调频信号x(t)＝sin(10πt²)通过本发明、HHT所求瞬时频率及真实频率的对比图；

图13为指数调频信号x(t)＝cos[2π·exp(3t)]通过本发明、HHT所求瞬时频率及真实频率的对比图；

图14为复合信号通过本发明、HHT所求瞬时频率及真实频率的对比图；

图15为本发明实际应用中时空二维联合采样示意图。

具体实施方式

一、AM-FM信号的物理模型

在振动理论中，频率是表示一个交变信号的基本变量，它是指单位时间内运动物体往复振动的次数。如图1所示，当物体M以角速度ω沿半径为r的圆周运动时，物体M在直径上的投影P是一个简谐振动，其位移为：

s(t)＝r cosθ(t)＝r cosωt (1-1)

当时间t经过2π/ω后，P往复振动一次。显然，当角速度ω越大，P往复运动越快。但投影P仅做的是在x轴、[-r，r]范围内的往复运动，并不能形成r cosωt这样的轨迹。当M沿着垂直于振动平面、过坐标原点的直线同时也做速度为v(t)的直线运动时，物体M的位置坐标是以时间t为参数的函数：

$\left\{\begin{array}{c}x=rcos\mathrm{\ω}t\\ y=rsin\mathrm{\ω}t\\ z={\∫}_0^{t}v\left(u\right)du\end{array}\right.---(1-2)$

当v(t)＝1，即单位速度时

$\left\{\begin{array}{c}x=rcos\mathrm{\ω}t\\ y=r\sin \mathrm{\ω}t\\ z={\∫}_0^{t}1du=t\end{array}\right.---(1-3)$

物体M在x-z平面的投影即为r cosωt。因此可以认为余弦信号是做方程(1-3)运动在x-z平面的投影。当物体M做方程(1-2)或方程(1-3)运动沿着z轴在x-y平面的投影时，即如图1所示。由此可以认为物体M做两种形式的运动，一种是沿着z轴的直线运动，另一种是绕着z轴的圆周运动。根据振动理论，频率量在圆周运动下才是有意义的，表达的是往复快慢。

以上分析中，角速度ω和半径r均为常数，这样形成的信号是平稳信号，其瞬时频率处处相等。但在实际中，物体M绕z轴转动的半径和角速度并非上述余弦信号物理模型中恒定半径r、均匀角速度ωt，而是a(t)和θ(t)，其运动方程为：

$\left\{\begin{array}{c}x=a\left(t\right)cos\mathrm{\θ}\left(t\right)\\ y=a\left(t\right)sin\mathrm{\θ}\left(t\right)\\ z={\∫}_0^{t}1du=1\end{array}\right.---(1-4)$

其在x-y平面投影如图2所示。

二、AM-FM类信号要满足的条件

方程(1-4)中的a(t)与θ(t)体现了非平稳信号随时间变化的根本特征，因而是分析非平稳信号的基础。但由图2可以直观的看出，欲使方程(1-4)能表征信号的交变特征，调幅a(t)相对于被调制信号cosθ(t)应为缓变信号，其取值不应影响cosθ(t)的单调性，即x(t)与cosθ(t)应具有同向单调性，亦即

x′(t)·(cosθ(t))′≥0 (2-1)

根据式(2-1)及上述的物理模型，调制信号a(t)应该大于零。因此，这里的AM-FM信号必须满足下列条件：

(1)a(t)＞0；

(2)θ(t)≥0，但不要求单调；

(3)a(t)的Fourier分析所得最高频率必须小于cosθ(t)的最低频率

(4)x′(t)·(cosθ(t))′≥0

三、AM-FM类信号瞬时频率的数学公式

与简单单调函数相比，AM-FM类信号表示振动模态信号，可以表示成：

x(t)＝a(t)cos[θ(t)] (3-1)

根据上述分析，a(t)是缓变信号，即若a(t)的频谱为A(f)，cos[θ(t)]是载波信号，其频谱为B(f)，则存在常数f₀，当|f|＞f₀时，A(f)＝0；当|f|＜f₀，B(f)＝0。

由上述分析可知，AM-FM信号可通过对过零点之处求取模态信号运动速度，即真实角速度时刻，而角速度与频率满足下式：

ω＝2π·f (3-2)

因此，瞬时频率的求解转化为振动模态信号瞬时速度的求解，而瞬时速度的求解可以利用一阶导数来求解。任给一个AM-FM信号，对式(3-1)求其一阶导数得：

x′(t)＝a′(t)cos[θ(t)]-θ′(t)a(t)sin[θ(t)] (3-3)

假设式(3-1)过零点时刻为t₀，那么可以确定一系列时间点t_0i，并且t_0i满足条件：

x(t_0i)＝0，i＝1，2，…，n (3-4)

由于包络函数a(t)≠0，因此只能是cos[θ(t_0i)]＝0，又因为cos²t_0i+sin²t_0i＝1，所以sin[θ(t_0i)]＝±1。再将t_0i代入式(3-3)，得

x′(t_0i)＝a′(t_0i)cos[θ(t_0i)]-θ′(t_0i)a(t_0i)sin[θ(t_0i)]＝±θ′(t_0i)a(t_0i) (3-5)

式(3-5)中θ′(t_0i)就是相位一阶导数，而相位一阶导数正是振动信号在t_0i时刻的瞬时速度，其中的±表示的是AM-FM信号相隔半个周期在相邻过零点时刻速度方向恰好相反，这正是质点过y轴时的速度方向体现。因此式(3-5)的两边同时除以2π·a(t_0i)，且取绝对值就可得AM-FM信号在时刻t_0i处的瞬时频率。另假设AM-FM信号所有的零点时刻为t₀₁，t₀₂，…，t_0n，则AM-FM信号在t₀₁，t₀₂，…，t_0n处的瞬时频率为f_i＝|x′_a(t_0i)/2πa(t_0i)|，i＝1，2，…，n，将这些在零点时刻得到的瞬时频率f_i用三次样条函数拟合即得整个时间段上的瞬时频率函数f(t)。

四、IMF信号是满足载波最低频率高于包络最高频率的信号

内蕴模态函数(Intrinsic mode function，简称IMF)是Hilbert-Huang变换的重要概念。IMF是经过经验模态分解方法(Empirical mode decompositionmethod，简称EMD)，对非平稳、非线性信号分解得到的一系列振动模态信号。IMF满足如下两个条件，(参见文献“HUANG N E.The empirical modedecomposition and Hilbert spectrum for nonlinear and non-stationary timeseries analysis[M].Proc R Soc London，1998.”)

(1)整个数据段内极值点的个数和零交叉点的个数必须相等或相差最多不能超过一个。

(2)在任何一点，由局部极大值点形成的包络线和由局部极小值点形成的包络线的平均值为零。

由第二个条件可以推知，局部极大值可以认为是采样所得，通过对这些采样值使用三次样条函数形成的上包络可以认为是对被采样函数的恢复；同理，下包络也是如此。而上、下包络之和为零，这说明上、下包络为同一个函数，仅相差一个符号。又根据Nyquist采样定理及恢复定理，采样频率f_s与被采样信号函数最高频率f_max必须满足：

f_s≥2f_max (4-1)

因此，可以近似认为IMF是最低载波频率与IMF包络最高频率满足式(4-1)，是一个特殊的AM-FM信号。

IMF信号可写成下式：

x(t)＝a(t)cos[θ(t)] (4-2)

其中a(t)是缓变信号。根据文献“Bedrosian E.A product theorem for hilberttransform[J].1963，Proc IEEE，51：868～869”，式(4-2)满足下式：

H[x_a(t)]＝a(t)sin[θ(t)] (4-3)

其中，H[·]表示Hilbert变换。

因此，对于IMF同样可以采用本发明的方法对其瞬时频率进行计算。

五、AM-FM信号移相处理

从公式(3-4)及图2可知，本发明方法是通过拟合过零点时刻(即过x轴时)一阶导数的方法获取整个瞬时频率函数f(t)的。同样，若能获得过y轴时刻的瞬时角速度，如图3所示，则拟合时数据点数将增加一倍。这实际上就是求与原信号x(t)＝a(t)cos[θ(t)]相正交的信号，即运动物体M向y轴投影的运动，移相本发明实施例中，采用Hilbert变换来进行移相

Hilbert变换定义式为：

$H\[x\left(t\right)\]=\hat{x}\left(t\right)=x\left(t\right)*\frac{1}{\mathrm{\π}t}=\frac{1}{\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}---(5-1)$

对于形如x(t)＝a(t)cos[θ(t)]的AM-FM信号，其中a(t)是缓变信号，即若a(t)的频谱为A(f)，cos[θ(t)]的频谱为B(f)，则存在常数f₀，当|f|＞f₀时，A(f)＝0；当|f|＜f₀，B(f)＝0。根据文献“Bedrosian E.A product theorem for hilberttransform[J].1963，Proc IEEE，51：868～869”，则有下列性质：

$H\[x\left(t\right)\]=\hat{x}\left(t\right)=a\left(t\right)\sin \[\mathrm{\θ}\left(t\right)\]---(5-2)$

通过上述分析可知拟合时，获取的数据点愈多则拟合愈准确。通过观察图2、图3可知，当运动物体M向图4中1、2、3、4、5轴分别投影时，就等价于AM-FM信号分别移相k＝0，1，…，n-1，其中n≥1。这样，通过求取n路过零点时刻的瞬时频率，将会使待拟合点数增加n-1倍，拟合准确度也大大提高。

为了获取图4中向1、2、3、4、5轴投影的信号，或是为了获取向与x轴夹角成φ的轴投影的信号，即获取移相-φ的信号，如下所示：

x_φ(t)＝a(t)cos[θ(t)-φ] (5-3)

将式(5-3)展开，得

x_φ(t)＝a(t)cos[θ(t)]cosφ+a(t)sin[θ(t)]sinφ (5-4)

又根据式(5-2)，式(5-4)整理得：

$\begin{array}{c}{x}_{\mathrm{\φ}}\left(t\right)=a\left(t\right)cos\[\mathrm{\θ}\left(t\right)\]cos\mathrm{\φ}+a\left(t\right)sin\[\mathrm{\θ}\left(t\right)\]sin\mathrm{\φ}\\ =x\left(t\right)cos\mathrm{\φ}+\hat{x}\left(t\right)sin\mathrm{\φ}\end{array}---(5-5)$

式(5-5)中的恰好就是原信号x(t)的Hilbert变换所得信号。由此，任意移相-φ的信号可以通过原信号x(t)与的加权和所得。

信号x(t)＝cos(πt²)分别移相为0，的结果如图5所示。

六、离散AM-FM瞬时频率实现

实际中，得到的都是AM-FM信号的离散数据，假设以周期T对式(3-1)进行抽样，得到离散序列x(n)，设其长度为N，则有：

x(n)＝a(n)cos[θ(n)] (6-1)

常用的离散微分有前向有限微分(FFD)，后向有限微分(BFD)和中心有限微分分别定义为式(6-2)，(6-3)，(6-4)：

x′(n)＝x(n+1)-x(n) (6-2)

x′(n)＝x(n)-x(n-1) (6-3)

$x^{\′}\left(n\right)=\frac{1}{2}\[x(n+1)-x(n-1)\]---(6-4)$

为寻找所有过零点时刻处的一阶导数，首先必须找到所有过零点时刻序列t₁，t₂，…，t_n。对任意相邻的离散序列对{x(n)，x(n+1)}的乘积有两种情况：

x(n)·x(n+1)＜0，则在过(n，x(n))和(n+1，x(n+1))之间必有一个过零点时刻t_i，采用式(6-2)求t_i处的一阶导数；

x(n)·x(n+1)＝0，则判断x(n)或x(n+1)等于0，采用式(6-4)求过零点时刻t_i处的一阶导数。

通过公式f_i＝|x′_a(t_i)/(2π·a(t_i))|，i＝1，2，…，n可求出过零点处瞬时频率。但信号x(t)＝a(t)cos[θ(t)]在过零点处的包络值是未知的，信号x(t)＝a(t)cos[θ(t)]的Hilbert变换结果为：

$\hat{x}\left(t\right)=a\left(t\right)\sin \[\mathrm{\θ}\left(t\right)\]---(6-5)$

分析x(t)，可知，x(t)过零点处的包络值恰好是的极值。因此，过零点瞬时频率为：

f_i＝|x′(t_i)/2πa(t_i)|，i＝1，2，…，n

其中a(t_i)是的t＝t_i处局部极大值；反之，过零点处的包络值恰好是x(t)的极值。

同理，分别计算原信号以及n-1个移相信号在过零点时刻的瞬时频率f_i，然后用三次样条函数拟合，即可得整个时间段上的瞬时频率函数f(t)。

七、瞬时频率计算的计算方法

由前述分析可知，本发明的瞬时频率计算方法可归纳如下：

(1)输入AM-FM信号x(t)＝a(t)cos[θ(t)]，其中a(t)为输入信号的幅值，θ(t)为输入信号的相位；

(2)求取输入信号x(t)＝a(t)cos[θ(t)]的局部极大值和极小值及发生时刻，记为(t_j，a(t_j))；

(3)根据步骤(2)得到的各不同时刻的a(t_j)值，采用样条插值函数拟合所述的a(t_j)绝对值，得到输入信号的包络函数a(t)；

(4)求取输入信号x(t)＝a(t)cos[θ(t)]的一阶导数x′(t)；

(5)求取x′(t)的局部极大、极小值及发生时刻，记为(t_i，x′(t_i))；

(6)根据公式计算t_i时刻输入信号的瞬时频率f(t_i)；

(7)对输入信号x(t)＝a(t)cos[θ(t)]进行移相处理，分别得到x₁(t)，x₂(t)，…x_n-1(t)，计算公式为：

$x_i\left(t\right)=x\left(t\right)cos\left(\frac{i\·\mathrm{\π}}{n}\right)+H\[x\left(t\right)\]\sin \left(\frac{i\·\mathrm{\π}}{n}\right),i=1,2,...n-1,$

其中H[·]为Hilbert变换，n≥1；

(8)对于步骤(7)得到的x₁(t)，x₂(t)，…x_n-1(t)，重复步骤(2)～(6)，得到(t_i，f(t_i))点序列；

(9)根据步骤(8)得到的(t_i，f(t_i))点序列，采用样条插值函数拟合所述的f(t_i)，得到输入信号x(t)＝a(t)cos[θ(t)]的瞬时频率函数f(t)。

其具体流程如图6所示。

八、实例说明

下面分别以正弦信号、典型调频信号、指数调频信号、复合信号为例进行说明。

1、正弦信号

x(t)＝sin(2π·t) t∈[0，6] (8-1)

式(8-1)是恒定频率f＝1。由第三节可知，当频率是恒值时，是均匀采样，图7很好地体现了均匀采样特性，且拟合效果也比较满意。

2、典型调频信号

x(t)＝sin(10πt²) t∈[0，1.8] (8-2)

3、指数调频信号

x(t)＝cos[2π·exp(3t)] t∈[0，1.0] (8-3)

式(8-2)频率表达式为f(t)＝(10πt²)′/2π＝10t，其最大频率为18，式(8-3)频率表达式为f(t)＝(2π·exp(3t))′/2π＝3exp(3t)，f∈[3，60.26]。同样由第三节可知，由于相位函数是非线性的，因此，通过信号过零点处求得的瞬时频率序列是非均匀采样所得。图8、图9中的“+”也很好地体现非均匀采样特点，且拟合效果也比较满意。

4、复合调频信号

$x\left(t\right)=\{\begin{array}{cc}\sin \left(2\mathrm{\π}\right(5{\left(1.4-t\right)}^2+2\left(1.4-t\right)\left)\right),& t\∈\[0,1.4\]\\ -\sin \[2\mathrm{\π}\·\exp \left(2\right(t-1.4\left)\right)\]& t\∈(1.4,2.8\]\end{array}---(8-4)$

式(8-4)是两段调频信号结合而成，其频率是先快到慢，然后在从慢到快。其频率表达式为

$f\left(t\right)=\left\{\begin{array}{cc}{\left(2\mathrm{\π}\right(5{\left(1.4-t\right)}^2+2\left(1.4-t\right)\left)\right)}^{\′}/2\mathrm{\π}=-10(1.4-t)-2& f\∈\[-16,-2\]\\ (2\mathrm{\π}\·\exp (2\left(t-1.4\right)\left)\right)/2\mathrm{\π}=2\exp (2\·(t-1.4\left)\right)& f\∈(2,32.89)\end{array}\right.$

由于取频率的绝对值，因此，频率是先从16变到2，而后再从2变到32.89。前半段频率是均匀变化，后半段是非均匀变化。由于两段信号相位函数均为非线性函数，体现在采样点上就是非均匀采样。图10中的“+”点也体现了非均匀采样。拟合时，只是在两段结合处出现了误差。

为了说明本发明的效果，将以上给出的四个调频信号，通过本方法、HHT方法分别求得其瞬时频率，并和调频信号真实频率进行比较，比较结果分别示于图11～图14中。从这4个图中可以看到，拟合成的曲线与真实频率曲线很接近，而采用HHT方法求得的瞬时频率曲线存在明显的Gibbs现象，尤其在信号开始和末尾段，在复合信号结合部Gibbs现象更是明显。采用HHT方法求取瞬时频率时，出现的这种Gibbs现象是由于其采用了离散Hilbert变换造成的，而采用在AM-FM信号过零点处求取导数方法获得的瞬时频率没有出现Gibbs现象。尽管HHT方法求取的瞬时频率曲线能反映真实瞬时频率的大概走势，但其较大的波动也带来了较大的误差，而过零点处求取导数方法获得的瞬时频率曲线不仅能反映真实瞬时频率的大概走势，而且无Gibbs现象，且比HHT方法误差更小。

九、实际应用

本发明的最佳实施应用分两种情况：空-时联合采样方法和分数阶Hilbert变换方法进行讨论。

1、空-时联合采样方法

当齿轮、滚动轴承、汽轮机等做旋转运动或往复机械运动出现故障时，振动信号中往往呈现调制特征，调制频率往往反映故障特征。此时，对于这类信号，最宜采用多路数据同时采集的方法。同时采集n路旋转信号，1～n-1路相当于对0路做AM-FM信号移相处理所得。直接采集1～n-1路信号可避免对AM-FM信号做移相处理带来的计算误差。如图15所示，在这4个位移传感器采集数据时，必须使其4个同时工作，即起始时刻要尽可能的相同，否则采集到的信号就不等价于对做AM-FM信号做移相处理所得的结果。只有起始时刻相同，才能得到正确的AM-FM不同移相的信号。

当分别获得同一个AM-FM信号的n路后，可对每一路分别用第二节给出的算法分别求出极大值序列、瞬时频率点序列，然后将n路极大值序列按时间先后顺序排列，拟合形成包络a(t)；再将n路瞬时频率点序列按时间先后顺序排成列，拟合成瞬时频率函数f(t)。

2、移相处理实现多路信号

当要处理的信号不是旋转体的信号时，此时就不能通过图15的方式获取AM-FM信号的n路表达了。此时只能得到形如x_a(t)＝a(t)cos[θ(t)]的信号，对此可以通过“第五节AM-FM信号移相处理”求取1～n-1路的信号。

当获得信号x(t)＝a(t)cos[θ(t)]分别移相k＝1，2，…，n-1的信号后，可按第七节求取瞬时频率的方法进行。

本发明说明书中未作详细描述的内容属本领域专业技术人员的公知技术。

Claims (2)

1.基于微分的AM-FM信号瞬时频率计算方法，其特征在于步骤如下：

(1)输入AM-FM信号x(t)＝α(t)cos[θ(t)]，其中a(t)为输入信号的幅值，θ(t)为输入信号的相位；所述的输入信号x(t)＝a(t)cos[θ(t)]应满足的要求为：a(t)＞0；θ(t)≥0，但不要求单调；a(t)的Fourier分析所得最高频率必须小于cosθ(t)的Fourier分析所得最低频率x′(t)·(cosθ(t))′≥0；

(2)求取输入信号x(t)＝a(t)cos[θ(t)]的局部极大值和极小值及发生时刻，记为(t_j，a(t_j))；

(3)根据步骤(2)得到的各不同时刻的a(t_j)值，采用样条插值函数拟合所述a(t_j)的绝对值，得到输入信号的包络函数a(t)；

(4)求取输入信号x(t)＝α(t)cos[θ(t)]的一阶导数x′(t)；

(5)求取x′(t)的局部极大值和极小值及发生时刻，记为(t_i，x′(t_i))；

(6)根据公式计算t_i时刻输入信号的瞬时频率f(t_i)；

(7)对输入信号x(t)＝a(t)cos[θ(t)]进行移相处理，分别得到x₁(t)，x₂(t)，…x_n-1(t)，其中n≥1；

(8)对于步骤(7)得到的x₁(t)，x₂(t)，…x_n-1(t)，重复步骤(2)～(6)，得到(t_i，f(t_i))点序列；

(9)根据步骤(8)得到的(t_i，f(t_i))点序列，采用样条插值函数拟合所述的f(t_i)，得到输入信号x(t)＝a(t)cos[θ(t)]的瞬时频率函数f(t)。

2.根据权利要求1所述的基于微分的AM-FM信号瞬时频率计算方法，其特征在于：所述步骤(7)中对输入信号x(t)＝a(t)cos[θ(t)]进行移相处理的公式为：

$x_i\left(t\right)=x\left(t\right)cos\left(\frac{i\·\mathrm{\π}}{n}\right)+H\[x\left(t\right)\]sin\left(\frac{i\·\mathrm{\π}}{n}\right),i=1,2,...n-1,$

其中为Hilbert变换，n≥1。