傳統(tǒng)的計(jì)算快速傅里葉變換的Cooley-Tukey算法效率極高，因其主要由蝶形運(yùn)算構(gòu)成，所以代碼形式也非常簡(jiǎn)單，只是需要將輸入或者輸出按照位反轉(zhuǎn)的方式重新排序。

這個(gè)重新排序的步驟并不是必須的。Clive Temperton于1991年在Self-Sorting In-Place Fast Fourier Transforms一文中給出了適用于混合基數(shù)的原地FFT算法，不需要對(duì)輸入或輸出重新排序。本文將介紹這種算法的原理并給出基數(shù)2（Radix-2）情況下的具體構(gòu)造和C++實(shí)現(xiàn)。作為FFT算法研究成果的集大成者，FFTW已應(yīng)用了這種算法。

FFT7071專欄?fourier.v.ariant.cn

預(yù)備：內(nèi)存中的矩陣轉(zhuǎn)置

設(shè)x[t]為一個(gè)長(zhǎng)度為M*N的向量，t也可表示為Ma+b。以M=5,N=3為例，x[0…14]={101, 102, 103,…, 115}，如果將x視作一個(gè)5行3列的矩陣，那么a列b行的矩陣元素即是x[Ma+b]：

a= 0, 1, 2 b=0 101 106 111 b=1 102 107 112 b=2 103 108 113 b=3 104 109 114 b=4 105 110 115

將這個(gè)矩陣轉(zhuǎn)置，不難發(fā)現(xiàn)轉(zhuǎn)置后的y[0…14]={101, 106, 111, 102,…, 110, 115}

a= 0, 1, 2, 3, 4 b=0 101 102 103 104 105 b=1 106 107 108 109 110 b=2 111 112 113 114 115

y與x的關(guān)系是y[Nb+a]=x[Ma+b]。這個(gè)轉(zhuǎn)置變換也可以用一個(gè)置換矩陣P表示：y=Px。

展開：DFT變換式

記長(zhǎng)度為

的信號(hào)為 ， 為 的離散傅里葉變換，并且設(shè)：

展開DFT變換得到以下表達(dá)式：

其中

， 。

利用

用y代換x并繼續(xù)展開單位根的冪，其中 ：

上式的求和等價(jià)于：

其中大括號(hào)內(nèi)的求和相當(dāng)于將y中的元素從地址0開始每相鄰的N個(gè)為一組總共M個(gè)長(zhǎng)度為N的DFT。如果將y看作M行N列的矩陣，這是對(duì)每一列的變換，變換的結(jié)果依次乘以

，這時(shí)剩下的一個(gè)求和相當(dāng)于對(duì)y的每一行單獨(dú)的DFT。

至此，長(zhǎng)度為MN的變換分解為了長(zhǎng)度為M和N的兩遍短變換。如果上式中不將第一遍DFT的結(jié)果乘以

，結(jié)果將是M*N的2維DFT。需要注意的是，變換y可以使上式的兩遍DFT均在原地進(jìn)行，如果變換的是x，為保持變換結(jié)果的順序正確，必須以轉(zhuǎn)置的形式寫回第一遍短變換的結(jié)果。

題外話：將中間結(jié)果寫到另一處存儲(chǔ)區(qū)x'，并且以x'為輸入做第二遍變換，結(jié)果寫回x，如此往復(fù)可以解決變換無法在原地進(jìn)行的問題，這即是Stockham FFT算法。但是如此一來FFT需要額外的等于x長(zhǎng)度的內(nèi)存，即需要額外O(N)的空間。因?yàn)橹脫Q群中的元素均可分解為2置換，也就是對(duì)換的乘積，置換矩陣也可做如此分解，將轉(zhuǎn)置操作變換成一系列對(duì)換，從而可在原地進(jìn)行，僅需要O(1)額外空間。然而轉(zhuǎn)置只能分解為數(shù)量巨大的對(duì)換，這種操作的效率不高。這也預(yù)示著，充分利用內(nèi)存中數(shù)據(jù)的對(duì)換，可以保持O(1)的額外空間需求，同時(shí)使FFT在原地進(jìn)行且順序正確。

展開：DFT變換矩陣和素因子算法

運(yùn)用上文得出的分解到DFT的矩陣形式：

實(shí)際變換的是y（也就是轉(zhuǎn)置的x），第一遍DFT作于相鄰的N個(gè)元素（每列），將結(jié)果逐個(gè)乘以一個(gè)旋轉(zhuǎn)因子，再變換間隔為N的每組元素（每行），這個(gè)過程對(duì)應(yīng)DFT矩陣的一種分解，也就是素因子分解FFT算法的基本構(gòu)造：

