加伯变换

加伯变换是窗函数为高斯函数的短时距傅立叶变换。

数学定义

将短时距傅立叶转换中的窗函数代入高斯函数，即可得下面的标准定义：

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