加伯变换 是窗函数 为高斯函数 的短时距傅立叶变换 。
将短时距傅立叶转换 中的窗函数 代入高斯函数 ,即可得下面的标准定义:
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.