离散之方波短时距傅立叶变换

短时距傅立叶变换(英文：short-time Fourier transform, STFT)是和傅立叶变换相关的一种数学转换关系，用于时间和频域之间的分析。 简单来说，在连续时间的例子中，一个函数可以先乘上仅在一段时间不为零的窗函数(window function)再进行一维的傅立叶变换。再将这个窗函数沿着时间轴挪移，并做傅立叶变换对时间(t)的积分。在一开始的连续的短时聚傅立叶变换(STFT)中，所表现的是从负无限大到正无限大，写成数学形式为：

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

$\mathbf {STFT} \{x()\}\equiv X(t,w)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-jw\tau }\,d\tau \qquad (1)$

其中$w(t)$是窗函数，而$x(t)$是代转换的讯号。$X(\tau ,\omega )$本质为$x(t)w(t-\tau )$的傅立叶变换，可以表示为时间和频率上的相位与强度。 但在实做上面，不太可能拥有无限长的时间来做分析，因此选定了一个范围内的时间来做分析，简单来说就是把w(t)使用一个在t=[-B B]，强度为1的方波来替换，写成数学形式表示为：

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

离散之方波-短时距傅立叶变换(Rec-STFT)

因讯号的分析上面，要做压缩或是简化时候，通常会使用离散的时间来考虑，因此，离散的短时距傅立叶变换被使用来分析离散的讯号。

第一种表示方式，称做为直观的改写，写成数学形式表示为：

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

对应(2)， $t=n\Delta _{t}\ ,f=m\Delta _{f}\ ,\tau =p\Delta _{t}\ ,B=Q\Delta _{t}\ ,Q=\left\lfloor {\frac {B}{\Delta _{t}}}\right\rfloor$

并满足取样定理$\Delta _{t}<{\frac {1}{2max(f)}}$，如此一来就能得到离散的短时间傅立叶变换，并由此来分析讯号，但若使用直观的改写来执行，虽然可以得到数据，但复杂度相对的复杂，所花费的时间过多，因此，便有简化的方式产生，由上式可知，是类似离散傅立叶变换，因此我们可以由离散的数学式子使用快速傅立叶变换(FFT)来改写。

第二种表示方式为，以FFT为基础的改写方式，写成数学形式表示为：

$X(n\Delta _{t},m\Delta _{f})=\Delta _{t}e^{j{\frac {2\pi (Q-n)m}{N}}}\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}}\qquad (4)$

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

由上式可知，在右半边的$\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}}$可以利用快速傅立叶变换(FFT)为N来得出结果，因此复杂度会比第一种直接改写的方式来的快速许多，但相对来说，就有几个限制条件需要符合，第一：$\Delta _{t}\Delta _{f}={\frac {1}{N}}$，N为整数，第二：$N\geq 2Q+1$，第三：取样定理，若满足以上条件，就能使用第二种型式的运算，接着因为$w(t)$使用方波讯号，在离散傅立叶变换(DFT)中会与下一个值的数值有关，因此可以更加简化DFT的运算量。

所以第三种形式为，基于快速傅立叶变换的递回公式，写成数学形式表示为：

$X(n\Delta _{t},m\Delta _{f})=X((n-1)\Delta _{t},m\Delta _{f})\ -\ x((n-Q-1)\Delta _{t})e^{-j{\frac {2\pi (n-Q-1)m}{N}}}\Delta _{t}\ +\ x((n+Q)\Delta _{t})e^{-j{\frac {2\pi (n+Q)m}{N}}}\Delta _{t}\qquad (5)$

但只有在$w(t)$为方波的情况时才能满足，但从原本的连续STFT来看时，所用的$w(t)$限制并不只在于方波，因此第三种方在实做不是很常见，接着因为可以从DFT来改写数学形式，因此使用CZT(chirp-Z转换)来做为使用的方法为：

$X(n\Delta _{t},m\Delta _{f})=e^{-j\pi m^{2}\Delta _{t}\Delta _{f}}\sum _{p=n-Q}^{n+Q}x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}e^{j\pi (p-m)^{2}\Delta _{t}\Delta _{f}}\Delta _{t}\qquad (6)$

若使用CZT来表示STFT时，虽然运算量平均会比第二，三种以FFT基础模型来的慢，但是可以放宽一个限制条件，$\Delta _{t}\Delta _{f}={\frac {1}{N}}$，N为整数，这也是一种可行的改写方式，由连续的STFT变成离散的STFT，再使用$w(t)$为方波的情况底下，可以有四种方式来计算出最后相同的结果，因此复杂度(运算量)的问题相对来说非常重要。

复杂度

有五种方式都可以得到出最后的结果，因此要来讨论复杂度(运算量)的高低，才能更快速的把讯号处理，假设n有T点，m有F点，首先来看第一种直接的改写，又称为暴力法，每个n点和m点都要做2Q+1次的运算，我们可以简单的看出其复杂度为：

$TF(2Q+1)\qquad (7)$

第二种以FFT为基础的模型，因为每个FFT所需乘法个数为 $(N/2)\log _{2}(N)$，加法个数为$N\log _{2}(N)$，其乘法运算较复杂为主要复杂度指标，因此每做一个FFT所需运算量为$(N/2)\log _{2}(N)$，并且有n个点，其复杂度为：

$T{\frac {N}{2}}\log _{2}(N)\qquad (8)$

