在信號處理 中,觀察信號的瞬時頻率 是很重要的課題。假設一實信號
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.