二維FFT程序?qū)崿F(xiàn)及其應(yīng)用

1、葉東明：二維FFT的程序?qū)崿F(xiàn)及其應(yīng)用二維FFT的程序?qū)崿F(xiàn)及其應(yīng)用Donjon理學(xué)院信息與計(jì)算科學(xué)專業(yè)2002屆2班摘要譜分析方法（Fourier變換方法）中，由逼近方法導(dǎo)出的代數(shù)方程組中系數(shù)矩陣基本上是滿的，求解計(jì)算量太大。進(jìn)一步分析發(fā)現(xiàn)求解一個(gè)N點(diǎn)數(shù)據(jù)離散Fourier變換的復(fù)數(shù)乘法計(jì)算量正比于。而FFT算法對(duì)于長(zhǎng)度為N的復(fù)序列，它的復(fù)數(shù)乘法計(jì)算量為,因此計(jì)算量的節(jié)省是巨大的。FFT算法的出現(xiàn)Fourier變換的研究和應(yīng)用的面貌出現(xiàn)根本轉(zhuǎn)變。人們開始重新考慮它的優(yōu)點(diǎn)，越來(lái)越多地將它用于解決各種各樣的實(shí)際問(wèn)題。本文分析了一維FFT算法及二維FFT算法，并給出了相應(yīng)的程序和實(shí)例。關(guān)鍵詞FFT程序

2、Fourier變換 2D-DFT 行列算法Contributing the program of fast Fourier transform for dimensional and its applicationAbstractThe spectrum analysis method (Fourier transformation method) center, by approaches in the algebra system of equations which the method derives the coefficient matrix basically is full,

3、the computation of solution is large. Further analyzes the discovery to solve N points according to be separated the Fourier transformation complex multiplication computation in proportion to. But the FFT algorithm regarding the length is the N duplicate sequence, its complex multiplication computat

4、ion is,Therefore computation saving is huge. FFT the algorithm appearance Fourier transformation research and the application appearance appears the radical transformation. The people start to reconsider its merit, are more and more many uses in it solving various actual problem. This article has an

5、alyzed the unidimensional FFT algorithm and the two-dimensional FFT algorithm, and has produced the corresponding procedure and the example.Keywords FFT program Fourier transform 2D-DFT Row arithmetic目錄0引言.11 預(yù)備知識(shí).11.1三角級(jí)數(shù)及其正交性.11.2周期函數(shù)的Fourier級(jí)數(shù).11.3Fourier積分.21.4Fourier變換.41.5離散Fourier變換.51.6二維Fou

6、rier變換.71.7二維離散Fourier變換.82DFT快速計(jì)算(FFT)算法分析及其程序?qū)嵗?82.1DFT的快速算法.82.2基2FFT算法分析.92.3復(fù)序列基2FFT算法.122.4二維復(fù)序列2DDFT的行列算法.14結(jié)論.20致謝.20參考文獻(xiàn).210引言變換常常可以簡(jiǎn)化問(wèn)題的分析和求解過(guò)程，在科學(xué)研究的許多領(lǐng)域，人們發(fā)現(xiàn)Fourier變換對(duì)于問(wèn)題的求解和簡(jiǎn)化特別有用。Fourier變換可以看作是時(shí)間域上的函數(shù)在頻率域上的表示，譜函數(shù)（Fourier變換）在頻率域中包含的信息和原時(shí)間域所包含的信息是等價(jià)的，不同的僅是信息的表述方式。譜分析方法的基本思想來(lái)源于經(jīng)典的Ritz-Gal

7、erkin方法。由于許多的問(wèn)題可應(yīng)用Fourier變換，因此計(jì)算數(shù)學(xué)很自然地將Fourier變換離散化產(chǎn)生所謂離散Fourier變換人（DFT），然后使用計(jì)算機(jī)求解。在譜分析方法中，大多數(shù)情況下譜函數(shù)在每一點(diǎn)上的值都與原函數(shù)在所有點(diǎn)上的值有聯(lián)系，因此由逼迫方法導(dǎo)出的代數(shù)方程組中系數(shù)矩陣基本上是滿的，求解計(jì)算量太大，即便使用高速計(jì)算機(jī)，其運(yùn)算時(shí)間還是太長(zhǎng)，是不可實(shí)現(xiàn)的。因此在相當(dāng)長(zhǎng)的時(shí)間內(nèi)，F(xiàn)ourier變換方法沒(méi)有得到應(yīng)有的重視、發(fā)展和應(yīng)用。人們不懈地尋找減少Fourier變換計(jì)算量的方法和技術(shù)，直至1965年，Cooler-Tukey發(fā)表的算法，才得到人們的認(rèn)同。人們稱之為“快速Fourie

