數(shù)字信號(hào)處理離散信號(hào)處理：第9章 離散傅里葉變換的計(jì)算

1、第9章 離散傅里葉變換的計(jì)算9.0 引言9.1 離散傅里變換的高效計(jì)算9.2 Goertzel算法9.3 按時(shí)間抽取的FFT算法9.4 按頻率抽取的FFT算法9.5 小結(jié)9.0 引言降低離散傅里葉變換的計(jì)算量方法：利用離散傅里葉變換系數(shù)的對(duì)稱性和周期性主要算法：Goertzel算法按時(shí)間抽取的FFT算法按頻率抽取的FFT算法里程碑：1965年 Cooley Tukery9.1 離散傅里變換的高效計(jì)算N點(diǎn)有限長(zhǎng)序列的DFT：反變換：計(jì)算復(fù)乘復(fù)加實(shí)乘實(shí)加復(fù)乘142復(fù)加12XkNN-14N4N-2DFT4N24N2+2NFFT直接計(jì)算：傅里葉變換系數(shù)的性質(zhì)對(duì)稱性（復(fù)共軛對(duì)稱）周期性（n,k的周期性）

2、利用對(duì)稱性可直接計(jì)算n 和 N-n乘法減少1/2，利用對(duì)稱性可大量降低計(jì)算量離散傅里變換快速算法的尋找過(guò)程1805年的高斯Runge(1905),Danielson & Lanczos(1942)指出DFT的計(jì)算量正比于N log N 而不是NN。1965年Cooley Tukery的論文指出，當(dāng)N為復(fù)合數(shù)時(shí)，可以進(jìn)行分解。該論文為DFT的快速算法的發(fā)現(xiàn)指出了一條行之有效的方法。此后出現(xiàn)了一大批快速算法：按時(shí)間抽取的FFT算法、按頻率抽取的FFT算法、素因子法等等。9.2 Goertzel算法利用了周期性定義一個(gè)中間序列計(jì)算復(fù)乘復(fù)加實(shí)乘實(shí)加復(fù)乘142復(fù)加12XkNN-14N4N-21144極點(diǎn)

3、：2次實(shí)乘 4次實(shí)加（每個(gè)樣本的計(jì)算量）2N次實(shí)乘 4N次實(shí)加零點(diǎn)：迭代的最后一次計(jì)算 4次實(shí)乘 4次實(shí)加計(jì)算量極點(diǎn)：2次實(shí)乘 4次實(shí)加（每個(gè)樣本的計(jì)算量）2N次實(shí)乘 4N次實(shí)加零點(diǎn)：迭代的最后一次計(jì)算 4次實(shí)乘 4次實(shí)加共需：2N2+4N次實(shí)乘 4N2+4N次實(shí)加利用共軛對(duì)稱一次可以算出兩個(gè)點(diǎn)的極點(diǎn)，零點(diǎn)復(fù)共軛。所以一次可以算出兩個(gè)點(diǎn)。N2+2N次實(shí)乘 2N2+2N次實(shí)加可以計(jì)算N點(diǎn)DFT其中的任意長(zhǎng)（M）其計(jì)算量MN9.3 按時(shí)間抽取的FFT算法2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT兩個(gè)2點(diǎn)DFT兩個(gè)2點(diǎn)DFT兩個(gè)2點(diǎn)DFT兩個(gè)2點(diǎn)DFT兩個(gè)4點(diǎn)D

4、FT兩個(gè)4點(diǎn)DFT兩個(gè)N/2點(diǎn)DFTX1(k).X2(k)x(n)一.DFT的計(jì)算工作量?jī)烧叩牟顒e僅在指數(shù)的符號(hào)和因子1/N. 通常x(n)和 都是復(fù)數(shù),所以計(jì)算一個(gè)X(k)的值需要N次復(fù)數(shù)乘法運(yùn)算,和 次復(fù)數(shù)加法運(yùn)算.那么,所有的X(k)就要N2次復(fù)數(shù)乘法運(yùn)算,N(N-1)次復(fù)數(shù)加法運(yùn)算.當(dāng)N很大時(shí),運(yùn)算量將是驚人的,如N=1024,則要完成1048576 次(一百多萬(wàn)次)運(yùn)算.這樣,難以做到實(shí)時(shí)處理.一個(gè)X(k)的值的工作量,如X(1)二.改進(jìn)的途徑 1. 的對(duì)稱性和周期性得:對(duì)稱性:周期性: 利用上述特性,可以將有些項(xiàng)合并,并將DFT分解為短序列,從而降低運(yùn)算次數(shù),提高運(yùn)算速度.1965

