快速傅立叶变换（Fast Fourier Transform）
FFT[算法]是由Cooley和Tukey在其论文《机器计算复傅里叶级数的方法》中首次提出，并基于复数DFT算法构建而成。
通过ComplexDFT来计算Real DFT
虽然Fast Fourier Transform (FFT)算法建立在Compute Complex Discrete Fourier Transform (CCDFT)的基础上,但其依然可以通过这一方法计算出Real Discrete Fourier Transform (RDCT).这是因为RDCT可以通过相对简便的方式转换为CCDT.通过查看图2-1中的信息图（Information Graph）,我们可以清晰地比较出两种变换的主要区别.在Real DFT中,信号的时间域由N个点构成,而频率域则分为实部信号与虚部信号两个部分,其长度均为N/2+1.相比之下,在CCDT中,时间域同样包含实部与虚部两个部分,但它们的长度分别为N.而在频率域方面,CCTF则将实部与虚部的长度扩展至N.
如图2-1所示，在时间域上ComplexDFT相对于RealDFT主要增加了引入了一个虚数分量这一关键操作；而频率域长度的变化则直接导致了这一技术差异的存在。
图2-1 Real DFT和ComplexDFT的区别
该区别的根本原因在于实数其其虚数部分为零；由此可知将实数值转换成复数值非常简单只需要附加一个系数设为零的虚数值分量即可；例如，在图2-1所示的ComplexDFT中若将时间域中的虚数值置零那么在频域中新增的部分也应置零这样一来图2-1中的RealDFT与ComplexDFT就会完全相等；当包含负频率时在DFT所对应的频域呈现周期性特性；具体而言在ComplexDFT中正频率范围是从0到N/2而负频率则分布在N/2+1到N-1区间内
相较于使用ComplexDFT来计算Real DFT的方法而言, 使用ComplexInverse DFT计算Real InverseDFT更为复杂. 其运算过程相对直接, 其中在频域中位于N/2+1至N-1范围内的各个频率分量的系数值与前半部分相对应的位置取反数值. 即运算过程中各个对应的频率分量值均为前半部分相应位置取反数值.
系数（N/2+1）=—系数（N/2-1）
系数（N/2+2）=—系数（N/2-2）
需要注意的是，在进行RealInverse DFT运算时，在0和N/2位置上没有对应的点存在。随后采用一个子程序来处理从N/2+1至N-1范围内的系数。具体而言，在进行RealInverse DFT运算时，在进行复数域DFT（complexDFT）之前需要执行以下步骤：首先将从0到N/2的所有系数复制至complexDFT中对应的区域，并利用预设算法对后续区域进行处理以完成整个变换过程。
6000'NEGATIVE FREQUENCY GENERATION
The function generates the complex frequency domain from the real one.
6020'Upon invocation of this subroutine, N% represents the count of data packets within this signal stream. Additionally, it includes information regarding those data packets.
6030'REX[ ] and IMX[ ] contain the real frequency domain in samples 0 toN/2.
Upon returning, REX[] and IMX[] hold the complex frequency domain in samples from 0 through N−1.
6050 '
6060FOR K% = (N%/2+1) TO (N%-1)
6070REX[K%] = REX[N%-K%]
6080IMX[K%] = -IMX[N%-K%]
6090NEXT K%
6100 '
6110RETURN
FFT的实现原理
FFT算法确实较为复杂，在本研究中不对细节展开讨论，并仅阐述其实现原理。在复数域中（即所谓的虚数域），时间域与频率域中的信号均是由N个复数值构成的序列数据集。每个复数值都包含一个实部与一个虚部两个基本数值参数来表示其特征信息。例如，在位置k=6处的复数值X[6]（即通常表示为X_6），其具体表示为实部Re(X_6)与虚部Im(X_6)两个部分共同作用的结果。
FFT算法的基本原理是将时间域中的一个由N个点组成的信号划分为N个仅含有一个点的子信号。接着分别计算这N个子信号在频域中的响应值。最后将这些频域响应综合生成完整的频域信号。
图2-2描述了一个包含12个点的示例信号在FFT中的分解过程。
图2-2 FFT中的分解(decompose)过程
图2-2中的过程看似复杂，在实际应用中可以通过如图2-3所示的位反转算法（bitreversal sorting）来具体完成。该算法通过对各点二进制位进行对称反转的操作，即可实现将N个点的信号分解为N个单点信号的过程。
图2-3位反转排序
该算法的下一步骤是计算这些N个单点信号在频域中的振幅值。这被认为是该算法中相对简单的步骤之一。每个单独信号自身的振幅量与其数值相等，在这一阶段无需进行任何计算或处理。
算法的最后一项步骤是将这N个频率域的数据点按照时间域分解时的逆序顺序进行整合(combine)，其中不宜采用位翻转法作为组合方式；这一过程在算法中占据核心地位并最为复杂
注
图2-4展示了如何将两个长度为4的频率域子波合并成一个长度为8的目标子波的过程。这一合成操作必须与随后的时间域分解操作严格互逆关系。例如，在时间域中我们有两个子波abcd和efgh，则要将其整合为一个包含8个点的目标子波则需依次执行以下两步：首先对这两个子波进行延展处理（即将其延长至8个点长），即将每个位置插入一个零值；然后将这两个延展后的序列相加即可得到新的综合序列X=x₁x₂...x₈.
图2-4FFT组合（synthesis）
当时间域用0稀释时，对应的频率域会复制自己。
在时间域中进行移位后再以零填充时，在其对应的频率域将乘以一个三角函数，并进而进行复制。
与abcd不同的是efgh的稀释方式，在对abcd进行稀释时得到的结果是a₀b₀c₀d₀，并且其中偶数位置上的数值均为零；而当对efgh执行稀释操作时会得到的结果是₀e₀f₀g₀h，并且其中奇数位置上的数值均为零。由此可见，在时间域中将一个信号向右移动一位相当于在频率域中对其进行某种变换
图2-5呈现了两个频率域中长度为4的信号组合的过程。左侧Odd-Four Point Frequency Spectrum是指时间域信号中奇数位置零后的频率域表示（如EFGH），而右侧Even-Four Point Frequency Spectrum则是指时间域信号中偶数位置零后的频率域表示（如ABCD）。
为了更清晰地呈现这一过程, 通过图2-6展示其中关键的两点. 其因形状酷似一只展翅高飞的蝴蝶, 因此这一图形被科学界将其命名为butterfly process.
图2-5 FFT组合过程
图2-6 Butterfly
如图2-7所示的是FFT变换的整体流程图,它主要由三个关键模块构成:1.通过位操作算法可以实现时域信号的分解过程;2.将时域信号分解后得到的一组独立样本直接转换到频域时,由于每个样本在频域中的幅度值与时域幅度值相等,因此无需任何计算量即可完成转换;3.本算法的关键内容则集中在第三模块,这也是整个流程的核心内容
图2-7FFT流程图
在图2-7中所述之图示中包含三种不同的循环结构：该外层循环对应于从图2-7所示lgN层级范围内的合成操作；中间一层则负责每一层级的具体结合步骤；而最内层的部分则执行蝴蝶操作（butterfly process）。
下面是FFT算法的一段Basic代码：
1000'THE FAST FOURIER TRANSFORM
1010'Upon entry, N% contains the number of points in the DFT, REX[ ] and
1020'IMX[ ] contain the real and imaginary parts of the input. Uponreturn,
1030'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 toN%-1.
1040'
1050PI = 3.14159265 'Set constants
1060NM1% = N%-1
1070ND2% = N%/2
1080M% = CINT(LOG(N%)/LOG(2))
1090J% = ND2%
1100'
1110FOR I% = 1 TO N%-2 'Bit reversal sorting
1120IF I% >= J% THEN GOTO 1190
1130TR = REX[J%]
1140TI = IMX[J%]
1150REX[J%] = REX[I%]
1160IMX[J%] = IMX[I%]
1170REX[I%] = TR
1180IMX[I%] = TI
1190K% = ND2%
1200IF K% > J% THEN GOTO 1240
1210J% = J%-K%
1220K% = K%/2
1230GOTO 1200
1240J% = J%+K%
1250NEXT I%
1260'
1270FOR L% = 1 TO M% 'Loop for each stage
1280LE% = CINT(2^L%)
1290LE2% = LE%/2
1300UR = 1
1310UI = 0
1320SR = COS(PI/LE2%) 'Calculate sine & cosine values
1330SI = -SIN(PI/LE2%)
1340FOR J% = 1 TO LE2% 'Loop for each sub DFT
1350JM1% = J%-1
1360FOR I% = JM1% TO NM1% STEP LE% 'Loop for each butterfly
1370IP% = I%+LE2%
1380TR = REX[IP%]*UR – IMX[IP%]*UI 'Butterfly calculation
1390TI = REX[IP%]*UI + IMX[IP%]*UR
1400REX[IP%] = REX[I%]-TR
1410IMX[IP%] = IMX[I%]-TI
1420REX[I%] = REX[I%]+TR
1430IMX[I%] = IMX[I%]+TI
1440NEXT I%
1450TR = UR
1460UR = TRSR - UISI
1470UI = TRSI + UISR
1480NEXT J%
1490NEXT L%
1500'
1510RETURN
更快的FFT算法
存在多种方法以加快FFT算法的速度进行优化与实现，并能达到最高可达约20%至40%的加速比的效果。具体而言，在时间域进行分解时，在每个信号仅包含四个点时停止分解流程以减少计算量；另一种优化策略则是通过将时间域虚数部分置零来实现对复数FFT算法的有效转换，并将其简化为实数形式以降低计算复杂度。下面给出实数逆FFT算法的具体伪代码实现：
4000'INVERSE FFT FOR REAL SIGNALS
4010'Upon entry, N% contains the number of points in the IDFT, REX[ ]and
4020'IMX[ ] contain the real and imaginary parts of the frequency domainrunning from
4030'index 0 to N%/2. The remaining samples in REX[ ] and IMX[ ] areignored.
4040'Upon return, REX[ ] contains the real time domain, IMX[ ] containszeros.
4050'
4060'
4070FOR K% = (N%/2+1) TO (N%-1)'Make frequency domain symmetrical
4080REX[K%] = REX[N%-K%]'(as in Table 12-1)
4090IMX[K%] = -IMX[N%-K%]
4100NEXT K%
4110'
4120FOR K% = 0 TO N%-1'Add real and imaginary parts together
4130REX[K%] = REX[K%]+IMX[K%]
4140NEXT K%
4150'
4160GOSUB 3000‘Calculateforward real DFT (TABLE 12-6)
4170'
4180FOR I% = 0 TO N%-1'Add real and imaginary parts together
4190REX[I%] = (REX[I%]+IMX[I%])/N%'and divide the time domain by N%
4200IMX[I%] = 0
4210NEXT I%
4220'
4230RETURN
图2-8描绘了FFT中所采用的对称性原理。a与b分别代表同一个时域信号，并且其虚数部分全部为零。c对应于频率域中的实部而d则对应于频率域中的虚部。其中c表现出偶对称特性而d则表现出奇对称特性。
图2-8DFT中实数部分的对称
图2-9DFT中虚数部分的对称
两者结构相似，在时间域实部a等于零的情况下，在虚部b不为零时，在频率域中对应的曲线c呈现奇对称特性、曲线d呈现偶对称特性。
两者结构相似，在时间域实部a等于零的情况下，在虚部b不为零时，在频率域中对应的曲线c呈现奇对称特性、曲线d呈现偶对称特性。
上文阐述了时间域中某一部分为零的情况。当时间域的实部与虚部均不为零时会出现什么情况？频率域可通过将两个或多个频谱相加得到。关键在于频率域具有奇偶对称性质的波谱可被精确划分为两部分。输入信号可被分解成两部分：其中N/2个奇数位信号放置于时间域信号之实部；另外N/2个偶数位信号则放置于虚部。这一操作使得长度为N的FFT变换得以简化至长度N/2的操作。经过上述步骤后会得到两个长度均为N/2的结果序列；通过FFT算法将它们组合起来即可获得RealFFT变换结果。下文将提供该算法对应的伪代码实现：
3000'FFT FOR REAL SIGNALS
3010'Upon entry, N% contains the number of points in the DFT, REX[ ]contains
3020'the real input signal, while values in IMX[ ] are ignored. Uponreturn,
3030'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 toN%-1.
3040'
3050NH% = N%/2-1'Separate even and odd points
3060FOR I% = 0 TO NH%
3070REX(I%) = REX(2*I%)
3080IMX(I%) = REX(2*I%+1)
3090NEXT I%
3100'
3110N% = N%/2'Calculate N%/2 point FFT
3120GOSUB 1000
3130N% = N%*2
3140'
3150NM1% = N%-1'Even/odd frequency domain decomposition
3160ND2% = N%/2
3170N4% = N%/4-1
3180FOR I% = 1 TO N4%
3190IM% = ND2%-I%
3200IP2% = I%+ND2%
3210IPM% = IM%+ND2%
3220REX(IP2%) = (IMX(I%) + IMX(IM%))/2
3230REX(IPM%) = REX(IP2%)
3240IMX(IP2%) = -(REX(I%) - REX(IM%))/2
3250IMX(IPM%) = -IMX(IP2%)
3260REX(I%) = (REX(I%) + REX(IM%))/2
3270REX(IM%) = REX(I%)
3280IMX(I%) = (IMX(I%) - IMX(IM%))/2
3290IMX(IM%) = -IMX(I%)
3300NEXT I%
3310REX(N%*3/4) = IMX(N%/4)
3320REX(ND2%) = IMX(0)
3330IMX(N%*3/4) = 0
3340IMX(ND2%) = 0
3350IMX(N%/4) = 0
3360IMX(0) = 0
3370'
3380PI = 3.14159265'Complete the last FFT stage
3390L% = CINT(LOG(N%)/LOG(2))
3400LE% = CINT(2^L%)
3410LE2% = LE%/2
3420UR = 1
3430UI = 0
3440SR = COS(PI/LE2%)
3450SI = -SIN(PI/LE2%)
3460FOR J% = 1 TO LE2%
3470JM1% = J%-1
3480FOR I% = JM1% TO NM1% STEP LE%
3490IP% = I%+LE2%
3500TR = REX[IP%]*UR - IMX[IP%]*UI
3510TI = REX[IP%]*UI + IMX[IP%]*UR
3520REX[IP%] = REX[I%]-TR
3530IMX[IP%] = IMX[I%]-TI
3540REX[I%] = REX[I%]+TR
3550IMX[I%] = IMX[I%]+TI
3560NEXT I%
3570TR = UR
3580UR = TRSR - UISI
3590UI = TRSI + UISR
3600NEXT J%
3610RETURN