8、r變換（FFT）”方法。FFT算法有效地解決了Fourier變換計(jì)算量太大的問(wèn)題，它的出現(xiàn)使得Fourier變換的研究和應(yīng)用的面貌出現(xiàn)根本轉(zhuǎn)變。人們開始重新考慮它的優(yōu)點(diǎn)，越來(lái)越多地將它用于解決各種各樣的實(shí)際問(wèn)題。本文首先給出了FFT用到一些基本理論知識(shí)，接著對(duì)一維及二維FFT算法進(jìn)行了分析，并給出了相應(yīng)的程序及其實(shí)例。1 預(yù)備知識(shí)1.1三角級(jí)數(shù)及其正交性設(shè)為任意實(shí)數(shù)，是長(zhǎng)度為的區(qū)間。顯然三角函數(shù)，為正整數(shù)，皆為周期為的周期函數(shù)。由計(jì)算易得： (1-1-1)由此三角函數(shù)系 中任意兩個(gè)不同函數(shù)的乘積在上的積分等于0。我們稱這個(gè)函數(shù)系在長(zhǎng)為區(qū)間上具有正交性。1.2周期函數(shù)的Fourier級(jí)數(shù)定義1.

9、2.1設(shè)以為周期的函數(shù)在上可積和絕對(duì)可積，則令：(1-2-1)從而得到三角級(jí)數(shù) (1-2-2)我們稱這個(gè)三角級(jí)數(shù)為關(guān)于三角函數(shù)系的Fourier級(jí)數(shù)，稱，為的Fourier系數(shù)，記為： (1-2-3)當(dāng)是周期為的可積和絕對(duì)可積的連續(xù)函數(shù)時(shí)，可展開為收斂的Fourier級(jí)數(shù)（其證明參考文獻(xiàn)8）:利用Euler公式由此得這時(shí)Fourier級(jí)數(shù)可改變成復(fù)數(shù)形式： (1-2-4)顯然 1.3Fourier積分假設(shè)在上可積和絕對(duì)可積，且在任何上展開為Fourier級(jí)數(shù)。 (1-3-1)其中 (1-3-2)這時(shí) (1-3-3)當(dāng)時(shí) 記這時(shí)(1-3-3)右端第2項(xiàng)為：其中當(dāng)時(shí)，形式上有 (1-3-4)上式可

10、以改寫為其中 (1-3-5)在(1-3-4)中，令 (1-3-6)由是連續(xù)函數(shù)（證明見(jiàn)參考文獻(xiàn)8），積分 (1-3-7)存在。定義1.3.1當(dāng)時(shí)，由(1-3-7)得到的積分公式稱為的Fourier積分，記為 (1-3-8)在的Fourier積分中，由于是的偶函數(shù)，由(1-3-8)得 (1-3-9)以積分存在，且為的奇函數(shù)，所以令，其Cauchy積分主值為0，即因此(1-3-10)如果在點(diǎn)滿足Fourier積分收斂條件，則(1-3-10)右端等于。由此得Fourier積分的復(fù)數(shù)形式(1-3-11)類似地，我們還可以得到Fourier積分另外一種復(fù)數(shù)形式(1-3-12)1.4Fourier變換應(yīng)用

