短时距傅里叶变换

短时距傅里叶变换（Short-time Fourier Transform, STFT）是傅里叶变换的一种变形，也称作加窗傅里叶变换（Windowed Fourier transform）或Time-dependent Fourier transform，用于决定随时间变化的信号局部部分的正弦频率和相位。实际上，计算短时距傅里叶变换的过程是将长时间信号分成数个较短的等长信号，然后再分别计算每个较短段的傅里叶变换。通常拿来描绘频域与时域上的变化，为时频分析中其中一个重要的工具。

与傅里叶变换在概念上的区别

将信号做傅里叶变换后得到的结果，并不能给予关于信号频率随时间改变的任何信息。以下的例子作为说明：

$x(t)={\begin{cases}\cos(440\pi t);&t<0.5\\\cos(660\pi t);&0.5\leq t<1\\\cos(524\pi t);&t\geq 1\end{cases}}$

傅里叶变换后的频谱和短时距傅里叶变换后的结果如下：

傅里叶变换后, 横轴为频率(赫兹)

短时距傅里叶变换, 横轴为时间(秒)，纵轴为频率(赫兹)

由上图可发现，傅里叶变换只提供了有哪些频率成分的信息，却没有提供时间信息；而短时傅里叶变换则清楚的提供这两种信息。这种时频分析的方法有利于频率会随着时间改变的信号，如音乐信号和语音信号等分析。

定义

连续短时傅里叶变换

简单来说，在连续时间的例子，一个函数可以先乘上仅在一段时间不为零的窗函数再进行一维的傅里叶变换。再将这个窗函数沿着时间轴挪移，所得到一系列的傅里叶变换结果排开则成为二维表象。数学上，这样的操作可写为：

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

另外也可用角频率来表示：

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

其中$w(t)$是窗函数，窗函数种类有很多种，会在稍后再做仔细讨论。$x(t)$是待变换的信号。$X(t,\omega )$是$w(t-\tau )x(\tau )$的傅里叶变换。 随着$t$的改变，窗函数在时间轴上会有位移。经$w(t-\tau )x(\tau )$后，信号只留下了窗函数截取的部分做最后的傅里叶变换，所得到的结果为一复数函数，代表着信号随时间与频率变化的大小与相位。

离散短时傅里叶变换

在离散时间的例子，数据会被切割成数个大量的帧，而每组帧通常会互相重叠，避免因切割方式造成边界的误差。而每组帧在各自进行傅里叶变换后所得的复数结果会再进行相加，可得到每个点时间与频率变化的大小与相位。数学上，这样的操作可写为：

$\mathbf {STFT} \{x[n]\}(m,\omega )\equiv X(m,\omega )=\sum _{n=-\infty }^{\infty }x[n]w[n-m]e^{-j\omega n}$

相同地，其中$w[n]$是窗函数，$x[n]$是待变换的信号。在这个例子里，m是离散的且ω是连续的，但大部分实际的应用当中，短时距傅里叶变换在电脑中都是以快速傅里叶变换进行计算(见实现方法的快速傅里叶变换)，而此时这两个参数都是离散且被量化的。

Sliding 离散傅里叶变换

当只想要得知特定少数的ω，或是短时距傅里叶变换每次窗函数移动m的值，则短时距傅里叶变换可以利用sliding DFT算法更有效地计算出来。

反短时距傅里叶变换

短时距傅里叶变换是可逆的，也就是说原本的信号可以借由反短时距傅里叶变换将短时距傅里叶变换后的信号还原。

其中最广为接受的反短时距傅里叶变换方法是重叠-相加之卷积法，此方法也促成了更多样的信号处理方法。

反短时距傅里叶变换，其数学类似傅里叶变换，但须消除窗函数的作用，首先必须先将窗函数的总面积规模化使得

$\int _{-\infty }^{\infty }w(\tau )\,d\tau =1.$

而从上也可轻易地得出

$\int _{-\infty }^{\infty }w(t-\tau )\,d\tau =1\quad \forall \ t$

和

$x(t)=x(t)\int _{-\infty }^{\infty }w(t-\tau )\,d\tau =\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau .$

连续傅里叶变换公式如下：

