短时距傅立叶变换
短时距傅立叶变换(Short-time Fourier Transform, STFT)是傅立叶变换的一种变形,也称作加窗傅里叶变换(Windowed Fourier transform)或Time-dependent Fourier transform,用于决定随时间变化的信号局部部分的正弦频率和相位。实际上,计算短时距傅立叶变换的过程是将长时间信号分成数个较短的等长信号,然后再分别计算每个较短段的傅立叶转换。通常拿来描绘频域与时域上的变化,为时频分析中其中一个重要的工具。
与傅立叶转换在概念上的区别.
将讯号做傅立叶变换后得到的结果,并不能给予关于信号频率随时间改变的任何资讯。以下的例子作为说明:
formula_1
傅立叶变换后的频谱和短时距傅立叶转换后的结果如下:
由上图可发现,傅立叶转换只提供了有哪些频率成份的资讯,却没有提供时间资讯;而短时傅立叶转换则清楚的提供这两种资讯。这种时频分析的方法有利于频率会随著时间改变的信号,如音乐信号和语音信号等分析。
定义.
连续短时傅立叶转换.
简单来说,在连续时间的例子,一个函数可以先乘上仅在一段时间不为零的窗函数再进行一维的傅立叶变换。再将这个窗函数沿著时间轴挪移,所得到一系列的傅立叶变换结果排开则成为二维表象。数学上,这样的操作可写为:
formula_2
另外也可用角频率来表示:
formula_3
其中formula_4是窗函数,窗函数种类有很多种,会在稍后再做仔细讨论。formula_5是待变换的讯号。formula_6是formula_7的傅立叶变换。 随著formula_8的改变,窗函数在时间轴上会有位移。经formula_7后,信号只留下了窗函数截取的部分做最后的傅立叶转换,所得到的结果为一复数函数,代表著信号随时间与频率变化的大小与相位。
离散短时傅立叶转换.
在离散时间的例子,资料会被切割成数个大量的帧,而每组帧通常会互相重叠,避免因切割方式造成边界的误差。而每组帧在各自进行傅立叶转换后所得的复数结果会再进行相加,可得到每个点时间与频率变化的大小与相位。数学上,这样的操作可写为:
formula_10
相同地,其中formula_11是窗函数,formula_12是待变换的讯号。在这个例子里,m是离散的且ω是连续的,但大部分实际的应用当中,短时距傅立叶转换在电脑中都是以快速傅立叶转换进行计算(见实现方法的快速傅立叶变换),而此时这两个参数都是离散且被量化的。
Sliding 离散傅立叶转换.
当只想要得知特定少数的ω,或是短时距傅立叶转换每次窗函数移动m的值,则短时距傅立叶转换可以利用sliding DFT演算法更有效地计算出来。
反短时距傅立叶转换.
短时距傅立叶转换是可逆的,也就是说原本的信号可以借由反短时距傅立叶转换将短时距傅立叶转换后的信号还原。
其中最广为接受的反短时距傅立叶转换方法是重叠-相加之折积法,此方法也促成了更多样的信号处理方法。
反短时距傅立叶转换,其数学类似傅立叶转换,但须消除窗函数的作用,首先必须先将窗函数的总面积规模化使得
formula_13
而从上也可轻易地得出
formula_14
和
formula_15
连续傅立叶转换公式如下:
formula_16
将formula_5进行上述的替换:
formula_18
formula_19
将积分顺序进行交换:
formula_20
formula_21
formula_22
因此傅立叶转换可以视为某种将formula_5所有的短时距傅立叶转换的相位同调部分进行相加。
而反傅立叶转换公式如下:
formula_24
因此 formula_5可以从formula_26被复原
formula_27
或
formula_28
与上面所列的窗函数的式子进行比较,可得
formula_29
对反傅立叶转换公式中的formula_26来说formula_31是不变的
formula_32
另外用角频率来表示:
formula_33
窗函数.
窗函数通常满足下列特性:
常见的窗函数有:方形、三角形、高斯函数等,而短时距傅立叶转换也因窗函数的不同而有不同的名称。而加伯转换,即为窗函数是高斯函数的短时距傅立叶转换,通常没有特别说明的短时距傅立叶转换,即为加伯转换。
非对称窗函数.
当在特殊应用时,窗函数特性的第一点可以不满足,如下图的非对称窗函数formula_38,其中formula_39。左图为窗函数原本的图形,而在计算短时距傅立叶变换时,需将窗函数转到formula_31轴上得出formula_41,换言之,欲得到的短时距傅立叶变换的结果需在formula_42的时间点才能算出,因此若formula_43愈小,即可愈快得结果,此种非对称窗函数可应用在地震波、碰撞侦测...等,需要即时处理的应用。
方形窗函数的短时距傅立叶转换.
概念.
右图即为方形窗函数的一个例子,其数学定义:
formula_44
可以随要分析的信号,来调整B的大小(即调整方形窗函数的宽度)。至于B的选择,将会在下面探讨。
短时傅立叶转换可以简化为
formula_45
反短时傅立叶转换可简化为
formula_46
特性.
其大部分的特性都与傅立叶转换的特性相对应
formula_47
formula_48
formula_49
若有一信号formula_50,formula_51分别为formula_52做方形窗函数短时 距傅立叶转换的结果,则formula_53。
formula_54
formula_55
1. 当formula_56,
formula_57
2. 当formula_58,
formula_59
formula_75
方形窗函数宽度formula_60的选取.
结果如右图所示,B越大则在频率变化处(t = 10, 20)附近的频率越不准确,即可能会有多个频率成分出现。但同时,其他时间点的能量则较集中;没有如B较小时,频率散开或模糊的情形。
上述也是其中一个小波转换及多解析度分析作为改进的方向,其中多解析度分析能在高频时有较好的时间轴解析,而在低频时能有较好的频率轴解析,此种组合较契合许多实际的应用。
时间轴与频率轴的解析度无法同时提升也与海森堡不确定性原理有关,即时间与频率的标准差乘积有所限制,而高斯函数恰好能符合不确定性原理的极值,也就是两者同时达到最好的解析度,而应用高斯函数的时频分析方法即为加伯转换,而在经过修改及多解析度分析后,成为了莫莱小波。
其他窗函数.
高斯窗函数.
概念.
高斯窗函数的短时距傅立叶转换又称为加伯转换。以下是高斯函数的数学定义,
formula_76
据此,短时傅立叶转换可以写为
formula_77
三角形窗函数.
概念.
三角形窗函数如右图所示,数学定义如下,
formula_78
formula_79
可使用在震幅改变的情况下,相对于方形窗函数,可更好的滤除杂讯。
海宁(Hanning/ Hann)窗函数.
概念.
海宁函数如右图所示,数学定义如下,
formula_80
相较于三角形窗函数,海宁窗函数更为贴近现实讯号的趋势,可进一步滤除杂讯。
汉明(Hamming)窗函数.
概念.
汉明窗函如右图所示,数学定义如下,
formula_81
跟海宁窗函数类似,但两端不为零。
海宁与汉明窗的区别.
窗函数有四个指标,分别为
因为汉明窗两端不能到零,而海宁窗两端为零。从以上频率响应来看,汉明窗可以有效减少靠近的旁办,但在较远的旁办泄漏比海宁窗严重。
如何决定窗函数.
可根据以下条件来选取窗函数,
在决定复杂度跟解析率后,可利用不同的窗函数达到更好的滤杂讯效果。
瑞利频率.
当Nyquist频率是能被有意义分析的频率最大值的限制,而瑞利频率则是能被有限频宽频的窗函数解析的频率最小值的限制。若给定一窗函数的长度是T秒,最低能被解析的频率即为1/T Hz。
瑞利频率在短时距傅立叶变化的应用中扮演重要的角色,像是在分析神经信号时。
频谱(Spectrogram).
Spectrogram即短时傅立叶转换后结果的绝对值平方,两者本质上是相同的,在文献上也常出现spectrogram这个名词。
formula_82
应用.
短时距傅立叶变换及其他工具经常用于分析音乐。
如右图所示,
音频工程师使用这种视觉来获取有关音频样本的信息。
此外,因频率会随时间而改变,短时距也可使用在以下情境,
若与时间无关,如卷积,照片等则不能使用短时距傅立叶变换来进行分析。而影片属于3D讯号,其短时距傅立叶产物为6D讯号,故也不适用。
短时距傅立叶变换实现方法.
从连续短时距傅立叶变化的定义出发
formula_83
令 formula_84 ,则上述式子时域可从连续转为离散
formula_85
若当formula_86
上式可改写为
formula_87
直接运算.
限制条件.
(1)要满足Nyquist criterion
formula_88
formula_89的频宽为formula_90。而formula_91的频宽则为formula_92,formula_93的频宽也为formula_92
因为在时域相乘相当于在频域做折积,因此formula_95的频宽为formula_96(通常formula_97会远大于formula_98,所以主要影响频宽的是formula_97)
formula_100
转换到离散形式(formula_101),其中formula_102
formula_103,由于无限大的上下限实务上做不到,所以尝试变成有限大的上下限。
假设formula_104 for formula_105
formula_106
formula_108
假设t-axis有T个取样点,f-axis有F个取样点,则我们总共要对TF个点做formula_109次的运算,因此可得复杂度为formula_110
优点:简单及有弹性(因为限制少)
缺点:复杂度较高
快速傅立叶变换.
限制条件.
(1)要满足Nyquist criterion
formula_111
推导.
标准的离散傅立叶转换式子为
formula_114
由直接运算得知如下公式
formula_115
因此为了让上式符合离散傅立叶转换的上下界,令formula_116代入上式即可得
formula_117
其中
formula_118
运算步骤.
假设formula_119
formula_120
步骤一:计算formula_121
步骤二:formula_122
步骤三:决定formula_123
步骤四:formula_124
步骤五:转换formula_125成formula_126
步骤六:设formula_127,并回到步骤三,直到formula_128
formula_129
借由取样定理可得知formula_130
假设formula_131及formula_132,则经由formula_133可得formula_134
formula_135及formula_136,则经由formula_137可得formula_138
步骤一:formula_139
步骤二:formula_140
步骤三:计算formula_141
步骤四:利用求得的formula_123计算快速傅立叶转换
formula_143
步骤五:转换formula_125到formula_126
formula_146
formula_148
因此可将上式改为formula_149,其中formula_150代表取m除以N的余数
步骤六:设定formula_127,回到步骤三直到formula_152
时间复杂度.
利用FFT计算formula_153,其中每次FFT的时间复杂度为
formula_154
总时间复杂度为formula_155
优缺点.
优点:与直接运算相比,复杂度较低
缺点:较多限制,包括
formula_156
使用快速傅立叶变换加上递回关系式.
限制条件.
(1)要满足Nyquist criterion
formula_157
(2)formula_112
(3)formula_113
(4)需为方形窗函数的短时距傅立叶转换
推导.
因为是方形窗函数
formula_160,因此原式可由此关系变成以下式子
formula_161
而由此可看出n和n-1有递回关系,如下
formula_162
(1)以FFT计算formula_163
其中formula_164
(2)利用递回关系式计算算formula_165
则formula_166
时间复杂度.
(1)FFT计算一次
formula_163
(2)利用递回关系,计算formula_169时的数值,因此共会执行T-1次递回,如下式
formula_166
每次递回都要计算formula_171及formula_172两个乘法(相当于2F的复杂度)
总时间复杂度 formula_174
优缺点.
优点:四种运算中,最低的复杂度formula_175
缺点:
使用Chirp-Z 转换.
限制条件.
(1)要满足Nyquist criterion
formula_157
推导.
令formula_177
即可由直接运算的式子导出Chirp_Z变换的式子,如下所示
formula_178
运算步骤.
Step1:formula_179 formula_180
Step2:formula_181
Step3:formula_182
时间复杂度.
当n为定值时
(1)假设formula_183 相乘时间复杂度为2Q+1
(2)令formula_184,则formula_185 convolution时间复杂度为 formula_186
(3)formula_187相乘时间复杂度为 F
因此,总时间复杂度为formula_188
虽然此实现方法和使用FFT计算的时间复杂度相同,但因为convolution相当于做三次FFT,因此实际操作时运算时间约为使用FFT计算的2~3倍
优缺点.
优点:只有一项限制:formula_189
缺点:与前四种相比,复杂度是中间的。
Unbalanced Sampling for STFT and WDF.
将直接法和快速傅立叶转换方法做修正
1.直接法.
formula_100
修正后 :formula_191
其中, formula_192 ,formula_193
假设formula_194 for formula_195,则上下限可借由以下推导而修正
formula_196
则上限可以写成formula_197,下限则以此类推
注:formula_198(输入讯号的取样间隔)
formula_199(在t轴上的输出讯号的取样间隔)
然而,formula_200是整数会是比较好的。
formula_201
则经由上述公式可求得S=441,代表经由unbalanced sampling,我们跟原本formula_202相比可减少441倍的取样点。
时间复杂度.
由于t轴的取样点少了S倍,因此跟原本的直接运算复杂度相比,只要把formula_203即可,如下:
复杂度:formula_204
2.快速傅立叶转换.
限制条件.
(1) formula_205
(3) formula_208,formula_209的频宽是 formula_210
i.e. formula_211 ,当 formula_212
过程.
formula_213
令 formula_214
formula_215 for formula_216
formula_217for formula_218
修正后:formula_219
运算步骤.
假设formula_220
formula_221
formula_222
步骤一:计算formula_223
步骤二:formula_224
步骤三:决定formula_123
步骤四:formula_124
步骤五:转换formula_227
步骤六:设定formula_127及返回步骤三,直到formula_229
复杂度.
formula_204
Non-Uniform formula_199.
(1) 先用比较大的formula_199
(2) 如果发现formula_233 和 formula_234 之间有很大的差异,则在formula_235,formula_236 之间选用比较小的取样区间formula_237
再用Unbalanced Sampling for STFT and WDF 中修正后的快速傅立叶转换方法算出 formula_241,formula_242formula_243
(3) 以此类推,如果 formula_244的差距还是太大,则再选用更小的取样间隔formula_245
若有一音乐信号总共有1.6秒,formula_249
formula_254,共29点
可以这样做的原因为:有些音乐讯号在和弦与和弦中间几乎没有变化,因此可以挑选较大的formula_199取样;和弦在变换时,频率会变化的较剧烈,因此变换和弦是需要用较多的取样点。借由此种non-uniform的取样,可以让我们大幅减少运算量,从最一开始的formula_256可看出我们的运算量大幅降低。
生成维基百科快照图片,大概需要3-30秒!