11、Fourier積分公式解決實(shí)際問(wèn)題時(shí)，常把它改寫成積分變換的形式。定義1.4.1若實(shí)變量的復(fù)值函數(shù) (1-4-1)在上可積和絕對(duì)可積，其中為實(shí)值函數(shù)，則含參變量的積分 (1-4-2)把函數(shù)變?yōu)楹瘮?shù)，稱為Fourier正變換，簡(jiǎn)稱為Fourier變換。稱為的Fourier變換或譜函數(shù)，稱為的Fourier逆變換或原函數(shù)。定義1.4.2對(duì)應(yīng)Fourier變換(1-4-2)由Fourier積分可得 (1-4-3)稱(1-4-3)為Fourier變換(1-4-2)的逆變換或稱為反演公式。逆變換使我們可以從一個(gè)時(shí)間函數(shù)的Fourier變換來(lái)確定這個(gè)時(shí)間函數(shù)。當(dāng)函數(shù)和函數(shù)由公式(1-4-2) 、(1-4-

12、3)相聯(lián)系則稱函數(shù)和為Fourier變換對(duì)，記為由Fourier積分(1-3-11)可以將Fourier變換及逆變換作為變換對(duì)寫成較為一般的形式： (1-4-4)其中，滿足，取+1或-1。在許多資料中還直接規(guī)定，取，這種取法雖然不會(huì)產(chǎn)生本質(zhì)錯(cuò)誤，但當(dāng)正逆變換配對(duì)時(shí)會(huì)引起一些混淆。解決這個(gè)矛盾的一個(gè)比較合理的方案是定義一個(gè)新的變量，表示頻率，F(xiàn)ourier變換對(duì)取為 (1-4-5)1.5離散Fourier變換如果想用計(jì)算機(jī)來(lái)計(jì)算和處理Fourier變換就需要將Fourier變換離散化。若自變量是定義在時(shí)間軸上的連續(xù)變量，我們則稱是連續(xù)時(shí)間函數(shù)。若自變量?jī)H在時(shí)間軸的離散點(diǎn)上取值，則稱為離散時(shí)間函數(shù)

13、，或稱為離散時(shí)間序列,即定義1.5.1對(duì)于任一屬于空間的離散時(shí)間序列稱 (1-5-1)為的離散時(shí)間Fourier變換DTFT(discrete time Fourier transform)。為的函數(shù)。由于故是周期為周期函數(shù)。對(duì)應(yīng)(1-5-1)可得DTFT逆變換公式為 (1-5-2)離散時(shí)間Fourier正逆變換也可以定義為 (1-5-3)其中是以1為周期的周期函數(shù)。DTFT對(duì)也記為。對(duì)于Fourier變換 (1-5-4)取充分大的有限區(qū)間，使得的Fourier積分在外充分小，以至可以近似地認(rèn)為將2N等分，則子區(qū)間長(zhǎng)度和分點(diǎn)為 (1-5-5)這時(shí)在每個(gè)小區(qū)間若用左矩形求積公式計(jì)算數(shù)值積分，并再

14、將離散化取 (1-5-6)這時(shí) (1-5-7)記 (1-5-8)則得到(1-5-9)利用數(shù)值積分的復(fù)化左矩形公式得到的式(1-5-9)給出了序列和序列之間的關(guān)系。定義1.5.2對(duì)長(zhǎng)度為N的復(fù)序列稱 (1-5-10)為序列的離散Fourier變換DFT(discrete Fourier transform)，或稱為DFT正變換。其中。 由DFT正變換(1-5-10)可得對(duì)應(yīng)的DFT逆變換為 (1-5-11)DFT正變換公式(1-5-10)和逆變換(1-5-11)給出了序列與序列之間的線性互逆關(guān)系，可記為1.6二維Fourier變換定義1.6.1對(duì)于兩個(gè)實(shí)變量的復(fù)值函數(shù)，稱含兩個(gè)參變量的積分 (1

15、-6-1)為的二維Fourier變換2DFT（2 dimensional Fourier transform或稱為二維譜函數(shù)。對(duì)應(yīng)地稱 (1-6-2)為(1-6-1)的二維Fourier變換的逆變換，又稱為反演公式。1.7二維離散Fourier變換定義1.7.1對(duì)二維復(fù)序列，稱 (1-7-1)為二維離散Fourier變換2DDFT（2dimensional discrete Fourier transform）。對(duì)應(yīng)地稱 (1-7-2)為二維離散Fourier變換（1-7-1）的逆變換。這里。2 DFT快速計(jì)算(FFT)算法分析及其程序?qū)嵗?.1 DFT的快速算法本文將長(zhǎng)度為N的序列變?yōu)殚L(zhǎng)度為