$X(\omega )=\int _{-\infty }^{\infty }x(t)e^{-j\omega t}\,dt.$

将$x(t)$进行上述的替换：

$X(\omega )=\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau \right]\,e^{-j\omega t}\,dt$

$=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-j\omega t}\,d\tau \,dt.$

将积分顺序进行交换：

$X(\omega )=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-j\omega t}\,dt\,d\tau$

$=\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-j\omega t}\,dt\right]\,d\tau$

$=\int _{-\infty }^{\infty }X(\tau ,\omega )\,d\tau .$

因此傅里叶变换可以视为某种将$x(t)$所有的短时距傅里叶变换的相位同调部分进行相加。

而反傅里叶变换公式如下：

$x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\omega )e^{+j\omega t}\,d\omega ,$

因此 $x(t)$可以从$X(\tau ,\omega )$被复原

$x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+j\omega t}\,d\tau \,d\omega .$

或

$x(t)=\int _{-\infty }^{\infty }\left[{\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+j\omega t}\,d\omega \right]\,d\tau .$

与上面所列的窗函数的式子进行比较，可得

$x(t)w(t-\tau )={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+j\omega t}\,d\omega .$

对反傅里叶变换公式中的$X(\tau ,\omega )$来说$\tau$是不变的

$x(t)=w(t_{1}-t)^{-1}\int _{-\infty }^{\infty }X(t_{1},f)e^{j2\pi ft}\,df;\ \ w(t_{1}-t)\neq 0$

另外用角频率来表示：

$x(t)={\frac {1}{2\pi }}w^{-1}(t_{1}-t)\int \limits _{-\infty }^{\infty }X(t_{1},w)e^{jwt}dw$

窗函数

窗函数通常满足下列特性:

1. $w(t)=w(-t)\,$，即为偶函数。
2. $max(w(t))=w(0)\,$，即窗函数的中央通常是最大值的位置。
3. $w(t_{1})\geq w(t_{2}),|t_{2}|\geq |t_{1}|$，即窗函数的值由中央开始向两侧单调递减。
4. $w(t)\cong 0,|t|\to \infty$，即窗函数的值向两侧递减为零。

常见的窗函数有：方形、三角形、高斯函数等，而短时距傅里叶变换也因窗函数的不同而有不同的名称。而加伯变换，即为窗函数是高斯函数的短时距傅里叶变换，通常没有特别说明的短时距傅里叶变换，即为加伯变换。

非对称窗函数

当在特殊应用时，窗函数特性的第一点可以不满足，如下图的非对称窗函数$w(t)$，其中$B_{1}\neq B_{2}$。左图为窗函数原本的图形，而在计算短时距傅里叶变换时，需将窗函数转到$\tau$轴上得出$w(t-\tau )$，换言之，欲得到的短时距傅里叶变换的结果需在$t+B_{1}$的时间点才能算出，因此若$B_{1}$愈小，即可愈快得结果，此种非对称窗函数可应用在地震波、碰撞侦测...等，需要即时处理的应用。

优缺点

• 优点：比起傅里叶变换更能观察出信号瞬时频率的信息。
• 缺点：计算复杂度高

方形窗函数的短时距傅里叶变换

概念

方形窗函数，B = 50，横轴为时间(秒)

右图即为方形窗函数的一个例子，其数学定义： $w(t)={\begin{cases}\ 1;&|t|\leq B\\\ 0;&|t|>B\end{cases}}$

可以随要分析的信号，来调整B的大小（即调整方形窗函数的宽度）。至于B的选择，将会在下面探讨。

短时傅里叶变换可以简化为

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

反短时傅里叶变换可简化为

$x(t)=\int _{-\infty }^{\infty }X(t_{1},f)e^{j2\pi ft}\,df;t-B<t_{1}<t+B$

特性

其大部分的特性都与傅里叶变换的特性相对应

• 积分特性

$\int _{-\infty }^{\infty }X(t,f)\,df=\int _{t-B}^{t+B}x(\tau )\int _{-\infty }^{\infty }e^{-j2\pi f\tau }\,df\,d\tau ={\begin{cases}\ x(0);&|t|\leq B\\\ 0;&|t|>B\end{cases}}$

