前言

在语音识别和话者识别中,需要将声音编码为数字信息,提取其中的特征进行处理,而最常用的语音特征便是梅尔倒谱系数MFCC。为了将频谱的包络与细节分离开来,以便分析声道共振特征和基音周期,需要把这种非线性问题转化为线性问题。MFCC预处理之后的第一步便是通过快速傅里叶变换(Fast Fourier Transform)将时域信号转变为信号频谱。而快速傅里叶变换是离散傅里叶变换(Discrete Fourier Transform)的改进算法,将其时间复杂度从 O ( n 2 ) O(n^{2}) O(n2)改进为 O ( n log ⁡ n ) O(n\log{n}) O(nlogn)

离散傅里叶变换 DFT

傅里叶变换能将满足一定条件的函数表示成正弦和/或余弦函数或它们的积分的线性组合。根据函数在时域的特性分别有如下几种变换:

变换时域频域
连续傅里叶变换连续,非周期性连续,非周期性
傅里叶级数连续,周期性离散,非周期性
离散时间傅里叶变换离散,非周期性连续,周期性
离散傅里叶变换离散,周期性离散,周期性

这四种傅立叶变换都是针对长度为无穷大的信号,但这对于计算机处理来说是不可能的。那么有没有针对长度有限的傅立叶变换呢?很遗憾没有,因为正余弦波被定义成从负无穷大到正无穷大,我们无法把一个长度无限的信号组合成长度有限的信号。
面对这种困难,解决方法是把长度有限的信号表示成长度无限的信号。例如,可以把信号无限地从左右进行延伸,延伸的部分用零来表示,这样这个信号就可以被看成是非周期性离散信号,我们可以用到离散时域傅立叶变换(DTFT)的方法。也可以把信号用复制的方法进行延伸,这样信号就变成了周期性离散信号,这时我们就可以用离散傅立叶变换方法(DFT)进行变换。
但是对于非周期性的信号,我们需要用无穷多不同频率的正弦曲线来表示,这对于计算机来说是不可能实现的。所以对于离散信号的变换只有离散傅立叶变换(DFT)才能被适用,对于计算机来说只有离散的和有限长度的数据才能被处理,其它的变换类型只有在数学演算中才能用到。

对于N点序列 { x [ n ] } 0 ≤ n ≤ N \left \{ x[n] \right \}_{0≤n≤N} {x[n]}0nN,它的离散傅里叶变换为:
图片描述
其中N表示傅里叶变换的点数,k表示傅里叶变换的第k个频谱,W称为旋转因子。通过欧拉公式进一步变换我们可以得到:
图片描述
由于余弦信号也是一种相位移动的正弦信号,因此变换后的第k点的实部和虚部对应的是一个相同频率的正弦信号,只是相位不同。随着k的变换正弦信号的频率也在发生变换,因此DFT将时序信号转化成不同频率的正弦信号的叠加。
DFT也可以用矩阵表示为 X = W N x X=W_Nx X=WNx,其中:
图片描述

快速傅里叶变换 FFT

快速傅里叶变换(Fast Fourier Transform),是利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法,于1965年由J.W.库利和T.W.图基提出。
对于多项式 f ( x ) = ∑ i = 0 n a i x i f(x)=\sum_{i=0}^{n}{a_i x^i} f(x)=i=0naixi f ( x ) = ∑ i = 0 n b i x i f(x)=\sum_{i=0}^{n}{b_i x^i} f(x)=i=0nbixi,定义其乘积 f g fg fg如下:
图片描述
显然,上式可以 O ( n 2 ) O(n^{2}) O(n2)的复杂度来计算乘积每一项的系数,但FFT可以 O ( n log ⁡ n ) O(n\log{n}) O(nlogn)的时间复杂度来计算这个乘积。

对于任意n次多项式 f ( x ) = ∑ i = 0 n a i x i f(x)=\sum_{i=0}^{n}{a_i x^i} f(x)=i=0naixi,除去这种最常见的表示方法,我们也可以根据其系数的有序排列 ( a 0 , a 1 , a 2 , ⋯ , a n ) (a_0,a_1,a_2,⋯,a_n) (a0,a1,a2,,an)来唯一确定 f ( x ) f(x) f(x),这就是多项式的系数表示。同时,由代数基本定理,n次多项式上n+1个不同的点能唯一确定一个多项式,因此 f ( x ) f(x) f(x)还可以被写作 { ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , ⋯ , ( x n , f ( n ) ) } \left \{ (x_0,f(x_0)),(x_1,f(x_1)),⋯,(x_n,f(_n)) \right \} {(x0,f(x0)),(x1,f(x1)),,(xn,f(n))},这就是多项式的点值表示。在 x 0 , x 1 , ⋯ , x n x_0,x_1,⋯,x_n x0,x1,,xn确定时,将f的点值表示记为 τ ( f ) τ(f) τ(f)
如果两个多项式 f ( x ) f(x) f(x) g ( x ) g(x) g(x)的点值表示法分别为:
图片描述
则二者乘积的点值表示法为:
图片描述
由上式可知,多项式的乘法在点值表示法下可以 O ( n ) O(n) O(n)的复杂度计算,所以如果能够在较低的时间复杂度内将系数表示转化为点值表示,再将点值表示转回系数表示,就能以较低的时间复杂度计算多项式的乘法。这也正是快速傅里叶变换的基本思路。