16、N的序列的DFT及其逆變換形式取為(2-1-1)(2-1-2)為了可以用統(tǒng)一的程序處理DFT及其逆變換，我們將(2-1-2)改寫為 (2-1-3)兩個(gè)只差一個(gè)常數(shù)因子，求和可以使用相同的程序。Danielson和Lanczos兩人提出一種將長(zhǎng)度為N的序列的DFT分為兩個(gè)長(zhǎng)度為的序列的DFT。其中一個(gè)序列由原序列的偶數(shù)項(xiàng)構(gòu)成，記其DFT為，而另一個(gè)序列由的奇數(shù)項(xiàng)構(gòu)成，記其DFT為。這時(shí)原序列的DFT為 (2-1-4) (2-1-5)如此將長(zhǎng)度為N的序列的DFT用它的偶序列和奇序列的DFT，即和表示。由于，長(zhǎng)度為可以分別計(jì)算。而長(zhǎng)度為N的序列可按(2-1-4),(2-1-5)求得。只要我們預(yù)先取，

17、當(dāng)我們將長(zhǎng)度為N的序列的DFT轉(zhuǎn)變?yōu)榍髢蓚€(gè)長(zhǎng)為的序列，的DFT后，我們可再對(duì)，分別作同樣的轉(zhuǎn)變，即轉(zhuǎn)變?yōu)榍箝L(zhǎng)度的序列的DFT，這樣該算法便可以遞歸地使用。2.2 基2FFT算法分析現(xiàn)在我們討論，為整數(shù)的FFT算法。將寫成位2進(jìn)制數(shù)，即 (2-2-1)其中只取0或1，。記 (2-2-2) (2-2-3)由的DFT定義得 (2-2-4)其中 (2-2-5)顯然的第一項(xiàng)的因子為1，只有第二項(xiàng)需要作一次復(fù)數(shù)乘法。所以計(jì)算需要次乘法。而的序號(hào)的左第一位為，將左移m-1位和的指數(shù)相符。繼續(xù)對(duì)分解得(2-2-6)其中 (2-2-7)計(jì)算需要次復(fù)數(shù)乘法，同時(shí)的序號(hào)的左兩位為，將其逆排后左移m-2位和的指數(shù)相符

18、。繼續(xù)分解則有 。若記 (2-2-8)則有 (2-2-9)對(duì)某一確定的，計(jì)算需要次復(fù)數(shù)乘法，同時(shí)的序號(hào)的左位為，將其逆排后左移位和的指數(shù)相符。當(dāng)完成步計(jì)算后得 (2-2-10)可以看到，的序號(hào)逆排后得恰好為的序號(hào)。計(jì)算每一需要次復(fù)數(shù)乘法，所以由求得共需要次復(fù)數(shù)運(yùn)算。即當(dāng)時(shí)，求長(zhǎng)度為N復(fù)序列的DFT基2FFT算法運(yùn)算總次數(shù)為 (2-2-11)從以上分析可得：基2FFT算法實(shí)現(xiàn)的關(guān)鍵在于下標(biāo)的組織；基2FFT算法的必要條件為，m為整數(shù)。對(duì)不滿足條件的序列，我們可用零擴(kuò)充的方法，即取最小的，滿足，則可令： (2-2-12)當(dāng)回到時(shí)間域后，我們可以將擴(kuò)充的零元素刪除。2.3 復(fù)序列基2FFT算法設(shè)是長(zhǎng)

19、度為N的復(fù)序列，其中。算法2.3.1基2FFT正變換 記：(1) 計(jì)算；(2) 對(duì)計(jì)算：a.記；b.逆排得；c.左移位得，由此查得；d.對(duì)計(jì)算(3) 對(duì)完成：a. 逆排得b. 令正變換完成，由得其DFT。算法2.3.2基2FFT逆變換置：(1) 計(jì)算(2) 置(3) 用算法2.3.1求的DFT，并記為(4) 對(duì)，令計(jì)算結(jié)束得的DFT逆變換。算法2.3.1和算法2.3.2的程序（fft1.c）在附錄中給出。fft1.c的主要函數(shù)fft1(double A2N,int ifft)程序流程圖如下：結(jié)束開始ifirst=0?Y調(diào)用初始化函數(shù)initial(),并置ifirst=1ifft= -1?N令

