加伯轉換 是窗函數 為高斯函數 的短時距傅立葉變換 。
將短時距傅立葉轉換 中的窗函數 代入高斯函數 ,即可得下面的標準定義:
G
x
(
t
,
f
)
=
∫
−
∞
∞
e
−
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
{\displaystyle G_{x}(t,f)=\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )\,d\tau }
以下是幾種常見的替代定義:
G
x
1
(
t
,
f
)
=
∫
−
∞
∞
e
−
π
(
τ
−
t
)
2
e
−
j
2
π
f
(
τ
−
t
2
)
x
(
τ
)
d
τ
{\displaystyle Gx_{1}(t,f)=\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f(\tau -{\frac {t}{2}})}x(\tau )\,d\tau }
G
x
2
(
t
,
f
)
=
2
4
∫
−
∞
∞
e
−
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
{\displaystyle Gx_{2}(t,f)={\sqrt[{4}]{2}}\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )\,d\tau }
G
x
3
(
t
,
ω
)
=
∫
−
∞
∞
e
−
(
τ
−
t
)
2
2
e
−
j
ω
τ
x
(
τ
)
d
τ
{\displaystyle Gx_{3}(t,\omega )=\int _{-\infty }^{\infty }e^{-{\frac {(\tau -t)^{2}}{2}}}e^{-j\omega \tau }x(\tau )\,d\tau }
G
x
4
(
t
,
ω
)
=
1
2
π
∫
−
∞
∞
e
−
(
τ
−
t
)
2
2
e
−
j
ω
(
τ
−
t
2
)
x
(
τ
)
d
τ
{\displaystyle Gx_{4}(t,\omega )={\sqrt {\frac {1}{2\pi }}}\int _{-\infty }^{\infty }e^{-{\frac {(\tau -t)^{2}}{2}}}e^{-j\omega (\tau -{\frac {t}{2}})}x(\tau )\,d\tau }
註:在文獻上可能會看到不同形式的加伯轉換,但本質上都是一樣的。
由於實作時,不能計算無限大的積分式子,所以根據高斯函數 會從兩側遞減的性質,我們可以將上式進一步化簡:
∵
{
e
−
π
a
2
<
0.00001
;
|
a
|
>
1.9143
e
−
a
2
/
2
<
0.00001
;
|
a
|
>
4.7985
{\displaystyle \because {\begin{cases}\ e^{-\pi a^{2}}<0.00001;&|a|>1.9143\\\ e^{-a^{2}/2}<0.00001;&|a|>4.7985\end{cases}}}
∴
{
G
x
(
t
,
f
)
≃
∫
t
−
1.9143
t
+
1.9143
e
−
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
G
x
4
(
t
,
ω
)
≃
1
2
π
∫
t
−
4.7985
t
+
4.7985
e
−
(
τ
−
t
)
2
2
e
−
j
ω
(
τ
−
t
/
2
)
x
(
τ
)
d
τ
{\displaystyle \therefore {\begin{cases}\ G_{x}(t,f)\simeq \int _{t-1.9143}^{t+1.9143}e^{-\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )\,d\tau \\\ Gx_{4}(t,\omega )\simeq {\sqrt {\frac {1}{2\pi }}}\int _{t-4.7985}^{t+4.7985}e^{-{\frac {(\tau -t)^{2}}{2}}}e^{-j\omega (\tau -t/2)}x(\tau )\,d\tau \end{cases}}}
其他窗函數 的短時距傅立葉變換 ,如利用方型窗函數的短時距傅立葉變換 ,無法同時兼顧時間軸和頻率軸的解像度;一者解像度提升,另一者解像度必定下降。但高斯函數由海森堡測不準原理 可得知,是最能同時讓兩軸兼顧解像度的窗函數(將於下面章節詳述)。
高斯函數 為傅立葉轉換 的特徵函數:
∫
−
∞
∞
e
−
π
t
2
e
−
j
2
π
f
τ
d
t
=
e
−
π
f
2
{\displaystyle \int _{-\infty }^{\infty }e^{-\pi t^{2}}e^{-j2\pi f\tau }dt=e^{-\pi f^{2}}}
∫
−
∞
∞
e
−
t
2
2
e
−
j
ω
t
d
t
=
e
−
f
2
2
{\displaystyle \int _{-\infty }^{\infty }e^{-{\frac {t^{2}}{2}}}e^{-j\omega t}dt=e^{-{\frac {f^{2}}{2}}}}
因此經過轉換後其性質不變。因此可讓加伯轉換後在時間軸和頻率軸的性質相互對稱。
上述提到,高斯函數 是最能兼顧時間與頻率解像度的窗函數。我們利用這個章節來詳細討論。
對於一個信號
x
(
t
)
{\displaystyle x(t)}
,當
|
t
|
→
∞
{\displaystyle |t|\to \infty }
,若
t
x
(
t
)
=
0
{\displaystyle {\sqrt {t}}x(t)=0}
,則
σ
t
σ
f
≥
1
/
4
π
{\displaystyle \sigma _{t}\sigma _{f}\geq 1/4\pi \,}
其中
σ
t
2
=
∫
(
t
−
μ
t
)
2
P
x
(
t
)
d
t
,
σ
f
2
=
∫
(
f
−
μ
f
)
2
P
X
(
f
)
d
f
{\displaystyle \sigma _{t}^{2}=\int (t-\mu _{t})^{2}P_{x}(t)dt,\sigma _{f}^{2}=\int (f-\mu _{f})^{2}P_{X}(f)df\,}
μ
t
=
∫
t
P
x
(
t
)
d
t
,
μ
f
=
∫
f
P
X
(
f
)
d
f
{\displaystyle \qquad \mu _{t}=\int tP_{x}(t)dt\quad \quad \quad \ ,\mu _{f}=\int fP_{X}(f)df}
P
x
(
t
)
=
|
x
(
t
)
|
2
∫
|
x
(
t
)
|
2
d
t
,
P
X
(
f
)
=
|
X
(
f
)
|
2
∫
|
X
(
f
)
|
2
d
f
{\displaystyle \qquad P_{x}(t)={\frac {|x(t)|^{2}}{\int |x(t)|^{2}dt}}\quad \;\;\;,P_{X}(f)={\frac {|X(f)|^{2}}{\int |X(f)|^{2}df}}}
由於兩者標準差相乘有下限,這個定理說明了我們沒有辦法同時精準量測時間和頻率,其中一者標準差下降(解像度上升),另一者標準差就上升(解像度下降)。
加伯轉換後的結果,橫軸是時間(秒),縱軸是頻率(赫茲)
當信號
x
(
t
)
{\displaystyle x(t)}
為高斯函數 時
x
(
t
)
=
e
π
t
2
,
X
(
f
)
=
e
−
π
f
2
{\displaystyle x(t)=e^{\pi t^{2}},X(f)=e^{-\pi f^{2}}}
套用以上函式求得變異數(其中由於高斯函數 為偶對稱函數,所以其
μ
t
=
0
{\displaystyle \mu _{t}=0}
)
∵
σ
t
2
=
∫
t
2
|
x
(
t
)
|
2
d
t
∫
|
x
(
t
)
|
2
d
t
=
(
1
2
)
5
2
π
1
2
=
1
4
π
{\displaystyle \because \sigma _{t}^{2}={\frac {\int t^{2}|x(t)|^{2}dt}{\int |x(t)|^{2}dt}}={\frac {({\frac {1}{2}})^{\frac {5}{2}}\pi }{\frac {1}{\sqrt {2}}}}={\frac {1}{4\pi }}}
藉由微積分公式可得:
∫
−
∞
∞
|
x
(
t
)
|
2
d
t
=
∫
−
∞
∞
e
−
2
π
t
d
t
=
1
2
{\displaystyle \int _{-\infty }^{\infty }|x(t)|^{2}dt=\int _{-\infty }^{\infty }e^{-2\pi t}dt={\frac {1}{\sqrt {2}}}}
及
∫
−
∞
∞
t
2
|
x
(
t
)
|
2
d
t
=
∫
−
∞
∞
t
2
e
−
2
π
t
2
d
t
=
Γ
(
3
/
2
)
(
2
π
)
3
2
{\displaystyle \int _{-\infty }^{\infty }t^{2}|x(t)|^{2}dt=\int _{-\infty }^{\infty }t^{2}e^{-2\pi t^{2}}dt={\frac {\Gamma (3/2)}{(2\pi )^{\frac {3}{2}}}}}
∴
σ
t
σ
f
=
1
4
π
{\displaystyle \therefore \sigma _{t}\sigma _{f}={\frac {1}{4\pi }}}
即高斯函數滿足測不準定理的最下限,所以是所有窗函數中能使時間和頻率兩者解像度都達到最高的函數。
變形的高斯函數同樣會滿足測不準原理的下限,如以下例子:
e
−
π
(
t
−
t
0
)
2
{\displaystyle e^{-\pi (t-t_{0})^{2}}}
:對概率分佈做位移,標準差不會改變。
A
e
−
π
t
2
{\displaystyle Ae^{-\pi t^{2}}}
:分子與分母同乘A,可消掉。因此標準差不會改變。
e
−
π
t
2
e
j
2
π
f
0
t
{\displaystyle e^{-\pi t^{2}}e^{j2\pi f_{0}t}}
:在時域乘上
e
j
2
π
f
0
t
{\displaystyle e^{j2\pi f_{0}t}}
相當於在頻域對頻率做位移,標準差一樣不會改變。
e
−
π
b
t
2
{\displaystyle e^{-\pi bt^{2}}}
:在時域做縮放,頻域會做相反的縮放,因此標準差也不會改變。
x
(
t
)
=
{
cos
(
2
π
t
)
;
t
<
10
cos
(
6
π
t
)
;
10
≤
t
<
20
{\displaystyle x(t)={\begin{cases}\cos(2\pi t);&t<10\\\cos(6\pi t);&10\leq t<20\end{cases}}}
右圖為即加伯轉換的結果,可以看出其時間和頻率都維持相當程度的解像度。
以下提供一個簡單的範例來比較加伯轉換以及利用方形窗函數的短時傅立葉轉換:
x
(
t
)
=
{
cos
(
10
π
t
)
;
t
<
10
cos
(
12
π
t
)
;
10
≤
t
<
20
cos
(
9
π
t
)
;
20
≤
t
<
30
{\displaystyle x(t)={\begin{cases}\cos(10\pi t);&t<10\\\cos(12\pi t);&10\leq t<20\\\cos(9\pi t);&20\leq t<30\end{cases}}}
方形窗函數短時傅立葉轉換(橫軸:時間, 縱軸:頻率)
加伯轉換(橫軸:時間, 縱軸:頻率)
從圖中可以發現方形窗函數的短時傅立葉轉換會有能量擴散的情形,而加伯轉換則是清晰的時頻圖。
由於高斯窗函數 的寬度可以由一常數做調整,因此我們將這個參數加入加伯轉換的數學式子中,讓轉換更加彈性,如下式:
G
x
(
t
,
f
)
=
σ
4
∫
−
∞
∞
e
−
σ
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
{\displaystyle G_{x}(t,f)={\sqrt[{4}]{\sigma }}\int _{-\infty }^{\infty }e^{-\sigma \pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )d\tau }
而根據前面章節所述。實作時,不能計算無限大的積分式子,所以根據高斯函數 會從兩側遞減的性質,我們可以將上式進一步化簡:
G
x
(
t
,
f
)
≃
σ
4
∫
t
−
1.9143
σ
t
+
1.9143
σ
e
−
σ
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
{\displaystyle G_{x}(t,f)\simeq {\sqrt[{4}]{\sigma }}\int _{t-{\frac {1.9143}{\sqrt {\sigma }}}}^{t+{\frac {1.9143}{\sqrt {\sigma }}}}e^{-\sigma \pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )d\tau }
根據傅立葉轉換的縮放公式,假設
w
(
t
)
=
e
−
π
σ
t
2
{\displaystyle w(t)=e^{-\pi \sigma t^{2}}}
,則傅立葉轉換後為
W
(
f
)
=
1
σ
e
−
π
f
2
σ
{\displaystyle W(f)={\frac {1}{\sqrt {\sigma }}}e^{-{\frac {\pi f^{2}}{\sigma }}}}
,使其能根據需求而調整時域解像度或頻域解像度
改變高斯函數的寬度,和改變方形窗函數短時距傅立葉變換的效果類似。若選取較大的
σ
{\displaystyle \sigma }
,時域的高斯窗函數較窄,則時域有較高的解像度,而頻域的高斯窗函數較寬,所以頻域的解像度會下降(通常用於需要時域解像度較高的應用,例如:音樂訊號);反之,若選取較小的
σ
{\displaystyle \sigma }
,時域的高斯窗函數較寬,則時域的解像度下降,而頻域的高斯窗函數較窄,所以頻域的解像度會上升(通常運用在需要頻域解像度較高的應用,例如:氣候)。雖然還是有兩軸之間的解像度的犧牲,但比起其他無法滿足測不準原理下限的窗函數,加伯轉換的兩軸還是能相對維持較高的解像度。
若應用於瞬時頻率改變較劇烈的應用,則可考慮使用窗寬度隨時間而變動的加伯轉換數學式子,如下
G
x
(
t
,
f
)
=
σ
(
t
)
4
∫
−
∞
∞
e
−
σ
(
t
)
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
{\displaystyle G_{x}(t,f)={\sqrt[{4}]{\sigma (t)}}\int _{-\infty }^{\infty }e^{-\sigma (t)\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )d\tau }
當瞬時頻率變動非常快時,使用較大的
σ
{\displaystyle \sigma }
值,使其時域解像度能較高;當瞬時頻率變動很慢時,使用較小的
σ
{\displaystyle \sigma }
值,使其頻域解像度能較高。
Direct Implementation [ 編輯 ]
X
(
t
,
f
)
=
∫
−
∞
∞
w
(
t
−
τ
)
x
(
τ
)
e
−
j
2
π
f
τ
d
τ
{\displaystyle X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau }
w
(
t
)
=
e
−
π
σ
t
2
{\displaystyle w(t)=e^{-\pi \sigma t^{2}}}
令
t
=
n
Δ
t
,
f
=
m
Δ
f
,
τ
=
p
Δ
t
{\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{t}}
可將式子改寫為離散形式:
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
−
∞
∞
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=-\infty }^{\infty }{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
w
(
t
)
≅
0
f
o
r
|
t
|
>
B
,
B
Δ
t
=
Q
{\displaystyle w(t)\cong 0\qquad for\left|t\right|>B,{\frac {B}{\Delta _{t}}}=Q}
w
(
(
n
−
p
)
Δ
t
)
≅
0
{\displaystyle w((n-p)\Delta _{t})\cong 0\qquad }
f
o
r
|
n
−
p
|
>
B
Δ
t
{\displaystyle for\left|n-p\right|>{\frac {B}{\Delta _{t}}}}
,
|
p
−
n
|
>
Q
{\displaystyle \left|p-n\right|>Q}
therefore,only when
−
Q
<
p
−
n
<
Q
{\displaystyle -Q<p-n<Q}
w
(
(
n
−
p
)
Δ
t
)
{\displaystyle w((n-p)\Delta _{t})}
is nonzero
可改寫為:
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
按照此式即可實現
e
−
π
σ
a
2
<
0.00001
{\displaystyle e^{-\pi \sigma a^{2}}<0.00001}
w
h
e
n
|
a
|
>
1.9143
{\displaystyle when\left|a\right|>1.9143}
Q
=
1.9143
σ
Δ
t
{\displaystyle Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}}
B
=
1.9143
σ
{\displaystyle B={\frac {1.9143}{\sqrt {\sigma }}}}
(1)
Δ
t
<
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
O(TFQ) T:時間取樣點數 F:頻率取樣點數 Q:
Q
=
1.9143
σ
Δ
t
{\displaystyle Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}}
優點:簡單實現,限制條件少
缺點:時間複雜度高
FFT-Based Method(快速傅立葉轉換)[ 編輯 ]
由Direct Implementation可得下式
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
令
q
=
p
−
(
n
−
Q
)
→
p
=
(
n
−
Q
)
+
q
{\displaystyle q=p-(n-Q)\to p=(n-Q)+q}
且離散傅立葉轉換標準式
Y
[
m
]
=
∑
n
=
0
N
−
1
y
[
n
]
e
−
j
2
π
m
n
N
{\displaystyle Y[m]=\sum \limits _{n=0}^{N-1}y[n]e^{-j{\frac {2\pi mn}{N}}}}
可將式子整理為:
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
e
j
2
π
(
Q
−
n
)
m
N
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{j{\textstyle {{2\pi \,(Q-n)m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}}
按照此式將
x
1
{\displaystyle {x_{1}}}
以fft()算出帶入即可實現
其中
x
1
(
q
)
=
w
(
(
Q
−
q
)
Δ
t
)
x
(
(
n
−
Q
+
q
)
Δ
t
)
{\displaystyle {x_{1}}\left(q\right)=w\left({(Q-q){\Delta _{t}}}\right)x\left({(n-Q+q){\Delta _{t}}}\right)}
,
0
≤
q
≤
2
Q
{\displaystyle 0\leq q\leq 2Q}
,
w
(
t
)
=
e
−
π
σ
t
2
{\displaystyle w(t)=e^{-\pi \sigma t^{2}}}
x
1
(
q
)
=
0
,
2
Q
<
q
≤
N
{\displaystyle {x_{1}}\left(q\right)=0,2Q<q\leq N}
Q
=
1.9143
σ
Δ
t
{\displaystyle Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}}
B
=
1.9143
σ
{\displaystyle B={\frac {1.9143}{\sqrt {\sigma }}}}
Matlab及python 皆可呼叫fft函式完成
Y
[
m
]
=
∑
n
=
0
N
−
1
y
[
n
]
e
−
j
2
π
m
n
N
{\displaystyle Y[m]=\sum \limits _{n=0}^{N-1}y[n]e^{-j{\frac {2\pi mn}{N}}}}
假設
t
=
n
0
Δ
t
,
(
n
0
+
1
)
Δ
t
,
⋯
⋯
,
(
n
0
+
T
−
1
)
Δ
t
{\displaystyle t=n_{0}\Delta _{t},(n_{0}+1)\Delta _{t},\cdots \cdots ,(n_{0}+T-1)\Delta _{t}}
f
=
m
0
Δ
f
,
(
m
0
+
1
)
Δ
f
,
⋯
⋯
,
(
m
0
+
F
−
1
)
Δ
f
{\displaystyle \,f=m_{0}\Delta _{f},(m_{0}+1)\Delta _{f},\cdots \cdots ,(m_{0}+F-1)\Delta _{f}}
step 1:計算
n
0
,
m
0
,
T
,
F
,
N
,
Q
{\displaystyle n_{0},m_{0},T,F,N,Q}
step 2:
n
=
n
0
{\displaystyle n=n_{0}}
step 3:決定
x
1
(
q
)
{\displaystyle x_{1}(q)}
step 4:
X
1
(
m
)
=
F
F
T
[
x
1
(
q
)
]
{\displaystyle X_{1}(m)=FFT[x_{1}(q)]}
step 5:轉換
X
1
(
m
)
{\displaystyle X_{1}(m)}
成
X
(
n
Δ
t
,
m
Δ
f
)
{\displaystyle X(n\Delta _{t},m\Delta _{f})}
step 6:設
n
=
n
+
1
{\displaystyle n=n+1}
and return to Step 3 until
n
=
n
0
+
T
+
1
{\displaystyle n=n_{0}+T+1}
(1)
Δ
t
<
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
(基本上任何實現方法都要避免贋頻效應)
(2)
Δ
t
Δ
f
=
1
N
{\displaystyle {\Delta _{t}}{\Delta _{f}}={\textstyle {1 \over {N}}}}
(3)
N
=
1
/
Δ
t
Δ
f
≥
2
Q
+
1
{\displaystyle N=1/{\Delta _{t}}{\Delta _{f}}\geq 2Q+1}
O
(
T
N
log
2
N
)
{\displaystyle O(TN{\log _{2}}N)}
優點:時間複雜度低
缺點:限制條件較直接實現法多
可改寫為:
由Direct Implementation可得下式
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
e
−
π
σ
a
2
<
0.00001
{\displaystyle e^{-\pi \sigma a^{2}}<0.00001}
w
h
e
n
|
a
|
>
1.9143
{\displaystyle when\left|a\right|>1.9143}
Q
=
1.9143
σ
Δ
t
{\displaystyle Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}}
B
=
1.9143
σ
{\displaystyle B={\frac {1.9143}{\sqrt {\sigma }}}}
令
e
x
p
(
−
j
2
π
m
p
Δ
t
Δ
f
)
=
e
x
p
(
−
j
π
p
2
Δ
t
Δ
f
)
e
x
p
(
j
π
(
p
−
m
)
2
Δ
t
Δ
f
)
e
x
p
(
−
j
π
m
2
Δ
t
Δ
f
)
{\displaystyle exp(-j2\pi \,mp{\Delta _{t}}{\Delta _{f}})=exp(-j\pi \,p^{2}{\Delta _{t}}{\Delta _{f}})exp(j\pi \,{(p-m)}^{2}{\Delta _{t}}{\Delta _{f}})exp(-j\pi \,m^{2}{\Delta _{t}}{\Delta _{f}})}
可將式子改寫為:
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
→
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
e
−
j
π
m
2
Δ
t
Δ
f
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
π
p
2
Δ
t
Δ
f
e
j
π
(
p
−
m
)
2
Δ
t
Δ
f
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}\to {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{-j\pi \,m^{2}{\Delta _{t}}{\Delta _{f}}}}\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j\pi \,p^{2}{\Delta _{t}}{\Delta _{f}}}}{e^{j\pi \,{(p-m)}^{2}{\Delta _{t}}{\Delta _{f}}}}}
按此式即可實現
Step1:
x
1
[
p
]
=
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
π
p
2
Δ
t
Δ
f
{\displaystyle x_{1}[p]=w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}}
n
−
Q
≤
p
≤
n
+
Q
{\displaystyle \quad \quad n-Q\leq p\leq n+Q}
Step2:
X
2
[
n
,
m
]
=
∑
p
=
n
−
Q
n
+
Q
x
1
[
p
]
c
[
m
−
p
]
c
[
m
]
=
e
j
π
m
2
Δ
t
Δ
f
{\displaystyle X_{2}[n,m]=\sum _{p=n-Q}^{n+Q}x_{1}[p]c[m-p]\quad \quad c[m]=e^{j\pi m^{2}\Delta _{t}\Delta _{f}}}
Step3:
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
e
−
j
π
m
2
Δ
t
Δ
f
X
2
[
m
,
n
]
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\Delta _{t}e^{-j\pi m^{2}\Delta _{t}\Delta _{f}}X_{2}[m,n]}
(1)
Δ
t
<
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
O
(
T
N
log
2
N
)
{\displaystyle O(TN{\log _{2}}N)}
優點:限制條件與Direct Implementation法一樣基本上沒有限制
缺點:時間複雜度與FFT-Based Method(快速傅立葉轉換)一樣
但由於加伯轉換無法使用Recursive Method(遞迴法)所以此不能算是缺點
加伯轉換的大部分的特性和方形窗函數短時距傅立葉轉換的特性都相似,有些特性甚至更加接近傅立葉轉換的特性。
當
k
≠
0
,
∫
−
∞
∞
G
x
(
t
,
f
)
e
j
2
π
k
t
f
d
f
=
e
−
π
(
k
−
1
)
2
t
2
x
(
k
t
)
{\displaystyle k\neq 0,\int _{-\infty }^{\infty }G_{x}(t,f)e^{j2\pi ktf}\,df=e^{-\pi (k-1)^{2}t^{2}}x(kt)}
當
k
=
0
,
∫
−
∞
∞
G
x
(
t
,
f
)
d
f
=
e
−
π
t
2
x
(
0
)
{\displaystyle k=0,\int _{-\infty }^{\infty }G_{x}(t,f)\,df=e^{-\pi t^{2}}x(0)}
當
k
=
1
,
∫
−
∞
∞
G
x
(
t
,
f
)
e
j
2
π
t
f
d
f
=
x
(
t
)
{\displaystyle k=1,\int _{-\infty }^{\infty }G_{x}(t,f)e^{j2\pi tf}\,df=x(t)}
(還原成原始信號)
若
y
(
t
)
=
x
(
t
−
t
0
)
{\displaystyle y(t)=x(t-t_{0})\,}
,則
G
y
(
t
,
f
)
=
G
x
(
t
−
t
0
,
f
)
e
−
j
2
π
f
t
0
{\displaystyle G_{y}(t,f)=G_{x}(t-t_{0},f)e^{-j2\pi ft_{0}}\,}
若
y
(
t
)
=
x
(
t
)
e
j
2
π
f
0
t
{\displaystyle y(t)=x(t)e^{j2\pi f_{0}t}\,}
,則
G
y
(
t
,
f
)
=
G
x
(
t
,
f
−
f
0
)
{\displaystyle G_{y}(t,f)=G_{x}(t,f-f_{0})\,}
若有一信號
h
(
t
)
=
α
x
(
t
)
+
β
y
(
t
)
{\displaystyle h(t)=\alpha x(t)+\beta y(t)\,}
,
G
h
(
t
,
f
)
,
G
x
(
t
,
f
)
,
G
y
(
t
,
f
)
{\displaystyle G_{h}(t,f),G_{x}(t,f),G_{y}(t,f)\,}
分別為
h
(
t
)
,
x
(
t
)
,
y
(
t
)
{\displaystyle h(t),x(t),y(t)\,}
做加伯轉換的結果,則
G
h
(
t
,
f
)
=
α
G
x
(
t
,
f
)
+
β
G
y
(
t
,
f
)
{\displaystyle G_{h}(t,f)=\alpha G_{x}(t,f)+\beta G_{y}(t,f)\,}
。
若
t
>
t
0
{\displaystyle t>t_{0}}
時
x
(
t
)
=
0
{\displaystyle x(t)=0}
,則
∫
−
∞
∞
|
G
x
(
t
,
f
)
|
2
d
f
<
e
−
2
π
(
t
−
t
0
)
2
∫
−
∞
∞
|
G
x
(
t
0
,
f
)
|
2
d
f
{\displaystyle \int _{-\infty }^{\infty }|G_{x}(t,f)|^{2}df<e^{-2\pi (t-t_{0})^{2}}\int _{-\infty }^{\infty }|G_{x}(t_{0},f)|^{2}df}
∫
−
∞
∞
|
G
x
(
t
,
f
)
|
2
d
f
=
∫
−
∞
∞
|
x
(
τ
)
|
2
d
τ
≃
∫
u
−
1.9143
u
+
1.9143
e
−
2
π
(
τ
−
u
)
2
|
x
(
τ
)
|
2
d
τ
{\displaystyle \int _{-\infty }^{\infty }|G_{x}(t,f)|^{2}\,df=\int _{-\infty }^{\infty }|x(\tau )|^{2}\,d\tau \simeq \int _{u-1.9143}^{u+1.9143}e^{-2\pi (\tau -u)^{2}}|x(\tau )|^{2}\,d\tau }
∫
−
∞
∞
∫
−
∞
∞
G
x
(
t
,
f
)
G
y
∗
(
t
,
f
)
d
f
d
t
=
∫
−
∞
∞
x
(
τ
)
y
∗
(
τ
)
d
τ
{\displaystyle \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }G_{x}(t,f)G_{y}^{*}(t,f)\,dfdt=\int _{-\infty }^{\infty }x(\tau )y^{*}(\tau )\,d\tau }
∫
−
∞
∞
∫
−
∞
∞
G
x
(
t
,
f
)
G
y
(
t
,
f
)
d
f
d
t
=
∫
−
∞
∞
x
(
τ
)
y
(
τ
)
d
τ
{\displaystyle \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }G_{x}(t,f)G_{y}(t,f)dfdt=\int _{-\infty }^{\infty }x(\tau )y(\tau )d\tau }
1. 當
x
(
t
)
=
δ
(
t
)
{\displaystyle x(t)=\delta (t)\,}
,
G
x
(
t
,
f
)
=
e
−
π
t
2
{\displaystyle G_{x}(t,f)=e^{-\pi t^{2}}}
2. 當
x
(
t
)
=
1
{\displaystyle x(t)=1\,}
,
G
x
(
t
,
f
)
=
e
j
2
π
f
t
e
π
f
2
{\displaystyle G_{x}(t,f)=e^{j2\pi ft}e^{\pi f^{2}}\,}
和方形窗函數短時距傅立葉轉換 不同的是,加伯轉換的結果對於時間和頻率軸較對稱,也比較沒有旁波(sidelobe);也印證了上述所說的,加伯轉換較能維持兩個軸的解像度。
最佳時間-頻率局部化特性
Gabor Transform 使用高斯窗函數,與其他常見窗函數(如Rectangle、Triangle、Hanning、Hamming)相比,滿足測不準定理的最小下限(Minimum Uncertainty Principle)。這意味着,高斯函數能夠在時間域和頻率域中同時提供最佳的解像度,避免信號特徵的模糊或失真。
高時間解像度:能捕捉信號的快速變化,對於瞬態信號(如語音中的短促音位或振動信號中的瞬時變化)尤為重要。
高頻率解像度:能精確分辨信號中的穩態頻率成分,特別適合於分析連續且平穩的周期信號。
算法穩健且實現簡單
Gabor Transform 基於傅里葉變換的數學理論,其結構清晰且實現相對簡單。現代數值計算技術(如快速傅里葉變換,FFT)的發展進一步提升了 Gabor Transform 的計算效率,使其能夠在高效實現的同時保持穩健性。
穩健性:由於其依賴於成熟的數學基礎,在實施中容易檢測和修正潛在錯誤。
實現便利性:現有的數學工具庫(如 MATLAB、Python 的 Scipy、Octave)提供了高度封裝的 Gabor Transform 函數,大幅降低了實現門檻,讓開發者能更專注於應用場景設計,而非底層算法調試。
廣泛的應用場景
語音識別:利用 Gabor Transform 提取時頻特徵,提升語音識別準確性,尤其在噪聲環境中。
圖像處理
紋理分析:有效捕捉圖像的方向與頻率特徵,用於紋理分類和圖像分割。
邊緣檢測:適用於醫學圖像和場景理解,改善邊緣檢測效果。
機械振動信號分析
故障檢測:檢測機械設備(如齒輪、軸承)異常頻率,應用於設備維護。
非平穩信號分析:用於旋轉機械的健康監測。
計算複雜度較高
Gabor Transform在高維數據(如圖像信號處理)中,計算複雜度可能顯著增大,特別是在多頻段或多方向分析時。每個窗函數的計算都需要執行一次傅立葉變換,這對於大數據或實時應用場景來說,可能成為性能瓶頸。
計算複雜度的挑戰
在圖像處理中,Gabor Transform通常需要對圖像的不同尺度和方向應用一組 Gabor 濾波器,以提取豐富的特徵資訊。這意味着需要對每個尺度和方向分別進行濾波操作,導致計算量隨着濾波器數量的增加而線性增長。此外,對於高解像度圖像,每次濾波操作都涉及大量像素的計算,進一步增加了計算負擔。
多頻段和多方向分析的負擔
為了捕捉圖像中不同方向和頻率成分的特徵,通常需要使用多個 Gabor 濾波器組,涵蓋各種頻率和方向。這種多頻段和多方向的分析策略雖然能夠提供豐富的特徵表示,但也導致了計算資源的高需求。特別是在實時應用中,如視頻處理或即時圖像分析,這種高計算量可能導致延遲,影響系統的即時性和效率。
解像度折衷的不可避免性
根據測不準定理,Gabor Transform 的時間和頻率解像度達到了理論的最佳折衷,但這也意味着:
受測不凖定理約束,當需要同時對信號的快速變化與細微頻率差異進行精確分析時,時間和頻率的解像度會有可能不足以同時滿足所有需求。
相較於 Gabor Transform, Wigner Distribution Function(WDF)等方法,因是對訊號的自相關函數做傅立葉轉換,可以超越測不準原理約束的下限,因此能提供更高的時頻解像度,尤其是對於結構複雜的信號。然而,WDF 的非線性特性容易引入交叉干擾項(cross-terms),而為了為了結合兩者的優點,Gabor Wigner Transform應運而生
Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2024.
Alan V. Oppenheim, Ronald W. Schafer, John R. Buck : Discrete-Time Signal Processing, Prentice Hall, ISBN 0-13-754920-2
S. Qian and D. Chen, Joint Time-Frequency Analysis: Methods and Applications, Chap. 5, Prentice Hall, N.J., 1996.
Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2020.
S.C.Pei and S.G.Huang, STFT with adaptive window width based on the chirp rate . IEEE Transactions on Signal Processing, vol. 60,issue 8,pp. 4065-4080,2012.