对于多项式 f ( x ) = ∑ i = 0 n − 1 a i x i ∈ X n − 1 f(x)=\sum_{i=0}^{n-1}{a_i x^i}∈X_{n-1} f(x)=i=0n1aixiXn1,不失一般性地,设 n = 2 s n=2^s n=2s s ∈ N s∈\mathbb{N} sN(可以将一个多项式等价地看作是次数更高的高次项系数均为零的多项式,故可以将n看作第一个等于或大于它的2的整数次幂)我们将n次单位根 ω n = cos ⁡ 2 π n + i sin ⁡ 2 π n ω_n=\cos \frac{2π}{n}+i\sin \frac{2π}{n} ωn=cosn2π+isinn2π 0 0 0 ~ n − 1 n-1 n1次幂分别代入 f ( x ) f(x) f(x),便得到其点值向量 { ( ω n 0 , f ( ω n 0 ) ) , ( ω n 1 , f ( ω n 1 ) ) , ⋯ , ( ω n n − 1 , f ( ω n n − 1 ) ) } \left \{ (ω_n^0,f(ω_n^0)),(ω_n^1,f(ω_n^1)),⋯,(ω_n^n-1,f(ω_n^{n-1})) \right \} {(ωn0,f(ωn0)),(ωn1,f(ωn1)),,(ωnn1,f(ωnn1))}
可以看到如果进行朴素代入,算法的时间复杂度依然是 O ( n 2 ) O(n^{2}) O(n2)。但利用单位根如下两条性质,我们可以将算法进一步简化。
图片描述
按下标的奇偶性将 f ( x ) f(x) f(x)分为两组,并令:
图片描述
则有 f ( x ) = f 1 ( x 2 ) + x f 2 ( x 2 ) f(x)=f_1(x^2)+xf_2(x^2) f(x)=f1(x2)+xf2(x2)。设 0 < k < n 2 0<k<\frac{n}{2} 0<k<2n k ∈ Z k∈\mathbb{Z} kZ,分别将 x = ω n k x=ω_n^k x=ωnk x = ω n k + n 2 x=ω_n^{k+\frac{n}{2}} x=ωnk+2n代入 f ( x ) f(x) f(x)可得:
图片描述
k k k k + n 2 k+\frac{n}{2} k+2n取遍了 [ 0 , n − 1 ] [0,n-1] [0,n1]上所有的n个整数。观察上述两式的结构,我们只需要求出 f 1 ( ω n 2 k ) f_1(ω_{\frac{n}{2}}^k) f1(ω2nk) f 2 ( ω n 2 k ) f_2(ω_{\frac{n}{2}}^k) f2(ω2nk)即可在 O ( n ) O(n) O(n)的时间复杂度内求出 f ( x ) f(x) f(x)的值,两者的规模都是原问题一半。因此可以采取递归求解,将问题转化为求出 f 1 ( ω n 4 k ) f_1(ω_{\frac{n}{4}}^k) f1(ω4nk) f 2 ( ω n 4 k ) f_2(ω_{\frac{n}{4}}^k) f2(ω4nk),直至最后转化为求 f 1 ( ω 1 k ) = f 2 ( ω 1 k ) = 1 f_1(ω_1^k)=f_2(ω_1^k)=1 f1(ω1k)=f2(ω1k)=1。时间复杂度为 T ( n ) = T ( n 2 ) + O ( n ) = O ( n log ⁡ ⁡ n ) T(n) = T(\frac{n}{2})+O(n) = O(n\log⁡n) T(n)=T(2n)+O(n)=O(nlogn)

Reference:
语音识别特征提取方法
从头到尾彻底理解傅里叶变换算法
离散傅里叶变换 - 维基百科
傅里叶变换 ~ 离散傅里叶变换(DFT)
一小时学会快速傅里叶变换(Fast Fourier Transform)
快速傅里叶变换(FFT)超详解

更多推荐