20、A1i= -A1iN用算法2.3.1求的DFTYifft= -1?NY置A0i=x0i，A1i=x1i置A0i=A0i/N，A1i=-A1i/N圖2-1 fft1()程序流程圖下面我們給出一個(gè)實(shí)例。例2.3.1 設(shè)，為求其Fourier變換，將其離散化得序列，再求DFT，最后作DFT逆變換。 解：取，這時(shí)取，對(duì)離散化得記得主函數(shù)fft1m1.c（附錄中給出）主要功能：子函數(shù)functf()（附件中給出）由產(chǎn)生，由求其DFT，再由求其DFT逆變換。并將結(jié)果，存入數(shù)據(jù)文件fft1m1.d。fft1m1.c中主函數(shù)的流程圖如圖2-2所示。2.4 二維復(fù)序列2DDFT的行列算法2DDFT的快速算法很多

21、，最常用的是行列算法。對(duì)二維復(fù)序列定義中間序列 (2-4-1)這是個(gè)關(guān)于長(zhǎng)度為的復(fù)序列的一維DFT。進(jìn)而得2DDFT正變換為 (2-4-2)這是個(gè)關(guān)于長(zhǎng)度為的復(fù)序列的一維DFT。這樣，我們就把計(jì)算開始調(diào)用functf(A)， 由f(x)得到離散的序列把作為參數(shù)調(diào)用fft1 (A,1)函數(shù)，由Fourier正變換得到。把作為參數(shù)調(diào)用fft1 (A,1)函數(shù)，由Fourier逆變換得到。結(jié)束輸出輸出輸出圖2-2 fft1m1.c主函數(shù)的流程圖個(gè)點(diǎn)的2DDFT轉(zhuǎn)化為計(jì)算個(gè)長(zhǎng)度為的序列和個(gè)長(zhǎng)度為的序列的一維DFT。以上計(jì)算2DDFT的方法就稱為行列算法。以上給出的是2DDFT行列算法的正變換，對(duì)應(yīng)二維

22、行列算法2DDFT逆變換可以由下面兩個(gè)過(guò)程來(lái)實(shí)現(xiàn)： (2-4-3) (2-4-4)為了能用統(tǒng)一的形式和程序處理，將(2-4-3)取共軛得 (2-4-5)即先將取共軛，關(guān)于對(duì)序列作一維DFT正變換得 (2-4-6)當(dāng)時(shí)，它可以用FFT正變換相關(guān)程序?qū)崿F(xiàn)。求得后，令(2-4-7)求得。從以上算法分析可知對(duì)和，與一維的的要求一樣，需滿足，當(dāng)和不滿足時(shí)，同樣采用零元素?cái)U(kuò)充方法。這時(shí)，我們就可以使用一維基2FFT的程序，從而實(shí)現(xiàn)2DDFT的快速算法。算法2.4.1二維復(fù)序列2DDFT行列算法正變換由求(1) 對(duì)關(guān)于使用一維FFT相關(guān)算法求長(zhǎng)度為的一維序列的DFT: (2-4-8)(2) 對(duì)關(guān)于使用一維F

23、FT相關(guān)算法求長(zhǎng)度為的一維序列的DFT: (2-4-9)算法2.4.2二維復(fù)序列2DDFT行列算法逆變換由求(1) 對(duì)關(guān)于使用一維FFT相關(guān)算法求長(zhǎng)度為的一維序列的DFT逆變換 (2-4-10)(2) 對(duì)關(guān)于使用一維FFT相關(guān)算法求長(zhǎng)度為的一維序列的DFT逆變換 (2-4-11)算法2.4.1和算法2.4.2的程序(fft2.c，rowco2.c)在附錄中給出,fft2.c包括初始化函數(shù)initial()、長(zhǎng)度為的復(fù)序列的FFT程序fft21()、長(zhǎng)度為的復(fù)序列的FFT程序fft22()、逆排序號(hào)的函數(shù)inverseq()。fft21()和fft22()的程序流程與流程圖2.1一樣，在此就不再