• 位移特性(时间轴方向的移动)

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

• 调变特性(频率轴方向的移动)

$\int _{t-B}^{t+B}\left(x(\tau )e^{j2\pi f_{0}\tau }\right)e^{-j2\pi f\tau }\,d\tau =X(t,f-f_{0})$

• 线性特性

若有一信号$h(t)=\alpha x(t)+\beta y(t)\,$，$H(t,f),X(t,f),Y(t,f)\,$分别为$h(t),x(t),y(t)\,$做方形窗函数短时 距傅里叶变换的结果，则$H(t,f)=\alpha X(t,f)+\beta Y(t,f)\,$。

• 能量积分特性

$\int _{-\infty }^{\infty }|X(t,f)|^{2}\,df=\int _{t-B}^{t+B}|x(\tau )|^{2}\,d\tau$

$\int _{-\infty }^{\infty }X(t,f)Y^{*}(t,f)\,df=\int _{t-B}^{t+B}x(\tau )y^{*}(\tau )\,d\tau$

• 特殊信号

1. 当$x(t)=\delta (t)\,$，

$X(t,f)={\begin{cases}\ 1;&|t|\leq B\\\ 0;&|t|>B\end{cases}}$

2. 当$x(t)=1\,$，

$X(t,f)=2Bsinc(2Bf)e^{-j2\pi ft}\,$

方形窗函数宽度$(B)$的选取

方形窗函数短时距傅里叶变换用不同窗函数宽度(B)的比较，横轴为时间(秒)，纵轴为频率(赫兹)

• 由上述特性中的特殊信号$x(t)=\delta (t)$来分析，信号只有在$t=0$的时候有值；若短时距傅里叶变换是理想的话，$X(t,f)$应该只有在$t=0$的时候有能量。但由上面的特性可发现，能量会出现在${\begin{smallmatrix}|t|\leq B\end{smallmatrix}}$中间。因此，若我们取较小的$B$，则可使结果趋近理想。
• 接着我们来分析$x(t)=1$，信号因为没有改变，应该为DC。若短时距傅里叶变换是理想的话，$X(t,f)$应该只有在$f=0$的时候有能量。但由上面的特性可发现，能量会沿着频率轴呈现sinc函数。若我们取较大的$B$，可使sinc函数沿着频率轴变窄，使得结果趋近理想。
• 综合以上说明，若我们使用较大的方形窗函数宽度$(B)$，则$X(t,f)$时间轴的分辨率会下降；频率轴的分辨率上升。若使用较小的$B$，则$X(t,f)$时间轴的分辨率会上升；频率轴的分辨率下降。我们以下面做为例子说明：

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

结果如右图所示，B越大则在频率变化处(t = 10, 20)附近的频率越不准确，即可能会有多个频率成分出现。但同时，其他时间点的能量则较集中；没有如B较小时，频率散开或模糊的情形。

上述也是其中一个小波变换及多分辨率分析作为改进的方向，其中多分辨率分析能在高频时有较好的时间轴解析，而在低频时能有较好的频率轴解析，此种组合较契合许多实际的应用。

时间轴与频率轴的分辨率无法同时提升也与海森堡不确定性原理有关，即时间与频率的标准差乘积有所限制，而高斯函数恰好能符合不确定性原理的极值，也就是两者同时达到最好的分辨率，而应用高斯函数的时频分析方法即为加伯变换，而在经过修改及多分辨率分析后，成为了莫莱小波。

优缺点

• 优点：方形窗函数的短时距傅里叶变换有许多可应用的数学特性，在数字的应用上所需的计算时间较少。
• 缺点：时频分析的表现较差

其他窗函数

高斯窗函数

概念

高斯窗函数的短时距傅里叶变换又称为加伯变换。以下是高斯函数的数学定义，

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

据此，短时傅里叶变换可以写为

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

优缺点

• 优点：可以在时间跟频率上有更好的平衡，得到较清楚的时频图。
• 缺点：因窗函数跟信号本身的乘法，计算时间跟复杂度都比较高。

三角形窗函数

三角形函数，横轴为时间，B=1/2

概念

