在信号处理 中,观察信号的瞬时频率 是很重要的课题。假设一实信号
x
(
t
)
{\displaystyle x(t)\,}
可写成指数信号的N项相加(有无穷多种表示法,以
N
{\displaystyle N\,}
小的为宜),即
x
(
t
)
=
∑
k
=
1
N
a
k
⋅
e
j
⋅
ϕ
k
(
t
)
{\displaystyle x(t)=\sum _{k=1}^{N}a_{k}\cdot e^{j\cdot \phi _{k}(t)}}
, 其中
a
k
{\displaystyle a_{k}\,}
为虚常数。
则瞬时频率 (以频率表示)
f
k
(
t
)
=
ϕ
k
′
(
t
)
2
π
{\displaystyle f_{k}(t)={\frac {\phi _{k}^{\prime }(t)}{2\pi }}}
, k=1,...,N
以频率来表示(单位为赫兹 ):
f
k
(
t
)
=
ϕ
k
′
(
t
)
2
π
{\displaystyle f_{k}(t)={\frac {\phi _{k}^{\prime }(t)}{2\pi }}}
, k=1,...,N
以角频率 来表示(单位为弧度每秒):
f
k
(
t
)
=
ϕ
k
′
(
t
)
{\displaystyle f_{k}(t)=\phi _{k}^{\prime }(t)}
, k=1,...,N
直观上,瞬时频率为相位的微分。对于自然界中的实数讯号,如何定义讯号的相位。Gabor提出解析讯号法(Analytic Signal Method),将实数讯号表示为对应的复数讯号,即可定义复数讯号的大小与相位,将实数讯号的瞬时频率求出。
实数讯号
s
(
t
)
{\displaystyle s\left(t\right)}
的解析讯号(Analytic Signal)
z
(
t
)
{\displaystyle z\left(t\right)}
定义为
解析函数的极座标表示 瞬时相位 瞬时频率
z
(
t
)
=
s
(
t
)
+
j
π
∫
−
∞
∞
s
(
τ
)
t
−
τ
d
τ
.
{\displaystyle z(t)=s(t)+{\frac {j}{\pi }}\int _{-\infty }^{\infty }{\frac {s(\tau )}{t-\tau }}\,d\tau .\,}
其中虚数项为实数讯号
s
(
t
)
{\displaystyle s\left(t\right)}
的希尔伯特转换 (Hilbert Transform),将它定义为
s
^
(
t
)
{\displaystyle {\widehat {s}}\left(t\right)}
。称作解析函数的理由是,此型式的复数函数满足柯西-里曼(Cauchy-Riemann)的可微分条件,称之为解析函数(Analytic Function)。因此,解析讯号
z
(
t
)
{\displaystyle z\left(t\right)}
可以表示为
z
(
t
)
=
s
(
t
)
+
j
s
^
(
t
)
=
m
(
t
)
⋅
e
j
θ
(
t
)
{\displaystyle z(t)=s(t)+j{\widehat {s}}\left(t\right)=m(t)\cdot e^{j\theta (t)}\,}
其中
m
(
t
)
=
s
2
(
t
)
+
s
^
2
(
t
)
{\displaystyle m(t)={\sqrt {s^{2}(t)+{\widehat {s}}^{2}(t)}}}
;
θ
(
t
)
=
arctan
(
s
^
(
t
)
s
(
t
)
)
{\displaystyle \theta (t)=\arctan \left({{\widehat {s}}(t) \over s(t)}\right)\,}
如果
s
(
t
)
{\displaystyle s\left(t\right)}
是没有任何限制条件的时间讯号,计算出来的瞬时频率可能不是正确的结果。对于平均值为零的局部对称讯号而言,前述定义的瞬时频率才具有物理意义。在1998年,黄锷 (Norden E. Huang)博士提出一个有效的演算法,将讯号先行分解成具有局部对称之分量,以正确地求得资料的瞬时频率。这个方法称为希尔伯特-黄转换 (Hilbert Huang Transform, HHT)。
以下简单的例子来说明,对于平均值为零的讯号,此瞬时频率的定义才具有物理意义。对于一个弦波讯号
s
(
t
)
{\displaystyle s\left(t\right)}
,
s
(
t
)
=
β
+
cos
(
t
)
{\displaystyle s\left(t\right)=\beta +\cos {(t)}}
考虑三种情况: (1)
β
=
0
{\displaystyle \beta =0}
(2)
0
<
β
<
1
{\displaystyle 0<\beta <1}
(3)
β
>
1
{\displaystyle \beta >1}
(1)
β
=
0
{\displaystyle \beta =0}
: 当弦波讯号平均值为零时,
z
(
t
)
{\displaystyle z\left(t\right)}
在复数平面上的描述是以座标原点为中心的单位圆,它的相位角
θ
(
t
)
{\displaystyle \theta (t)}
则是以座标原点为中心,反时针方向 呈线性递增,其图形为斜率1的直线,而瞬时频率是一个常数值。
(2)
0
<
β
<
1
{\displaystyle 0<\beta <1}
:
z
(
t
)
{\displaystyle z\left(t\right)}
在复数平面上仍然是一个单位圆,但其圆心从原点偏移了
β
{\displaystyle \beta }
个单位,其相角
θ
(
t
)
{\displaystyle \theta (t)}
不再呈现线性递增,瞬时频率出现震荡的现象,而不是应有的常数值。
(3)
β
>
1
{\displaystyle \beta >1}
: 因为
β
{\displaystyle \beta }
值超过单位圆的半径,因此
z
(
t
)
{\displaystyle z\left(t\right)}
的圆心在单位圆之外。如此相位角
θ
(
t
)
{\displaystyle \theta (t)}
在[
−
π
{\displaystyle -\pi }
/
2
{\displaystyle 2}
,
π
{\displaystyle \pi }
/
2
{\displaystyle 2}
]震荡,瞬时频率出现负值,与原讯号的特性有极大的差别。
因为在目前许多数位信号处理 的应用上都与信号的频谱或信号的频宽有很大的关系。 若能确实地侦测信号的瞬时频率,则通道频宽可以被可适性(adaptive)的决定,如此一来能更有效地利用系统资源,提高系统效能。
调变(modulation)
多工方式(multiplexing)
滤波器的设计(filter design)
信号压缩(data compression)
信号分析(signal analysis)
信号辨识(signal identification)
语音信号处理(acoustical signal processing)
制作系统的模型(system modling)
雷达系统的分析(rader system analysis)
取样(sampling)
当瞬时频率为常数即
ϕ
(
t
)
{\displaystyle \phi (t)\,}
为一阶时间函数,使用傅立叶变换 做信号分析。 由于从傅立叶变换 中是无法观察出信号频谱随著时间改变的变化。 故只有当瞬时频率为常数,不是时间的函数时,便可使用傅立叶变换 做信号分析。
当瞬时频率不为常数即
ϕ
(
t
)
{\displaystyle \phi (t)\,}
为高阶时间函数,使用时频分析 做信号分析。 从时频分析 可观察出信号频率随著时间变化的改,这是傅立叶变换 无法做到的。 因此当瞬时频率为时间的函数,使用时频分析 做信号分析,如下图,可以确切地观察到信号瞬时频率的变化。
x
(
t
)
=
{
c
o
s
(
π
t
)
t
≤
10
c
o
s
(
3
π
t
)
10
<
t
≤
20
c
o
s
(
2
π
t
)
t
>
20
{\displaystyle x(t)={\begin{cases}cos(\pi t)\ \ \ t\leq 10\\cos(3\pi t)\ \ \ 10<t\leq 20\ \ \ \\cos(2\pi t)\ \ \ t>20\end{cases}}}
当信号只有一个成分,且曲线是周期性的,具有固定的振幅、频率和相,可以使用希尔伯特转换,计算信号的相位求瞬时频率。
瞬时频率
=
1
2
π
d
d
t
θ
,
θ
=
t
a
n
−
1
x
H
(
t
)
x
(
t
)
{\displaystyle \,={\frac {1}{2\pi }}{\frac {d}{dt}}\theta ,\quad \theta =tan^{-1}{\frac {x_{H}(t)}{x_{(}t)}}}
x
H
(
t
)
=
1
π
∫
−
∞
∞
x
(
τ
)
t
−
τ
d
τ
{\displaystyle x_{H}(t)={\frac {1}{\pi }}\int _{-\infty }^{\infty }{\frac {x(\tau )}{t-\tau }}d\tau \,}
举例:
c
o
s
(
2
π
f
t
)
→
H
i
l
b
e
r
t
s
i
n
(
2
π
f
t
)
{\displaystyle cos(2\pi ft)\;\xrightarrow {Hilbert} \;sin(2\pi ft)}
θ
=
t
a
n
−
1
(
t
a
n
(
2
π
f
t
)
)
=
2
π
f
t
{\displaystyle \theta =tan^{-1}(tan(2\pi ft))=2\pi ft}
1
2
π
d
d
t
θ
=
f
{\displaystyle {\frac {1}{2\pi }}{\frac {d}{dt}}\theta =f}
当信号为复数函数、非sinusoid曲线、有多个成分或
θ
{\displaystyle \theta }
有多个解,此时可以先将信号分成sinusoid-like成分和趋势,再运用Hilbert transform或短时距傅立叶变换 (STFT) 求零交叉的数量。
1. 使用经验模态分解 (EMD) 将信号分解为一系列本征模态函数 (IMFs)的振荡模态和趋势
2. 对每个IMF算瞬时频率:
使用Hilbert transform
计算STFT,瞬时频率 =
number of zero-crossings between
t
−
B
and
t
+
B
4
B
{\displaystyle {\frac {{\text{number of zero-crossings between}}\;t-B\;{\text{and}}\;t+B}{4B}}}
直接计算零交叉,零交叉 =
number of periods
2
{\displaystyle {\frac {\text{number of periods}}{2}}}
Jian-Jiun Ding, class lecture of Time Frequency Analysis and Wavelet transform, Graduate Institute of Communication Engineering, National Taiwan University, Taipei, Taiwan, 2007.
Leon Coen, Time-Frequency Analysis, Prentice Hall, 1995.
陈韦佑, "以希尔伯特-黄转换法为GPS接收机抑制调频干扰", 国立台湾大学电机工程研究所硕士论文, 2007.