24、描述，現(xiàn)在只給出rowco2.c的函數(shù)rowcolumn()的程序流程圖，如圖2-3：開始調(diào)用初始化函數(shù)initial()置j=0置B20j=A0ij,B21j=A1ijj=0,1,N2-1調(diào)用函數(shù)fft22(B2,ifft)置A0ij =B20j,A1ij= B21jj=0,1,N2-1i<N1?i+Y置i=0置B10i=A0ij,B11i=A1iji=0,1,N1-1調(diào)用函數(shù)fft21(B1,ifft)置A0ij =B10i,A1ij= B11ii=0,1,N1-1j<N2?j+YN結(jié)束圖2-3 rowcolumn()的程序流程圖下面我們給出一個(gè)實(shí)例。例2.4.1 設(shè)作離散化，

25、取N1512，N2512，。用行列算法求二維Fourier變換，然后再求其逆變換。解記。將二維Fourier變換離散化為求二維序列的DFT。主函數(shù)程序（fft2m.c）在附錄中給出，所得的數(shù)據(jù)存在文件fft2m.d中。fft1m1.c中主函數(shù)的流程圖如圖2-4所示。結(jié)束開始調(diào)用functf(A)， 由f(x,y)得到離散的序列把作為參數(shù)調(diào)用rowcolumn (A,1)，由Fourier正變換得到。把作為參數(shù)調(diào)用rowcolumn (A,-1)，由Fourier逆變換得到。輸出輸出輸出圖2-4 fft1m1.c中主函數(shù)的流程圖結(jié)論本文主要分析了FFT算法的基本思想，在對(duì)算法分析的基礎(chǔ)上給出其相

26、應(yīng)的程序，最后給出一個(gè)實(shí)例。但由于篇幅有限，本文只給出一維及二維FFT算法及其程序，且對(duì)二維FFT算法只給出了行列算法，對(duì)更高維的DFT變換及其它FFT算法有待進(jìn)一步學(xué)習(xí)研究。致謝首先，非常感謝梁宗旗老師在論文工作中給予的幫助、悉心指導(dǎo)和嚴(yán)格要求。同時(shí)也感謝所有老師，在四年中，他們嚴(yán)謹(jǐn)?shù)闹螌W(xué)態(tài)度和生活中細(xì)致的關(guān)懷使我在學(xué)習(xí)、生活上受益非淺，同時(shí)將激勵(lì)著我在今后的工作中奮進(jìn)。 此外，本論文的順利完成，得到了許建生同學(xué)的幫助，提出很多寶貴的意見(jiàn)，在此表示感謝。 最后，再次向?qū)ξ以谛W(xué)習(xí)期間，所有幫助、支持和鼓勵(lì)過(guò)我的老師，同學(xué)以及我的家人、朋友們呈上我最真摯的謝意！參考文獻(xiàn)1 蔣長(zhǎng)錦等快速傅里葉變

27、換及其C程序合肥：中國(guó)科學(xué)技術(shù)大學(xué)出版社，20042李榮華等微分方程數(shù)值解法北京：高等教育出版社，19973譚浩強(qiáng)C程序設(shè)計(jì)北京：清華大學(xué)出版社，19994翟瑞彩，謝偉松數(shù)值分析天津：天津大學(xué)出版社，20035李開泰，黃艾香等.有限元方法及其應(yīng)用.西安：西安交通大學(xué)出版社，1992 6 徐長(zhǎng)發(fā)，李紅. 實(shí)用偏微分方程數(shù)值解法. 武漢:華中科技大學(xué)出版社，2000.9 7 徐萃薇計(jì)算方法北京：高等教育出版，1985.18 陳傳璋，金福臨，朱學(xué)炎等.數(shù)學(xué)分析（下冊(cè)）.北京：高等教育出版社，1983.119E.O. Brigham. The fast Fourier transform. Prent

