加伯轉換

加伯轉換是窗函數為高斯函數的短時距傅立葉變換。

數學定義

將短時距傅立葉轉換中的窗函數代入高斯函數，即可得下面的標準定義：

$G_{x}(t,f)=\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )\,d\tau$

以下是幾種常見的替代定義：

1. $Gx_{1}(t,f)=\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f(\tau -{\frac {t}{2}})}x(\tau )\,d\tau$
2. $Gx_{2}(t,f)={\sqrt[{4}]{2}}\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )\,d\tau$
3. $Gx_{3}(t,\omega )=\int _{-\infty }^{\infty }e^{-{\frac {(\tau -t)^{2}}{2}}}e^{-j\omega \tau }x(\tau )\,d\tau$
4. $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$

• 註：在文獻上可能會看到不同形式的加伯轉換，但本質上都是一樣的。

由於實作時，不能計算無限大的積分式子，所以根據高斯函數會從兩側遞減的性質，我們可以將上式進一步化簡：

$\because {\begin{cases}\ e^{-\pi a^{2}}<0.00001;&|a|>1.9143\\\ e^{-a^{2}/2}<0.00001;&|a|>4.7985\end{cases}}$

$\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}}$

為何選擇高斯函數作為窗函數

1. 其他窗函數的短時距傅立葉變換，如利用方型窗函數的短時距傅立葉變換，無法同時兼顧時間軸和頻率軸的解析度；一者解析度提升，另一者解析度必定下降。但高斯函數由海森堡測不準原理可得知，是最能同時讓兩軸兼顧解析度的窗函數(將於下面章節詳述)。
2. 高斯函數為傅立葉轉換的特徵函數：

• $\int _{-\infty }^{\infty }e^{-\pi t^{2}}e^{-j2\pi f\tau }dt=e^{-\pi f^{2}}$
• $\int _{-\infty }^{\infty }e^{-{\frac {t^{2}}{2}}}e^{-j\omega t}dt=e^{-{\frac {f^{2}}{2}}}$

因此經過轉換後其性質不變。因此可讓加伯轉換後在時間軸和頻率軸的性質相互對稱。

由測不準原理了解高斯函數的性質

上述提到，高斯函數是最能兼顧時間與頻率解析度的窗函數。我們利用這個章節來詳細討論。

• 海森堡測不準原理:

對於一個信號 $x(t)$，當$|t|\to \infty$，若${\sqrt {t}}x(t)=0$，則

$\sigma _{t}\sigma _{f}\geq 1/4\pi \,$

其中 $\sigma _{t}^{2}=\int (t-\mu _{t})^{2}P_{x}(t)dt,\sigma _{f}^{2}=\int (f-\mu _{f})^{2}P_{X}(f)df\,$

$\qquad \mu _{t}=\int tP_{x}(t)dt\quad \quad \quad \ ,\mu _{f}=\int fP_{X}(f)df$

$\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)$為高斯函數時

$x(t)=e^{\pi t^{2}},X(f)=e^{-\pi f^{2}}$

套用以上函式求得變異數(其中由於高斯函數為偶對稱函數，所以其$\mu _{t}=0$)

$\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 }}$

藉由微積分公式可得：$\int _{-\infty }^{\infty }|x(t)|^{2}dt=\int _{-\infty }^{\infty }e^{-2\pi t}dt={\frac {1}{\sqrt {2}}}$

及$\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}}}}$

$\therefore \sigma _{t}\sigma _{f}={\frac {1}{4\pi }}$

即高斯函數滿足測不準定理的最下限，所以是所有窗函數中能使時間和頻率兩者解析度都達到最高的函數。

• 變形的高斯函數同樣會滿足測不準原理的下限，如以下例子：

1. $e^{-\pi (t-t_{0})^{2}}$：對機率分布做位移，標準差不會改變。
2. $Ae^{-\pi t^{2}}$：分子與分母同乘A，可消掉。因此標準差不會改變。
3. $e^{-\pi t^{2}}e^{j2\pi f_{0}t}$：在時域乘上$e^{j2\pi f_{0}t}$相當於在頻域對頻率做位移，標準差一樣不會改變。
4. $e^{-\pi bt^{2}}$：在時域做縮放，頻域會做相反的縮放，因此標準差也不會改變。