5、年,庫(kù)利(cooley)和圖基(Tukey)首先提出FFT算法.對(duì)于N點(diǎn)DFT,僅需(N/2)log2N 次復(fù)數(shù)乘法運(yùn)算.例如N=1024=210 時(shí)，需要(1024/2)log2 210 =512*10=5120次。5120/1048576=4.88% ,速度提高20倍1.先將 按n的奇偶分為兩組作DFT,設(shè)N=2L ,不足時(shí),可補(bǔ)些零。這樣有:按時(shí)間抽取(DIT)的FFT算法 庫(kù)利-圖基算法一.算法原理(基2FFT)(一)N/2點(diǎn)DFT因此，n為偶數(shù)時(shí):n為奇數(shù)時(shí):由于: 所以,上式可表示為:(n為偶數(shù)) (n為奇數(shù)) x1(0)=x(0) x1(1)=x(2) N/2點(diǎn) x1(2)=x(

6、4) DFT x1(3)=x(6) x2(0)=x(1) x2(1)=x(3) N/2點(diǎn) x2(2)=x(5) DFT x2(3)=x(7) X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN2WN1WN0WN3X(0)X(1)X(2)X(3)由于 （周期性）,所以: 同理, 這就是說(shuō),X1(k),X2(k)的后一半,分別 等于其前一半的值。3.X(k)的后一半的確定又由于 ,所以 可見(jiàn),X(k)的后一半，也完全由X1(k),X2(k）的前一半所確定。*N點(diǎn)的DFT可由兩個(gè)N/2點(diǎn)的DFT來(lái)計(jì)算。 x1(0)=x(0) x1(1)=x(2) N/2點(diǎn) x1(2

7、)=x(4) DFT x1(3)=x(6) x2(0)=x(1) x2(1)=x(3) N/2點(diǎn) x2(2)=x(5) DFT x2(3)=x(7) X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN2WN1WN0WN3-1-1-1-1X(4)X(5)X(6)X(7)實(shí)現(xiàn)上式運(yùn)算的流圖稱作蝶形運(yùn)算4.蝶形運(yùn)算（N/2個(gè)蝶形) 前一半 后一半(前一半)(后一半)1 1 11-1由X1(k)、X 2(k)表示X(k)的運(yùn)算是一種特殊的運(yùn)算- 碟形運(yùn)算 x1(0)=x(0) x1(1)=x(2) N/2點(diǎn) x1(2)=x(4) DFT x1(3)=x(6) x2(0

8、)=x(1) x2(1)=x(3) N/2點(diǎn) x2(2)=x(5) DFT x2(3)=x(7) X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN2WN1WN0WN3-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)(3)N/2個(gè)蝶形運(yùn)算的運(yùn)算量:復(fù)乘次數(shù):復(fù)加次數(shù):(1)N/2點(diǎn)的DFT運(yùn)算量:復(fù)乘次數(shù):復(fù)加次數(shù):5.計(jì)算工作量分析按奇、偶分組后的計(jì)算量：(2)兩個(gè)N/2點(diǎn)的DFT運(yùn)算量:復(fù)乘次數(shù):復(fù)加次數(shù):復(fù)加:復(fù)乘:總共運(yùn)算量: *但是,N點(diǎn)DFT的復(fù)乘為N2;復(fù)加N(N-1);與分解后相比可知,計(jì)算工作點(diǎn)差不多減少一半。