其中W為下標(biāo)相應(yīng)長(zhǎng)度的DFT矩陣；P為x[Ma+b]轉(zhuǎn)置到y(tǒng)[Nb+a]的置換矩陣；I為M或N階的單位矩陣；D為M*N階對(duì)角矩陣，第bN+d行對(duì)角線上的值為exp(2πibd/(MN))。?是矩陣的Kronecker積，定義如下：

例如：

Kronecker積滿足結(jié)合律：

滿足“混合乘積”性質(zhì)：

以

為例：

繼續(xù)分解

：

運(yùn)用“混合乘積”的性質(zhì)將上式拆分為矩陣積：

對(duì)于長(zhǎng)度為

的DFT矩陣分解，設(shè)：

可以得到：

這種分解正是時(shí)間抽取（DIT）基數(shù)2（Radix-2）的Cooley-Tukey算法，下文中只考慮此種FFT，頻率抽取（DIF）在附錄中討論。

已知T對(duì)應(yīng)Cooley-Tukey算法每次迭代在整個(gè)輸入向量上的所有蝶形運(yùn)算，上式中的

為2行8列到8行2列的矩陣轉(zhuǎn)置，作用是將x[2b+a]的值變換到x[8a+b]，其中 ， 。觀察8a+b和2b+a的二進(jìn)制表示：2^3 2^2 2^1 2^0 [b] [b] [b] [a] = 2b+a [a] [b] [b] [b] = 8a+b

可知

的作用是將 的地址t二進(jìn)制位向右環(huán)移1位。 的作用是保持t的最高1位不變，t的余下三位向右環(huán)移1位，因此經(jīng)過所有的 變換：2^3 2^2 2^1 2^0 [k3] [k2] [k1] [k0] [k0] [k3] [k2] [k1] - P16 [k0] [k1] [k3] [k2] - P8 [k0] [k1] [k2] [k3] - P4

很明顯T之前所有

矩陣的乘積是輸入數(shù)據(jù)翻轉(zhuǎn)對(duì)應(yīng)地址二進(jìn)制位的置換矩陣。

對(duì)換：蝶形運(yùn)算和對(duì)換

設(shè)DFT的長(zhǎng)度為N，則x[t]中的t在二進(jìn)制下有N位，定義

為對(duì)換 和 的置換矩陣， 由 對(duì)換二進(jìn)制位中低位i和高位N-i-1得到。可以用一系列Q的乘積取代P的乘積：

的效果同樣完全反轉(zhuǎn)二進(jìn)制位：2^3 2^2 2^1 2^0 [k3] [k2] [k1] [k0] [k0] [k2] [k1] [k3] - Q03 [k0] [k1] [k2] [k3] - Q12

在N=2M或2M+1的情況下：

轉(zhuǎn)置P無法簡(jiǎn)單地在原地計(jì)算，而Q僅包含數(shù)量較少的對(duì)換，因此可以在原地完成。目前為止以T和Q組成的DFT矩陣W分解仍然表示先重排數(shù)據(jù)再開始蝶形運(yùn)算的迭代，如果將T和Q結(jié)合起來，就能省略重排數(shù)據(jù)的操作，但是這要求T和Q可交換。為了證明這一點(diǎn)，首先需要求出Q的表達(dá)式。

觀察

翻轉(zhuǎn)二進(jìn)制最高位和最低位的操作：0000 0000 0001 -> 1000 0010 0010 0011 -> 1010 0100 0100 0101 -> 1100 0110 0110 0111 -> 1110 1000 -> 0001 1001 1001 1010 -> 0011 1011 1011 1100 -> 0101 1101 1101 1110 -> 0111 1111 1111

可以發(fā)現(xiàn)

的作用是將前一半數(shù)中的奇數(shù) 和 對(duì)換。因此：

這里

是 階置換矩陣，對(duì)換 和 ， 由 交換二進(jìn)制的最高位和最低位得到。

在為二進(jìn)制數(shù)添加前綴的操作下

的作用是不變的：00000 00000 10000 10000 00001 -> 01000 10001 -> 11000 00010 00010 10010 10010 00011 -> 01010 10011 -> 11010 00100 00100 10100 10100 00101 -> 01100 10101 -> 11100 00110 00110 10110 10110 00111 -> 01110 10111 -> 11110 01000 -> 00001 11000 -> 10001 01001 01001 11001 11001 01010 -> 00011 11010 -> 10011 01011 01011 11011 11011 01100 -> 00101 11100 -> 10101 01101 01101 11101 11101 01110 -> 00111 11110 -> 10111 01111 01111 11111 11111

因此對(duì)于N位二進(jìn)制數(shù)：

為二進(jìn)制數(shù)添加后綴的操作使

作用于全部后綴，可得出：

設(shè)

，現(xiàn)在可以將 展開：

因此對(duì)于所有的

， 和 可交換。這使 可以寫為：

已知一次蝶形運(yùn)算的迭代

的作用是：將兩個(gè)長(zhǎng)度為 的DFT結(jié)果作為奇偶兩部分，合并為長(zhǎng)度為 的DFT結(jié)果，這樣的奇偶對(duì)共有 個(gè)。