• 以下提供一個簡單的例子來做模擬，

$x(t)={\begin{cases}\cos(2\pi t);&t<10\\\cos(6\pi t);&10\leq t<20\end{cases}}$

右圖為即加伯轉換的結果，可以看出其時間和頻率都維持相當程度的解析度。

高斯窗函數與方形窗函數比較

以下提供一個簡單的範例來比較加伯轉換以及利用方形窗函數的短時傅立葉轉換：

$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)={\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)\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^{-\pi \sigma t^{2}}$，則傅立葉轉換後為$W(f)={\frac {1}{\sqrt {\sigma }}}e^{-{\frac {\pi f^{2}}{\sigma }}}$，使其能根據需求而調整時域解析度或頻域解析度
• 改變高斯函數的寬度，和改變方形窗函數短時距傅立葉變換的效果類似。若選取較大的$\sigma$，時域的高斯窗函數較窄，則時域有較高的解析度，而頻域的高斯窗函數較寬，所以頻域的解析度會下降(通常用於需要時域解析度較高的應用，例如：音樂訊號);反之，若選取較小的$\sigma$，時域的高斯窗函數較寬，則時域的解析度下降，而頻域的高斯窗函數較窄，所以頻域的解析度會上升(通常運用在需要頻域解析度較高的應用，例如：氣候)。雖然還是有兩軸之間的解析度的犧牲，但比起其他無法滿足測不準原理下限的窗函數，加伯轉換的兩軸還是能相對維持較高的解析度。
• 若應用於瞬時頻率改變較劇烈的應用，則可考慮使用窗寬度隨時間而變動的加伯轉換數學式子，如下

$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$

當瞬時頻率變動非常快時，使用較大的$\sigma$值，使其時域解析度能較高;當瞬時頻率變動很慢時，使用較小的$\sigma$值，使其頻域解析度能較高。

實現方法及注意事項

Direct Implementation

$X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau$

$w(t)=e^{-\pi \sigma t^{2}}$

• Discrete Form:

令$t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{t}$ 可將式子改寫為離散形式:

${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)\cong 0\qquad for\left|t\right|>B,{\frac {B}{\Delta _{t}}}=Q$ $w((n-p)\Delta _{t})\cong 0\qquad$ $for\left|n-p\right|>{\frac {B}{\Delta _{t}}}$ ,$\left|p-n\right|>Q$ therefore,only when $-Q<p-n<Q$ $w((n-p)\Delta _{t})$ is nonzero 可改寫為:

${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^{-\pi \sigma a^{2}}<0.00001$ $when\left|a\right|>1.9143$ $Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}$ $B={\frac {1.9143}{\sqrt {\sigma }}}$

限制

• 避免贋頻效應(aliasing effect)

(1)${\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}$

時間複雜度

O(TFQ) T:時間取樣點數 F:頻率取樣點數 Q:$Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}$

優缺點

優點:簡單實現,限制條件少

缺點:時間複雜度高

FFT-Based Method(快速傅立葉轉換)

由Direct Implementation可得下式

${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)\to p=(n-Q)+q$且離散傅立葉轉換標準式$Y[m]=\sum \limits _{n=0}^{N-1}y[n]e^{-j{\frac {2\pi mn}{N}}}$ 可將式子整理為:

${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}}$以fft()算出帶入即可實現

其中${x_{1}}\left(q\right)=w\left({(Q-q){\Delta _{t}}}\right)x\left({(n-Q+q){\Delta _{t}}}\right)$ ,$0\leq q\leq 2Q$,$w(t)=e^{-\pi \sigma t^{2}}$

${x_{1}}\left(q\right)=0,2Q<q\leq N$

$Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}$

$B={\frac {1.9143}{\sqrt {\sigma }}}$

• Matlab及python 皆可呼叫fft函式完成$Y[m]=\sum \limits _{n=0}^{N-1}y[n]e^{-j{\frac {2\pi mn}{N}}}$

• 演算法

假設$t=n_{0}\Delta _{t},(n_{0}+1)\Delta _{t},\cdots \cdots ,(n_{0}+T-1)\Delta _{t}$

$\,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$

step 2：$n=n_{0}$

step 3：決定$x_{1}(q)$

step 4：$X_{1}(m)=FFT[x_{1}(q)]$