三角形窗函数如右图所示，数学定义如下，

$w(t)=max(1-\left\vert t\right\vert ,0)$

$w(t)={\begin{cases}1-\left\vert t\right\vert ,&\left\vert t\right\vert <2B;\\0,&{\text{otherwise.}}\end{cases}}$

可使用在震幅改变的情况下，相对于方形窗函数，可更好的滤除噪声。

海宁(Hanning/ Hann)窗函数

海宁函数

概念

海宁函数如右图所示，数学定义如下，

$w(t)={\begin{cases}0.5+0.5cos(\pi t/B),&{\text{when }}\left\vert t\right\vert \leq B\\0,&{\text{otherwise }}\end{cases}}$

相较于三角形窗函数，海宁窗函数更为贴近现实信号的趋势，可进一步滤除噪声。

汉明(Hamming)窗函数

汉明函数

概念

汉明窗函如右图所示，数学定义如下，

$w(t)={\begin{cases}0.54+0.46cos(\pi t/B),&{\text{when }}\left\vert t\right\vert \leq B\\0,&{\text{otherwise }}\end{cases}}$

跟海宁窗函数类似，但两端不为零。

海宁与汉明窗的区别

窗函数有四个指标，分别为

• 泄露指数 (Leakage Factor)
• 主办宽度 (Mainlobe width)
• 旁办衰减 (Sidelobe attenuation)
• 旁办滚降率 (Sidelobe roll-off rate)

  

  方形窗函数宽度(B)与STFT清晰率的取舍，横轴为时间(秒)，纵轴为频率(赫兹)

因为汉明窗两端不能到零，而海宁窗两端为零。从以上频率响应来看，汉明窗可以有效减少靠近的旁办，但在较远的旁办泄漏比海宁窗严重。

如何决定窗函数

可根据以下条件来选取窗函数，

• 复杂度，方形复杂度较低
• 解析率，以方形为例，越宽的主办可以得到更清楚的时频图，却会把噪声也一同显示，反之则得到不清晰的时频图

在决定复杂度跟解析率后，可利用不同的窗函数达到更好的滤噪声效果。

瑞利频率

当Nyquist频率是能被有意义分析的频率最大值的限制，而瑞利频率则是能被有限带宽频的窗函数解析的频率最小值的限制。若给定一窗函数的长度是T秒，最低能被解析的频率即为1/T Hz。

瑞利频率在短时距傅里叶变化的应用中扮演重要的角色，像是在分析神经信号时。

频谱(Spectrogram)

Spectrogram即短时傅里叶变换后结果的绝对值平方，两者本质上是相同的，在文献上也常出现spectrogram这个名词。

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

应用

应用短时距傅里叶变换分析声音信号

短时距傅里叶变换及其他工具经常用于分析音乐。

如右图所示，

1. 水平轴为频率，左侧为最低频率，右侧为最高频率
2. 条形高度(混和颜色表示)表示该频带内的频率幅度
3. 深度表示时间

音频工程师使用这种视觉来获取有关音频样本的信息。

此外，因频率会随时间而改变，短时距也可使用在以下情境，

• 信号采样 (signal sampling)，
• 调变 (modulation)，
• 生物信号 (biomedical signals)，等等

若与时间无关，如卷积，照片等则不能使用短时距傅里叶变换来进行分析。而影片属于3D信号，其短时距傅里叶产物为6D信号，故也不适用。

短时距傅里叶变换实现方法

从连续短时距傅里叶变化的定义出发

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

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

若当$\left|t\right|>B,w(t)\cong 0\qquad {\frac {B}{\Delta _{t}}}=Q$

上式可改写为

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

---

直接运算

限制条件

(1)要满足Nyquist criterion

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

$x(\tau )$的带宽为$\Omega _{x}$。而$w(\tau )$的带宽则为$\Omega _{w}$，$w(t-\tau )$的带宽也为$\Omega _{w}$

因为在时域相乘相当于在频域做卷积，因此$x(\tau )w(t-\tau )$的带宽为$\Omega _{x}+\Omega _{w}$(通常$\Omega _{x}$会远大于$\Omega _{w}$，所以主要影响带宽的是$\Omega _{x}$)