9、 由于N=2 L ,所以 N/2仍為偶數(shù),可以進(jìn)一步把每個(gè)N/2點(diǎn)的序列再按其奇偶部分，分解為兩個(gè)N/4的子序列。例如，n為偶數(shù)時(shí)的 N/2點(diǎn)分解為:(二) N/4點(diǎn)DFT進(jìn)行N/4點(diǎn)的DFT,得到(偶中偶)(偶中奇)從而可得到前N/4的X1(k)后N/4的X1(k)為(奇中偶)(奇中奇)同樣對(duì)n為奇數(shù)時(shí)，N/2點(diǎn)分為兩個(gè)N/4點(diǎn) 的序列進(jìn)行DFT，則有: 例如,N=8時(shí)的DFT可分解為四個(gè)N/4的DFT, 具體步驟如下:(1) 將原序列x(n)的“偶中偶”部分:構(gòu)成N/4點(diǎn)DFT，從而得到X3(0), X3(1)。(2) 將原序列x(n)的“偶中奇”部分:構(gòu)成N/4點(diǎn)DFT，從而得到X4(0

10、), X4(1)。(3) 將原序列x(n)的“奇中偶”部分:構(gòu)成N/4點(diǎn)DFT，從而得到X5(0), X5(1)。（4）將原序列x(n)的“奇中奇”部分:構(gòu)成N/4點(diǎn)DFT，從而得到X6(0), X6(1)。（5）由X3(0),X3(1),X4(0),X4(1)進(jìn)行碟形運(yùn)算,得到X1(0),X1(1),X1(2),X1(3)。（6）由X5(0),X5(1),X6(0),X6(1)進(jìn)行碟形運(yùn)算,得到X2(0),X2(1),X2(2),X2(3)。 (0)= (0)= (0) N/4 (1)= (2)= (4) DFT (0)= (1)= (2) N/4 (1)= (3)= (6) DFT (0)

11、= (0)= (1) N/4 (1)= (2)= (5) DFT (0)= (1)= (3) N/4 (1)= (3)= (7) DFT3 13 14 14 15 25 26 26 2 02NN02NNWWWW0123NNNN-1-1-1-2-1-1WWWWX (0)X (1)X (0)X (1)X (0)X (1)X (0)X (1)33445566X (0)X (1)X (2)X (3)X (0)X (1)X (2)X (3)11122221X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)（7）由X1(0),X1(1),X1(2),X1(3),X2(0), X2(1),X2

12、(2),X2(3)進(jìn)行碟形運(yùn)算,得到X(0), X(1),X(2),X(3),X(4),X(5),X(6),X(7)。 這樣,又一次分解,得到四個(gè)N/4點(diǎn)DFT, 兩級(jí)蝶形運(yùn)算,其運(yùn)算量有大約減少一半 即為N點(diǎn)DFT的1/4。 對(duì)于N=8時(shí)DFT,N/4點(diǎn)即為兩點(diǎn)DFT,因此 亦即, (0) (4) (2) (6) (1) (5) (3) (7)WN0WN0WN0W0N-1-1-1-1X (0)X (1)X (0)X (1)X (0)X (1)X (0)X (1)33445566WN0WN2WN0WN2-1-1-1-1X (0)X (1)X (2)X (3)X (0)X (1)X (2)X (

13、3)11121222WWWWN0N1N2N3-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)因此,8點(diǎn)DFT的FFT的運(yùn)算流圖如下這種FFT算法,是在時(shí)間上對(duì)輸入序列的次序是屬于偶數(shù)還是屬于奇數(shù)來(lái)進(jìn)行分解的,所以稱作按時(shí)間抽取的算法(DIT)。 二.運(yùn)算量 由上述分析可知,N=8需三級(jí)蝶形運(yùn)算 N=2 =8,由此可知,N=2L 共需L級(jí)蝶形運(yùn)算,而且每級(jí)都由N/2個(gè)蝶形運(yùn)算 組成,每個(gè)蝶形運(yùn)算有一次復(fù)乘,兩次復(fù)加。 因此,N點(diǎn)的FFT的運(yùn)算量為復(fù)乘: mF =（N/2）L=（N/2）log2 N復(fù)加: aF =N L=N log2 N 由于計(jì)算機(jī)的乘法運(yùn)算比加法

14、運(yùn)算所需的時(shí)間多得多，故以乘法作為比較基準(zhǔn)。 (0)=X0(0) X1(0) X2(0) X3(0)=X(0) (4)=X0(1) X1(1) X2(1) X3(1)=X(1) (2)=X0(2) X3(2)=X(2) (6)=X0(3) X3(3)=X(3) (1)=X0(4) X1(4) X2(4) X3(4)=X(4) (5)=X0(5) X3(5)=X(5) (3)=X0(6) X3(6)=X(6) (7)=X0(7) X1(7) X2(7) X3(7)=X(7)輸入數(shù)據(jù)、中間運(yùn)算結(jié)果和最后輸出均用同一存儲(chǔ)器。WWWWN0N0N0N0-1-1-1-1WWWWN0N2N0N2-1-1-1