step 5：轉換$X_{1}(m)$成$X(n\Delta _{t},m\Delta _{f})$

step 6：設$n=n+1$ and return to Step 3 until $n=n_{0}+T+1$

限制

• 避免贋頻效應(aliasing effect)

(1)${\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}$(基本上任何實現方法都要避免贋頻效應)

(2)${\Delta _{t}}{\Delta _{f}}={\textstyle {1 \over {N}}}$

(3)$N=1/{\Delta _{t}}{\Delta _{f}}\geq 2Q+1$

時間複雜度

$O(TN{\log _{2}}N)$

優缺點

優點:時間複雜度低

缺點:限制條件較直接實現法多

Chirp Z Transform

可改寫為: 由Direct Implementation可得下式

${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^{-\pi \sigma a^{2}}<0.00001$ $when\left|a\right|>1.9143$ $Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}$ $B={\frac {1.9143}{\sqrt {\sigma }}}$

令$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}\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)\Delta _{t})x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}$ $\quad \quad n-Q\leq p\leq n+Q$

Step2:$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\Delta _{t},m\Delta _{f})=\Delta _{t}e^{-j\pi m^{2}\Delta _{t}\Delta _{f}}X_{2}[m,n]$

限制

• 避免贋頻效應(aliasing effect)

(1)${\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}$

時間複雜度

$O(TN{\log _{2}}N)$

優缺點

優點:限制條件與Direct Implementation法一樣基本上沒有限制

缺點:時間複雜度與FFT-Based Method(快速傅立葉轉換)一樣

• 但由於加伯轉換無法使用Recursive Method(遞迴法)所以此不能算是缺點

特性

加伯轉換的大部分的特性和方形窗函數短時距傅立葉轉換的特性都相似，有些特性甚至更加接近傅立葉轉換的特性。

• 積分特性

當$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,\int _{-\infty }^{\infty }G_{x}(t,f)\,df=e^{-\pi t^{2}}x(0)$

當$k=1,\int _{-\infty }^{\infty }G_{x}(t,f)e^{j2\pi tf}\,df=x(t)$ (還原成原始信號)

• 位移特性

若$y(t)=x(t-t_{0})\,$，則$G_{y}(t,f)=G_{x}(t-t_{0},f)e^{-j2\pi ft_{0}}\,$

• 調變特性

若$y(t)=x(t)e^{j2\pi f_{0}t}\,$，則$G_{y}(t,f)=G_{x}(t,f-f_{0})\,$

• 線性特性

若有一信號$h(t)=\alpha x(t)+\beta y(t)\,$，$G_{h}(t,f),G_{x}(t,f),G_{y}(t,f)\,$分別為$h(t),x(t),y(t)\,$做加伯轉換的結果，則$G_{h}(t,f)=\alpha G_{x}(t,f)+\beta G_{y}(t,f)\,$。

• 功率衰減特性

若$t>t_{0}$時$x(t)=0$，則

$\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$

• 能量積分特性

$\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$

$\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }G_{x}(t,f)G_{y}^{*}(t,f)\,dfdt=\int _{-\infty }^{\infty }x(\tau )y^{*}(\tau )\,d\tau$

• 能量總和特性

$\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)=\delta (t)\,$，$G_{x}(t,f)=e^{-\pi t^{2}}$

2. 當$x(t)=1\,$，$G_{x}(t,f)=e^{j2\pi ft}e^{\pi f^{2}}\,$

和方形窗函數短時距傅立葉轉換不同的是，加伯轉換的結果對於時間和頻率軸較對稱，也比較沒有旁波(sidelobe)；也印證了上述所說的，加伯轉換較能維持兩個軸的解析度。

優缺點

• 優點: 時頻圖較清晰
• 缺點: 計算複雜度比方形窗函數短時傅立葉變換來的高，因為需做窗函數內與信號本身的乘法。

參見

• 閔可夫斯基空間
• 柯西不等式
• 三角不等式
• 完備空間

參考書目、資料來源

1. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2011.
2. Alan V. Oppenheim, Ronald W. Schafer, John R. Buck : Discrete-Time Signal Processing, Prentice Hall, ISBN 0-13-754920-2
3. S. Qian and D. Chen, Joint Time-Frequency Analysis: Methods and Applications, Chap. 5, Prentice Hall, N.J., 1996.
4. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2020.
5. 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.