推导

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

变换到离散形式($t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{t}$)，其中$\Delta _{t}={\frac {1}{f_{s}}}$

$X(n\Delta _{t},m\Delta _{f})=\sum _{p=-\infty }^{\infty }w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}$，由于无限大的上下限实务上做不到，所以尝试变成有限大的上下限。

假设$w(t)\cong 0$ for $|t|>B,{\frac {B}{\Delta _{t}}}=Q$

$X(n\Delta _{t},m\Delta _{f})=\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}$

• 对于缩放的加伯变换，$Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}$

时间复杂度

$TF(2Q+1)\to O(TFQ)$

假设t-axis有T个采样点，f-axis有F个采样点，则我们总共要对TF个点做$(2Q+1)$次的运算，因此可得复杂度为$TF(2Q+1)$

优缺点

优点：简单及有弹性(因为限制少)

缺点：复杂度较高

---

快速傅里叶变换

限制条件

(1)要满足Nyquist criterion

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

(2)${\Delta _{t}}{\Delta _{f}}={\textstyle {1 \over {N}}}$ (N可为任意整数)

(3) $N\geq 2Q+1$ (做N点傅里叶变换，输入必要<=N)

推导

标准的离散傅里叶变换式子为

$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)=\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$代入上式即可得

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

其中 ${\begin{cases}{x_{1}}\left(q\right)=w\left({(Q-q){\Delta _{t}}}\right)x\left({(n-Q+q){\Delta _{t}}}\right),&{\mbox{for}}{\rm {0}}\leq q\leq 2{\rm {Q}}\\{x_{1}}\left(q\right)=0,&{\mbox{for}}{\rm {2}}Q<q\leq N\end{cases}}$

运算步骤

假设$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}$

步骤一：计算$n_{0},m_{0},T,F,N,Q$

步骤二：$n=n_{0}$

步骤三：决定$x_{1}(q)$

步骤四：$X_{1}(m)=FFT[x_{1}(q)]$

步骤五：变换$X_{1}(m)$成$X(n\Delta _{t},m\Delta _{f})$

步骤六：设$n=n+1$，并回到步骤三，直到$n=n_{0}+T+1$

• 范例

${\begin{cases}x(t)=\cos {(2\pi t)},&{\mbox{when }}t<10\\x(t)=\cos {(6\pi t)},&{\mbox{when }}10\leq t<20\\x(t)=\cos {(4\pi t)},&{\mbox{when }}t\geq 20\\\end{cases}}$

借由采样定理可得知$\Delta _{t}<{\frac {1}{6}}$

假设$f=-5\sim 5$及$\Delta _{f}=0.1$，则经由$f=m\Delta _{f}$可得$m=-50\sim 50$

$\;t=0\sim 30$及$\Delta _{t}=0.1$，则经由$t=n\Delta _{t}$可得$n=0\sim 300$

步骤一：$n_{0}=0,m_{0}=-50,T=301,F=101,N={\frac {1}{\Delta _{t}\Delta _{f}}}=100,Q={\frac {B}{\Delta _{t}}}=10$

步骤二：$n=n_{0}=0$

步骤三：计算$x_{1}(q)(q=0\sim 99)$

步骤四：利用求得的$x_{1}(q)$计算快速傅里叶变换 $X_{1}[m]=\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}}$

步骤五：变换$X_{1}(m)$到$X(n\Delta _{t},m\Delta _{f})$

$X(n\Delta _{t},m\Delta _{f})=X_{1}[m]\Delta _{t}e^{j{\frac {2\pi (Q-n)m}{N}}}$

• 注：若是于程式中执行，要注意m可能为负数，所以需要利用到周期性性质$X_{1}[m]=X_{1}[m+N]$

$X_{1}[-50]=X_{1}[50],X_{1}[-49]=X_{1}[51],\cdots \cdots ,X_{1}[-1]=X_{1}[99]$

因此可将上式改为$X(n\Delta _{t},m\Delta _{f})=X[((m))_{N}]e^{j{\frac {2\pi (Q-n)m}{N}}}$，其中$((m))_{N}$代表取m除以N的余数