15、-1WWWWNNNN0同址計(jì)算 設(shè)用m(m=1,2, ,L)表示第m列;用k,j表示蝶形輸入數(shù)據(jù)所在的（上/下）行數(shù)(0,1,2, ,N-1);這時(shí)任何一個(gè)蝶形運(yùn)算可用下面通用式表示， 即 由運(yùn)算流圖可知,一共有N個(gè)輸入/出行,一共有l(wèi)og2 N=L列（級(jí)）蝶形運(yùn)算(基本迭代運(yùn)算). 所以,當(dāng)m=1時(shí),則有（前兩個(gè)蝶形） 當(dāng)m=2時(shí),則有（前兩個(gè)蝶形） 當(dāng)m=3時(shí),則有（前兩個(gè)蝶形） 可見(jiàn),在某列進(jìn)行蝶形運(yùn)算的任意兩個(gè)節(jié)點(diǎn)(行)k和j的節(jié)點(diǎn)變量 就完全可以確定蝶形運(yùn)算的結(jié)果 ,與其它行(節(jié)點(diǎn))無(wú)關(guān)。 這樣,蝶形運(yùn)算的兩個(gè)輸出值仍可放回蝶形運(yùn)算的兩個(gè)輸入所在的存儲(chǔ)器中,即實(shí)現(xiàn)所

16、謂原位運(yùn)算。每一組(列)有N/2個(gè)蝶形運(yùn)算,所以只需N個(gè)存儲(chǔ)單元，可以節(jié) 省存儲(chǔ)單元。 2 倒位序規(guī)律 由圖可知,輸出X(k)按正常順序排列在存儲(chǔ)單元,而輸入是按順序:這種順序稱作倒位序,即二進(jìn)制數(shù)倒位。n =00n =10n =01n =11n =01n =1101010101 (n2)x(000) 0 x(100) 4x(010) 2x(110) 6x(001) 1x(101) 5x(011) 3x(111) 7（偶）（奇）這是由奇偶分組造成的,以N=8為例說(shuō)明如下: 3.倒位序?qū)崿F(xiàn) 輸入序列先按自然順序存入存儲(chǔ)單元,然后經(jīng)變址運(yùn)算來(lái)實(shí)現(xiàn)倒位序排列設(shè)輸入序列的序號(hào)為n,二進(jìn)制為(n2 n1

17、 n0 )2 ,倒位序順序用 表示,其倒位序二進(jìn)制為(n0 n1 n2 )2 ,例如，N=8時(shí)如下表： 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 4 2 0 1 0 0 1 0 2 3 0 1 1 1 1 0 6 4 1 0 0 0 0 1 1 5 1 0 1 1 0 1 5 6 1 1 0 0 1 1 3 7 1 1 1 1 1 1 7 自然順序n 二進(jìn)制n n n 倒位序二進(jìn)制n n n 倒位順序n2 1 0 0 1 2變址處理方法A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(8)x(0) x(1) x(2) x(3) x(4) x(5) x(

18、6) x(7)x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)存儲(chǔ)單元自然順序變址倒位序4.蝶形運(yùn)算兩節(jié)點(diǎn)的距離:2m-1 其中,m表示第m列,且m =1, ,L 例如N=8=23 ,第一級(jí)(列)距離為21-1=1, 第二級(jí)(列)距離為22-1=2， 第三級(jí)(列)距離為23-1=4。 (0) (4) (2) (6) (1) (5) (3) (7)WN0WN0WN0W0N-1-1-1-1X (0)X (1)X (0)X (1)X (0)X (1)X (0)X (1)33445566WN0WN2WN0WN2-1-1-1-1X (0)X (1)X (2)X (3)X (

19、0)X (1)X (2)X (3)11121222WWWWN0N1N2N3-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)因此,8點(diǎn)DFT的FFT的運(yùn)算流圖如下 考慮蝶形運(yùn)算兩節(jié)點(diǎn)的距離為2m-1,蝶形運(yùn)算可表為 Xm(k)=Xm-1(k)+Xm-1(k+2m-1)WNr Xm(k+2m-1)=Xm-1(k)-Xm-1(k+2m-1)WNr 由于N為已知,所以將r的值確定即可。5.WNr 的確定(僅給出方法) 為此,令k=(n2n1n0)2,再將k=(n2n1n0)2 左移(L-m)位,右邊位置補(bǔ)零,就可得到(r)2的值，即(r)2 =(k)22L-m 。 例如