第三种FFT为基础的递回公式，只要先算出一个n点的FFT，其他的点数用递回的方式就能算出，因此后面是剩下的T-1个点对应到每个F点上的运算量，其复杂度为：

${\frac {N}{2}}log_{2}(N)\ +\ 2(T-1)F\qquad (9)$

第四种是将STFT中移动窗型的时域信号所做的离算傅立叶变换(DFT)进行递回运算，先在时域信号的前面插入窗函数长度减一的0并且假设并且做DFT，因为一序列的0做DFT也会等于0，因此就可以省去初始进行FFT或是DFT的运算，随后在使用递回公式即可，由DFT的公式可知:

假设起始为一个N的DFT如下式:

$X_{0}[k]=\sum _{k=0}^{N-1}x[n]e^{-j2\pi nk/N}\qquad (10)$

则当窗函数移动一格时改为做底下的运算:

$X_{1}[k]=\sum _{k=0}^{N-1}x[n+1]e^{-j2\pi nk/N}\qquad (11)$

那么如果想要直接从第10式的运算结果推得第11式，比较两条式子可知:

$X_{1}[k]=(X_{0}[k]-x[0]+x[N])e^{j2\pi k/N}\qquad (12)$

由上式可知，只要将前面做完DFT的结果进行递回运算并且乘上DFT之回旋因子可得下一级之结果，因此递回的DFT算法其复杂度为:

$2TF\qquad (13)$

由于此法不是使用FFT作为用算基底，因此不受窗函数长度之限制，在运算上指示进行相同的动作，因此也式芯片设计上做容易被实作。

第三种级第四种方法所使用的递回公式式源自于DFT回旋因子(Twiddle factor)之循环性，因此所使用的窗函数也必须满足循环性，目前所采用的窗函数中只有矩形函数(Rectangular window)满足这个特性，若是采没有循环性之窗函数(高斯窗，三角窗，Hamming窗等)，则递回法皆不适用。

最后一种使用CZT来改写的方式，可以套用CZT的复杂度来计算在$x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}$部分有2Q+1点，在$e^{j\pi (p-m)^{2}\Delta _{t}\Delta _{f}}$部分有2Q+F点，所以要做4Q+F点的FFT，并且做两次，还有$x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}$部分附加的2Q+1的乘法和外面F点的乘法，对于每个n点都要做一次，因此其复杂度为：

$T({\frac {4Q+F}{2}}\log _{2}(4Q+F)*2\ +\ 2Q\ +\ 1\ +\ F)\ \approx \ T((4Q+F)\log _{2}(4Q+F))\qquad (14)$

短时距傅立叶变换的窗函数

短时距傅立叶转换(STFT)可以视为在一段时域信号上用一移动的的窗函数取出片段进行傅立叶转换的结果，可以获得该时间点的频谱分析，然而此法并无法得到该时间点的顺时频率，随着所采用的窗函数的不同，会在分析结果的时域及频域上有不同的分辨率及效果。

对称型窗函数

在作时频时频分析时，如果目标是一段完整著时域信号，则通常采用对称型的窗函数，一般而言对称型的窗函数必须拥有满足以下的条件:

1.$w(t)=w(-t)\qquad (15)$

2.$max(w(t))=w(0),w(t_{1})\geq w(t_{2})\quad \forall \quad |t_{2}|>|t_{1}|\qquad (16)$

3.$w(t)\thickapprox 0\quad \forall t>>1\qquad (17)$

以常用的矩形窗为例，窗函数的大小会决定时域及频域上有不同的影响，假如取较小的窗函数，更能接近时域上的真实性，因此时域上的分辨率更好，但是作傅立叶转换的长度降低了，就造成了频域上的分辨率降低。

反之，若是取较大的矩形窗，有足够的长度来作DFT，虽然频域上分辨率会变好，但时域上就会相应降低，如下图:

除了矩形窗以外，高斯窗也是常用的窗函数，能够同时在时频域上都有不错的表现。以高斯函数作为窗函数的短时距傅立叶转换，又称为"加伯转换(Gabor Transform)"

非对称窗型函数

在做时频分析时，有时也必须面对及时处理的情况，也就是没有完整的时域信号，只有截止目前的信号，例如地震波的分析或是语音的分析，则传统的对称型窗函数则不能达到需求，在实务上若是分析必须使用未来的信号称为"非因果的(non-causal)"，这样的系统没有办法被实作，解决的方法时将信号进行一段的延迟，再进行分析，这样一来便可实作，但是必须等同时窗函数右边的宽度的延迟作为代价，所以在作及时分析的时候，有时会将窗函数的一边(通常是未来信号的那一边)缩短来达到短延迟及时分析的效果，如下图:

参见

• 正交转换
• 傅立叶级数
• 连续傅立叶变换
• 离散时间傅立叶变换
• 离散傅立叶变换
• 傅立叶分析
• 拉普拉斯转换
• 小波转换

参考

• Jian-Jiun Ding, Time frequency analysis and wavelet transform class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2009.
• Jian-Jiun Ding, Time frequency analysis and wavelet transform class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2018.
• Z. Kollar, F. Plesznik and S. Trumpf, "Observer-Based Recursive Sliding Discrete Fourier Transform [Tips & Tricks]," in IEEE Signal Processing Magazine, vol. 35, no. 6, pp. 100-106, Nov. 2018.doi: 10.1109/MSP.2018.2853196
• http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=8498081&isnumber=8496891