步骤六：设定$n=n+1$，回到步骤三直到$n=n_{0}+T-1$

时间复杂度

利用FFT计算$\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}$，其中每次FFT的时间复杂度为 $N{\log _{2}}N$

总时间复杂度为$TN{\log _{2}}N\to O(TN{\log _{2}}N)$

优缺点

优点：与直接运算相比，复杂度较低

缺点：较多限制，包括 ${\begin{cases}\Delta _{t}<{\frac {1}{2(\Omega _{x}+\Omega _{w})}}\\\Delta _{t}\Delta _{f}={\frac {1}{N}}\\N\geq 2Q+1\\\end{cases}}$

---

使用快速傅里叶变换加上递回关系式

限制条件

(1)要满足Nyquist criterion

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

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

(3)$N\geq 2Q+1$

(4)需为方形窗函数的短时距傅里叶变换

推导

因为是方形窗函数 ${w}\left((n-p){\Delta _{t}}\right)=1$，因此原式可由此关系变成以下式子

${X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{{x}\left({p{\Delta _{t}}}\right)}{{w}\left((n-p){\Delta _{t}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}\to {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{{x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}$

而由此可看出n和n-1有递回关系，如下

${X}\left({(n-1){\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-1-Q}^{n-1+Q}{{x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}={X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)+x((n-Q-1)\Delta _{t})-x((n+Q+1)\Delta _{t})$

(1)以FFT计算${X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{j{\textstyle {{2\pi \,(Q-n_{0})m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}\qquad {n_{0}}=min(n)$

其中${\begin{cases}{x_{1}}\left(q\right)=x\left({(n-Q+q){\Delta _{t}}}\right),&{\mbox{for}}{\rm {0}}\leq q\leq 2{\rm {Q}}\\{x_{1}}\left(q\right)=0,&{\mbox{for}}{q>2{\rm {Q}}}\end{cases}}$

(2)利用递回关系式计算算${X}\left({{n}{\Delta _{t}},m{\Delta _{f}}}\right),\qquad n={n_{0}}+1\backsim max(n)$

则${X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={X}\left({(n-1){\Delta _{t}},m{\Delta _{f}}}\right)-{x}\left((n-Q-1){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n-Q-1)m} \over N}}}}{\Delta _{t}}+{x}\left((n+Q){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n+Q)m} \over N}}}}{\Delta _{t}}$

时间复杂度

(1)FFT计算一次 ${X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{j{\textstyle {{2\pi \,(Q-n_{0})m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}\qquad {n_{0}}=min(n)$

• 时间复杂度：$O(N\log _{2}N)$

(2)利用递回关系，计算$n=n_{0}+1$时的数值，因此共会执行T-1次递回，如下式

${X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={X}\left({(n-1){\Delta _{t}},m{\Delta _{f}}}\right)-{x}\left((n-Q-1){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n-Q-1)m} \over N}}}}{\Delta _{t}}+{x}\left((n+Q){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n+Q)m} \over N}}}}{\Delta _{t}}$

每次递回都要计算${x}\left((n-Q-1){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n-Q-1)m} \over N}}}}{\Delta _{t}}$及${x}\left((n+Q){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n+Q)m} \over N}}}}{\Delta _{t}}$两个乘法(相当于2F的复杂度)

• 时间复杂度：$2F(T+1)\to O(TF)$

总时间复杂度 $2(T-1)F+N{\log _{2}}N\to O(FT)$

优缺点

优点：四种运算中，最低的复杂度$O(TF)$

缺点：

1. 只适用于方形窗函数的短时傅里叶变换
2. 由于递回的关系，会有累加误差。所以只要当中有小错误，误差会累积到最后，造成无可预期的错误
3. 不能用在不平衡的采样点

---

使用Chirp-Z 变换

限制条件

(1)要满足Nyquist criterion

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

推导

令$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}})$

即可由直接运算的式子导出Chirp_Z变换的式子，如下所示

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

时间复杂度

当n为定值时

(1)假设$x_{1}[p]=w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}\to$ 相乘时间复杂度为2Q+1