28、iceHall, Englewood cliffs, N. J. ,1974.10 H. J. Nussbaumer. Fast Fourier transform and convolution algorithms. Springer-verlag BerlinHeidelbeg New Yorkm, 1981. 附錄附錄1：fft1.c程序文件/*/#include <stdio.h>#include <stdlib.h>#include <math.h>#ifndef DN#define DN 4#endif#define eps 1.0e-9voi

29、d initial(void);void fft1(double A2N,int ifft);unsigned long inverseq(unsigned long int i1,unsigned long int imax);void fft1r2(double A2N,double B2N,int ifft);void fft1r1(double A2DN,int ifft);static int ifirst=0,sign=-1;static unsigned M,jup16,kup16;static double pi,wrN,wiN;/*initial:*/void initial

30、(void)unsigned long int i,j,ni;double unit,tt;ni=N;M=0;while(ni>1)ni=ni/2;M=M+1;ni=1;for(j=0;j<M;j+) ni=ni*2;if(ni-N)!=0) printf("N isn't compatible!nStrike any key to exit!n"); getchar();exit(1);kup0=N/2;/* kup1=2(M-l-1)*/jup0=1;for(j=1;j<M;j+)kupj=kupj-1/2;jupj=jupj-1*2;pi=4

31、.0*atan(double)1.0);unit=2.0*pi/(double)N);tt=0.0;for(i=0;i<N;i+) wri=cos(tt);wii=sign*sin(tt);tt=tt+unit;/* ifft=1<->FFT* ifft=-1<->IFFT*/ void fft1(double A2N,int ifft)unsigned long int l,k,i,j,jl0,jl0k,jl1k,reseq,ij;double wrc,wis ,ar,ai,x2N;if(ifirst=0)initial(),ifirst=1;if(ifft=-

32、1)for(i=0;i<N;i+) A1i=-A1i;for(l=0;l<M;l+) for(j=0;j<jupl;j+)reseq=inverseq(j,l);ij=kupl*reseq;jl0=2*kupl*j;wrc=wrij;wis=wiij;for(k=0;k<kupl;k+)jl0k=jl0+k;jl1k=jl0+kupl+k;ar=A0jl1k*wrc-A1jl1k*wis;ai=A0jl1k*wis+A1jl1k*wrc;A0jl1k=A0jl0k-ar;A1jl1k=A1jl0k-ai;A0jl0k=A0jl0k+ar;A1jl0k=A1jl0k+ai

33、;for(i=0;i<N;i+) j=inverseq(i,M);x0i=A0j;x1i=A1j;for(i=0;i<N;i+)A0i=x0i;A1i=x1i;if(ifft=-1)for(i=0;i<N;i+)A0i=A0i/(double)N);A1i=-A1i/(double)N);/*/unsigned long int inverseq(unsigned long int i1,unsigned long int imax) unsigned long int j,ii,i0,inv;i0=i1;inv=0;for(j=0;j<imax;j+)ii=i0%2;

34、inv=2*inv+ii;i0=(i0-ii)/2;return(inv);/*/附錄2：實(shí)例2.3.1的主函數(shù)程序fft1m1.c#define N 1024#include "fft1.c"void functf(double A2N);void main(void)unsigned long int i;static double A2N;FILE *fp;fp=fopen("ftt1m1.d","w");system("cls");functf(A);fprintf(fp,"The origina

35、l data,Ak:n");printf("The original data,Ak;n");for(i=0;i<N;i+)/fprintf(fp,"%4lu:%15.9f %16.9e in",i,A0i,A1i);fprintf(fp,"%15.9f",A0i);fprintf(fp,"n");fft1(A,1);fprintf(fp,"FFT,direct transform(Ak->xj),xj:n");printf("FFT,direct transfo

36、rm(Ak->xj),xj:n");for(i=0;i<N;i+)/fprintf(fp,"%4lu:%15.9f %16.9e in",i,A0i,A1i);fprintf(fp,"%15.9f",A0i);fprintf(fp,"n");fft1(A,-1);fprintf(fp,"iFFT,direct transform(x->A),Ak:n");printf("FFT,direct transform(x->A),Ak:n");for(i=0;i<