令

中單個(gè)蝶形運(yùn)算的偶、奇兩個(gè)輸入分別是 和 ， 會(huì)將 與 對(duì)換。取 ，則 將與 對(duì)換。在 的情況下， 與 分別是另一蝶形運(yùn)算的偶、奇輸入。

可見

中，輸入數(shù)據(jù)地址相差 的兩個(gè)蝶形運(yùn)算是成對(duì)的，第一個(gè)蝶形運(yùn)算的奇數(shù)項(xiàng)輸入與第二個(gè)蝶形運(yùn)算的偶數(shù)項(xiàng)輸入對(duì)換。如下圖所示：

下圖中畫出來N=16的基數(shù)2變換，輸入和輸出的順序均是正確的；圖中用顏色標(biāo)出了某些蝶形運(yùn)算，使蝶形運(yùn)算的配對(duì)更清晰。注意其中成對(duì)蝶形運(yùn)算的4個(gè)輸入輸出均在原地，并且與傳統(tǒng)Cooley-Tukey算法相比沒有計(jì)算量的差別。

下圖是作為對(duì)比的傳統(tǒng)Cooley-Tukey算法。

附錄：本文算法的C++實(shí)現(xiàn)

/* copyright 2020, github.com/zhxxch, all rights reserved. */ #include <complex> #include <iterator> #include <cmath> #include <cassert> #include <cstddef> template<typename iter_t> #if __cplusplus > 201703L requires std::random_access_iterator<iter_t> #endifinline void fft_in_place(iter_t arr_begin,iter_t arr_end, const bool is_inverse) {using cplx_t = typename std::iterator_traits<iter_t>::value_type;using real_t = typename cplx_t::value_type;const size_t length= std::distance(arr_begin, arr_end);assert("requires length = pow(2,n)"&& (length & (length - 1)) == 0);constexpr real_t pi= 3.141592653589793238462643383;size_t sub_ft_size = 1;size_t num_sub_ft = length / sub_ft_size;size_t num_sub_ft_pair = num_sub_ft / 2;for(; sub_ft_size < num_sub_ft_pair;sub_ft_size *= 2, num_sub_ft /= 2,num_sub_ft_pair /= 2) {for(size_t coupled_group_pos = 0;coupled_group_pos < length;coupled_group_pos += 2 * num_sub_ft_pair) {for(size_t sub_ft_pos = coupled_group_pos;sub_ft_pos< coupled_group_pos + num_sub_ft_pair;sub_ft_pos += 2 * sub_ft_size) {for(size_t i = sub_ft_pos, nth_pow = 0;i < sub_ft_pos + sub_ft_size;i++, nth_pow += num_sub_ft_pair) {const cplx_t W = exp((is_inverse ? 1. : -1.)* cplx_t(0, 2 * nth_pow * pi)/ (real_t)length);const cplx_t parit00= arr_begin[i];const cplx_t parit01= arr_begin[num_sub_ft_pair+ i]* W;const cplx_t parit10= arr_begin[i + sub_ft_size];const cplx_t parit11= arr_begin[num_sub_ft_pair + i+ sub_ft_size]* W;arr_begin[i] = parit00 + parit01;arr_begin[i + sub_ft_size]= parit00 - parit01;arr_begin[num_sub_ft_pair + i]= parit10 + parit11;arr_begin[num_sub_ft_pair + i+ sub_ft_size]= parit10 - parit11;}}}}for(; sub_ft_size < length; sub_ft_size *= 2,num_sub_ft /= 2, num_sub_ft_pair /= 2) {for(size_t sub_ft_pos = 0; sub_ft_pos < length;sub_ft_pos += 2 * sub_ft_size) {for(size_t i = sub_ft_pos, nth_pow = 0;i < sub_ft_pos + sub_ft_size;i++, nth_pow += num_sub_ft_pair) {const cplx_t parit1= arr_begin[i + sub_ft_size]* exp((is_inverse ? 1. : -1.)* cplx_t(0, 2 * nth_pow * pi)/ (real_t)length);const cplx_t parit0 = arr_begin[i];arr_begin[i] = parit0 + parit1;arr_begin[i + sub_ft_size]= parit0 - parit1;}}} }

使用方法-FFT

fft_in_place(arr.begin(), arr.end(), 0);

使用方法-IFFT

fft_in_place(arr.begin(), arr.end(), 1);

arr的長(zhǎng)度必須是2的冪。

附錄：時(shí)間抽取與頻率抽取

以

為例，頻率抽取的FFT算法中：

換為使用Q表達(dá)的形式則為：

因此Q作用于蝶形運(yùn)算的輸出。

總結(jié)

以上是生活随笔為你收集整理的基4fft算法的蝶形图_原地且自动整序的FFT算法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。