(2)令$c[m]=e^{j\pi m^{2}\Delta _{t}\Delta _{f}}$，则$\sum \limits _{p=n-Q}^{n+Q}x_{1}[p]c[m-p]\to$ convolution时间复杂度为 $3N{\log _{2}}N$

(3)${\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}}}}\to$相乘时间复杂度为 F

因此，总时间复杂度为$T(2Q+1+F+3N{\log _{2}}N)\to O(TN{\log _{2}}N)$

虽然此实现方法和使用FFT计算的时间复杂度相同，但因为convolution相当于做三次FFT，因此实际操作时运算时间约为使用FFT计算的2~3倍

优缺点

优点：只有一项限制：$\Delta _{t}<{\frac {1}{2(\Omega _{x}+\Omega _{w})}}$

缺点：与前四种相比，复杂度是中间的。

---

Unbalanced Sampling for STFT and WDF

将直接法和快速傅里叶变换方法做修正

1.直接法

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

修正后 ：$X(n\Delta _{t},m\Delta _{f})=\sum _{p=nS-Q}^{nS+Q}w((nS-p)\Delta _{\tau })x(p\Delta _{\tau })e^{-j2\pi pm\Delta _{\tau }\Delta _{f}}\Delta _{\tau }$

其中， $t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{\tau },B=Q\Delta _{\tau }$ ,$S={\frac {\Delta _{t}}{\Delta _{\tau }}},\Delta _{t}\neq \Delta _{\tau }$

假设$w(t)\approxeq 0$ for $|t|>B$，则上下限可借由以下推导而修正

$\int _{t+B}^{t-B}\to \int _{n\Delta _{t}+Q\Delta _{\tau }}^{n\Delta _{t}-Q\Delta _{\tau }}$ 则上限可以写成$n\Delta _{t}+Q\Delta _{\tau }==nS\Delta _{\tau }+Q\Delta _{\tau }=\Delta _{\tau }(nS+Q)$，下限则以此类推

注：$\Delta _{\tau }$(输入信号的采样间隔)

$\Delta _{t}$(在t轴上的输出信号的采样间隔)

然而，$S={\frac {\Delta _{t}}{\Delta _{\tau }}}$是整数会是比较好的。

• 假设一声音信号：

${\begin{cases}\Delta _{\tau }={\frac {1}{44100}}\\\Delta _{t}={\frac {1}{100}}\end{cases}}$ 则经由上述公式可求得S=441，代表经由unbalanced sampling，我们跟原本$\Delta _{t}=\Delta _{\tau }={\frac {1}{44100}}$相比可减少441倍的采样点。

时间复杂度

由于t轴的采样点少了S倍，因此跟原本的直接运算复杂度相比，只要把$T\to {\frac {T}{S}}$即可，如下：

复杂度：$O({\frac {T}{S}}N\log _{2}N)$

2.快速傅里叶变换

限制条件

(1) $\Delta _{\tau }\Delta _{f}={\frac {1}{N}}$

(2) $N={\frac {1}{\Delta _{\tau }\Delta _{f}}}>2Q+1$ : ($\Delta _{\tau }\Delta _{f}$只要是整数的倒数即可)

(3) $\Delta _{\tau }<{\frac {1}{2\Omega }}$，$w(\tau -t)x(\tau )$的带宽是 $\Omega$

i.e. $|FT\{w(\tau -t)x(\tau )\}|=|X(t,f)|\approx 0$ ，当 $|f|>\Omega$

过程

$X(n\Delta _{t},m\Delta _{f})=\sum _{p=nS-Q}^{nS+Q}w((nS-p)\Delta _{\tau })x(p\Delta _{\tau })e^{-j{\tfrac {2\pi pm}{N}}}\Delta _{\tau }$

令 $q=p-(nS-Q)\longrightarrow p=(nS-Q)+q$

$x_{1}(q)=w((Q-q)\Delta _{\tau })x((nS-Q+q)\Delta _{\tau })$ for $0\leq q\leq 2Q$

$x_{1}(q)=0\qquad \qquad \qquad \qquad \qquad \qquad \quad \quad$for $2Q<q<N$

修正后：$X(n\Delta _{t},m\Delta _{f})=\Delta _{\tau }e^{j{\tfrac {2\pi (Q-nS)m}{N}}}\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\tfrac {2\pi qm}{N}}}$