37、;N;i+)/fprintf(fp,"%4lu:%15.9f %16.9e in",i,A0i,A1i);fprintf(fp,"%15.9f",A0i);fprintf(fp,"n");fclose(fp); printf("OK!n");/*/void functf(double A2N)unsigned long int i;double dt,tt,A0,A1;dt=0.1;tt=(double)N*dt;A0=exp(-tt);for(i=0;i<N;i+)tt=(double)i*dt;A1=ex

38、p(-tt);A0i=(A1+A0/A1)*dt;A1i=0.0;/*/附錄3：程序文件fft2.c#include<stdio.h>#include<stdlib.h>#include<math.h>#define eps 1.0e-10unsigned int inverseq(unsigned int i1,unsigned int imax);void fft21(double A2N1,int ifft);void fft22(double A2N2,int ifft);void initial(void);static int ifirst=0;

39、sign=-1;static unsigned M1,M2,ijup16,ikup16,jjup16,jkup16;static double pi,wriN1,wiiN1,wrjN2,wijN2;/*initial: */void initial(void)unsigned int i,j,ni;double unit,tt;ni=N1;M1=0;while(ni>1) ni=ni/2;M1=M1+1;ni=1;for(j=0;j<M1;j+) ni=ni*2;if(ni-N1)!=0) printf("N1 isn't compatible!nStrike a

40、ny key to exit!n"); getchar();exit(1);ni=N2;M2=0;while(ni>1)ni=ni/2;M2=M2+1;ni=1;for(j=0;j<M2;j+) ni=ni*2;if(ni-N2)!=0)printf("N2 isn't compatible!n Strike any key to exit!n");getchar();exit(1);ikup0=N1/2;/* kupl=2(M1-l-1)*/ijup0=1;for(j=1;j<M1;j+)ikupj=ikupj-1/2;ijupj=ij

41、upj-1*2;pi=4.0*atan(double)1.0);unit=2.0*pi/(double)N1);tt=0.0;for(i=0;i<N1;i+) wrii=cos(tt);wiii=sin(tt);tt=tt+unit;jkup0=N2/2;/* jkupl=2(M2-l-1)*/jjup0=1;for(j=1;j<M2;j+)jkupj=jkupj-1/2;jjupj=jjupj-1*2;unit=2.0*pi/(double)N2);tt=0.0;for(i=0;i<N2;i+) wrji=cos(tt);wiji=sin(tt);tt=tt+unit;/*

42、row:ifft=1<->FFT,ifft=-1<->IFFT*/*/void fft21(double A2N1,int ifft)unsigned int l,k,i,j,jl0,jl0k,jl1k,reseq,ij;double wr,wi,ar,ai,x2N1;if(ifft=-1) for(i=0;i<N1;i+) A1i=-A1i;if(ifirst=0)initial();ifirst=1;for(l=0;l<M1;l+) for(j=0;j<ijupl;j+)reseq=inverseq(j,l);ij=ikupl*reseq;jl0=

43、2*ikupl*j;wr=wriij;wi=sign*wiiij;for(k=0;k<ikupl;k+) jl0k=jl0+k;jl1k=jl0+ikupl+k;ar=A0jl1k*wr-A1jl1k*wi;ai=A0jl1k*wi+A1jl1k*wr;A0jl1k=A0jl0k-ar;A1jl1k=A1jl0k-ai;A0jl0k=A0jl0k+ar;A1jl0k=A1jl0k+ai; for(i=0;i<N1;i+) j=inverseq(i,M1);x0i=A0j;x1i=A1j;for(i=0;i<N1;i+)A0i=x0i;A1i=x1i;if(ifft=-1)fo

44、r(i=0;i<N1;i+)A0i=A0i/(double)N1);A1i=-A1i/(double)N1);/*column:ifft=1<->FFT,ifft=-1<->IFFT*/*/void fft22(double A2N2,int ifft)unsigned int l,k,i,j,jl0,jl0k,jl1k,reseq,ij;double wr,wi,ar,ai,x2N2;if(ifirst=0)initial();ifirst=1;if(ifft=-1) for(i=0;i<N2;i+) A1i=-A1i;for(l=0;l<M2;l+) for(j=0;j<jjupl;j+)reseq=inverseq(j,l);ij=jkupl*reseq;jl0=2*jkupl*j;wr=wrjij;wi=sign*wijij;for(k=0;k<jkupl