20、 N=8=23 (1)k=2 , m=3 的r值,k=2=(010)2 左移L-m=3-3=0 , r=(010)2=2; (2)k=3 ,m=3 的r值; k= 3=(011)2,左移0位,r=3； (3)k=5 ,m=2的值; k=5=(101)2 左移L-m=1位 r=(010)2 =2。 6.存儲(chǔ)單元 存輸入序列 (n),n=0,1, ,N-1,計(jì)N 個(gè)單元; 存放系數(shù) ,r=0,1, ,（N/2）-1,需N/2個(gè)存儲(chǔ)單元； 共計(jì)(N+N/2)個(gè)存儲(chǔ)單元。9.3.2其他形式它與輸入是亂序，輸出順序比較，看出：相同點(diǎn)：運(yùn)算量一樣；不同點(diǎn)：第一是數(shù)據(jù)存入方式不同；第二是取用系數(shù)的順序不同。

21、x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(4)X(2)X(6)X(1)X(3)X(5)X(7)（1）它失掉了“原位運(yùn)算”的性質(zhì)。（2）為了變換N點(diǎn)數(shù)據(jù)，至少需要2N個(gè)復(fù)數(shù)單元。在內(nèi)存比較緊張時(shí)，而對(duì)輸入數(shù)據(jù)整序并不困難時(shí)一般不用。為了爭(zhēng)取速度，才采取這種變體。x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)DIF的FFT算法(桑德圖基算法) 一.算法原理 1.N點(diǎn)DFT的另一種表達(dá)形式9.4 按頻率抽取的FFT算法由于 故因此 X(k)可表為 2.N點(diǎn)DFT按k的奇偶分組可分為兩個(gè)

22、N/2的DFT 當(dāng)k為偶數(shù),即 k=2r時(shí),(-1)k =1； 當(dāng)k為奇數(shù),即 k=2r+1 時(shí) (-1)k =-1 。 這時(shí)X(k)可分為兩部分： k為偶數(shù)時(shí)： 可見(jiàn),上面兩式均為N/2的DFT。k為奇數(shù)時(shí)：3.蝶形運(yùn)算-14.N=8時(shí),按k的奇偶分解過(guò)程 先蝶形運(yùn)算，后DFT:-1-1-1-1WWWWNNNN0123 5.仿照DIT的方法 再將N/2點(diǎn)DFT按k的奇偶分解為兩個(gè)N/4點(diǎn)的DFT,如此進(jìn)行下去,直至分解為2點(diǎn)DFT。 (0) X(0) (1) X(4) (2) X(2) (3) X(6) (4) X(1) (5) X(5) (6) X(3) (7) X(7)-1-1-1-1W

23、WWWNNNN0123-1-1-1-1WWWWNNNN0202-1-1-1-1WWWWNNNN0000例如 N=8時(shí)DIF的FFT流圖如下：二.原位運(yùn)算 每級(jí)(列)都是由N/2個(gè)蝶形運(yùn)算構(gòu)成,即 -1WNr三.蝶形運(yùn)算兩節(jié)點(diǎn)的距離 一般公式為2L-m =N/2m 例如 N=23 =8 ： (1)m=1 時(shí)的距離為 8/2=4； (2)m=2 時(shí)的距離為 8/4=2; (3)m=3 時(shí)的距離為 8/8=1。四. 的計(jì)算 由于DIF蝶形運(yùn)算的兩節(jié)點(diǎn)的距離為 N/2m ,所以蝶形運(yùn)算可表為：r的求法： k=(n2 n1 n0 ) ,左移m-1位,右邊空出補(bǔ)零,得(r)2 ,亦即(r)2 =(k)2 2m-1 .例如,N=8:(1)m=1,k=2, k=(010)2 左移0位,(r)2=(010)2=2;(2)m=2,k=1, k=(001)2 左移1位,(r)2=(010)2=2;(3)m=2,k=5, k=(101)2 左移1位,(r)2=(010)2=2 . 1.相同點(diǎn) (1)進(jìn)行原位運(yùn)算； (2)運(yùn)算量相同,均為（N/2) Log2N 次復(fù)乘,N Log2N次復(fù)加。 2.不同點(diǎn) (1)DIT輸入為倒位序,輸出為自然順序；DIF正