运算步骤

假设$t=c_{0}\Delta _{t},(c_{0}+1)\Delta _{t},\cdots ,(c_{0}+C-1)\Delta _{t}=c_{0}S\Delta _{\tau },(c_{0}S+S)\Delta _{\tau },\cdots ,[c_{0}S+(C-1)S]\Delta _{\tau }$

$\quad \;\;f=m_{0}\Delta _{f},(m_{0}+1)\Delta _{f},\cdots ,(m_{0}+F-1)\Delta _{f}$

$\quad \;\;\tau =n_{0}\Delta _{\tau },(n_{0}+1)\Delta _{\tau },\cdots ,(n_{0}+T-1)\Delta _{\tau }$

步骤一：计算$c_{0},m_{0},n_{0},C,F,T,N,Q$

步骤二：$n=c_{0}$

步骤三：决定$x_{1}(q)$

步骤四：$X_{1}(m)=FFT[x_{1}(q)]$

步骤五：变换$X_{1}(m)\to X(n\Delta _{t},m\Delta _{f})$

步骤六：设定$n=n+1$及返回步骤三，直到$n=c_{0}+C-1$

复杂度

$O({\frac {T}{S}}N\log _{2}N)$

Non-Uniform $\Delta _{t}$

(1) 先用比较大的$\Delta _{t}$

(2) 如果发现$|X(n\Delta _{t},m\Delta _{f})|$ 和 $|X((n+1)\Delta _{t},m\Delta _{f})|$ 之间有很大的差异，则在$n\Delta _{t}$，$(n+1)\Delta _{t}$ 之间选用比较小的采样区间$\Delta _{t1}$

($\Delta _{\tau }<\Delta _{t1}<\Delta _{t}$，${\frac {\Delta _{t}}{\Delta _{t1}}}$ 和 ${\frac {\Delta _{t1}}{\Delta _{\tau }}}$皆为整数)

再用Unbalanced Sampling for STFT and WDF 中修正后的快速傅里叶变换方法算出 $X(n\Delta _{t}+\Delta _{t1},m\Delta _{f})$，$X(n\Delta _{t}+2\Delta _{t1},m\Delta _{f})$$X((n+1)\Delta _{t}-\Delta _{t1},m\Delta _{f})$

(3) 以此类推，如果 $|X(n\Delta _{t}+k\Delta _{t1},m\Delta _{f})|,|X((n+1)\Delta _{t}+(k+1)\Delta _{t1},m\Delta _{f})|$的差距还是太大，则再选用更小的采样间隔$\Delta _{t2}$

($\Delta _{\tau }<\Delta _{t2}<\Delta _{t1}$，${\frac {\Delta _{t1}}{\Delta _{t2}}}$ 和 ${\frac {\Delta _{t2}}{\Delta _{\tau }}}$皆为整数)

• 比较

若有一音乐信号总共有1.6秒，$\Delta _{\tau }={\frac {1}{44100}}$

1. 选择$\Delta _{t}=\Delta _{\tau }$，则共有$44100*1.6+1=70561$点
2. 选择$\Delta _{t}=0.01=441\Delta _{\tau }$，则共有$100*1.6+1=161$点
3. t随时间不同有不同的选择，如下

$t=0,0.05,0.1,0.15,0.2,0.4,0.45,0.46,0.47,0.48,0.49,0.5,0.55,0.6,0.8,0.85,0.9,0.95,0.96,0.97,0.98,0.99,1,1.05,1.1,1.15,1.2,1.4,1.6$，共29点

可以这样做的原因为：有些音乐信号在和弦与和弦中间几乎没有变化，因此可以挑选较大的$\Delta _{t}$采样;和弦在变换时，频率会变化的较剧烈，因此变换和弦是需要用较多的采样点。借由此种non-uniform的采样，可以让我们大幅减少运算量，从最一开始的$70561\to 29$可看出我们的运算量大幅降低。

参见

• 闵可夫斯基空间
• 柯西不等式
• 三角不等式
• 完备空间

参考书目、数据来源

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. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2017.
4. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2020.
5. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2022.