傅里葉和互相關(guān)算法及c程序解析

1、1.快速傅里葉變換（FFT）1.1葉變換簡(jiǎn)介快速傅里有限長(zhǎng)序列可以通過(guò)離散傅里葉變換（DFT）將其頻域也離散化成有限長(zhǎng)序列. 但其計(jì)算量太大，很難實(shí)時(shí)地處理問(wèn)題，因此引出了快速傅里葉變換（FFT）. 1965 年，Cooley和Tukey提出了計(jì)算離散傅里葉變換（DFT）的快速算法，將DFT的 運(yùn)算量減少了兒個(gè)數(shù)量級(jí)。從此，對(duì)快速傅里葉變換（FFT）算法的研究便不斷 深入，數(shù)字信號(hào)處理這門(mén)新興學(xué)科也隨FFT的出現(xiàn)和發(fā)展而迅速發(fā)展。根據(jù)對(duì)序 列分解與選取方法的不同而產(chǎn)生了FFT的多種算法，基本算法是基-2DIT和基 -2DIFo FFT在離散傅里葉反變換、線性卷積和線性相關(guān)等方面也有重要應(yīng)用。

2、快速傅氏變換（FFT）,是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的 奇、偶、虛、實(shí)等特性，對(duì)離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對(duì)傅氏 變換的理論并沒(méi)有新的發(fā)現(xiàn)，但是對(duì)于在訃算機(jī)系統(tǒng)或者說(shuō)數(shù)字系統(tǒng)中應(yīng)用離散 傅立葉變換，可以說(shuō)是進(jìn)了一大步?？焖俑盗⑷~變換作為一種數(shù)學(xué)方法，已經(jīng)廣泛地應(yīng)用在兒乎所有領(lǐng)域的頻譜 分析中，而且經(jīng)久不衰，因?yàn)樾盘?hào)處理方法沒(méi)有先進(jìn)和落后之分，只有經(jīng)典和現(xiàn) 代之別，在實(shí)際系統(tǒng)中用得最好的方法就是管用的方法。換句話說(shuō)，信號(hào)處理方 法與應(yīng)用背景和口的的貼近程度是衡量信號(hào)處理方法優(yōu)劣的唯一標(biāo)準(zhǔn)。FFT是快 速傅利葉變換（Fast FourierTransform簡(jiǎn)稱(chēng)FFT

3、）的英文縮寫(xiě)，它在為今科技世 界中的應(yīng)用相當(dāng)活躍，無(wú)論是在時(shí)間序列分析領(lǐng)域中，還是在我國(guó)剛剛興起的生 物頻譜治療的研究與應(yīng)用中，都有著重要的作用。同時(shí)，它乂是軟件實(shí)現(xiàn)數(shù)字濾波 器的必備組成部分之一。FFT算法的基本思想：利用DFT系數(shù)的特性，合并DFT運(yùn)算中的某些項(xiàng)，把 長(zhǎng)序列的DFT短序列的DFT,從而減少其運(yùn)算量。FFT算法分類(lèi)：時(shí)間抽選法DIT： Decimation-In-Time :頻率抽選法DIF：Dec imat i onTn - Frequency。 在此以按時(shí)間抽選的基-2FFT算法為例進(jìn)行說(shuō)明。1. 2算法原理N為2的整數(shù)幕的FFT算法稱(chēng)基-2FFT算法。設(shè)N點(diǎn)序列勁），n

4、 =0.1.2,-n-i按n的奇偶分成兩個(gè)長(zhǎng)為N/2的序列r =0,1,2,，x,(r) = x(2r) x2(r) = x(2r + l)由于則兀G）的DFT：其中:N7JV-1JV-1X 伙)=x(n)比丫“為偶數(shù) n為奇數(shù)N N-*1T1= x(20W嚴(yán) + j(2 尸+ 1)W嚴(yán)”r-0r-4)亡后)臨+吧$2(訓(xùn)蟲(chóng)r-0= Xl(/c) + WX2 伙)/U)=0,12，N-1(1) TOC o 1-5 h z X&)=乞%0)脇 2 =疇r-0r-07_,7_,X?伙)=2+1)晞2 =J2(門(mén)晞2r-()r-0按式（l）計(jì)算得到的只是X的前一半項(xiàng)數(shù)的結(jié)果，即（OWkW?1）。由

5、系數(shù)的周期性可推岀X（燈的后一半值，即比2=一恍乂因?yàn)?伙），X2（k）是以什為周期的，可推出：4AX + N/2） = 5）叫：嚴(yán)” =5“）w為=X）1)(4)同理可得,X2(k + N/2)=X)(5)r-0r-0由此可得N點(diǎn)對(duì)應(yīng)的DFTX伙）的計(jì)算式子:X伙)7+比億,=012.X 伙 + NI2) = X(k)W：X(k)上一式計(jì)算X（燈的前一半值，下一式計(jì)算X（k）的后一半值。因此只要求岀01區(qū)間內(nèi)各個(gè)整數(shù)k值所對(duì)應(yīng)的X和X2（k）,便可求出0到N1點(diǎn)全部的X值。明顯節(jié)約了運(yùn)算量。式中因子W；在復(fù)數(shù)乘法中起一個(gè)旋轉(zhuǎn)的作用，稱(chēng)為旋 轉(zhuǎn)因子。公式（6）的運(yùn)算可以歸結(jié)為兩個(gè)復(fù)數(shù)&、b求得

6、復(fù)數(shù)n+MVj和d叫。 用信號(hào)流程圖的方法可以簡(jiǎn)單的表示為如圖2所示。這樣的運(yùn)算稱(chēng)為蝶形運(yùn)算, 在FFT算法中占有核心的地位。圖2時(shí)間抽選蝶形運(yùn)算流圖顯然，每個(gè)蝶形運(yùn)算對(duì)應(yīng)于一次復(fù)數(shù)乘法和兩次復(fù)數(shù)加法運(yùn)算。如果用DFT 方法直接計(jì)算出X”）和X”），共需與次復(fù)數(shù)乘法運(yùn)算，再作N/2次蝶形運(yùn)算，乂 需N/2次復(fù)數(shù)乘法和N次復(fù)數(shù)加法。這樣算出N點(diǎn)DFT共需要（T+n） /2次復(fù)數(shù) 乘法和（N2/2）+N次復(fù)數(shù)加法。當(dāng)N較大時(shí)，同直接計(jì)算N點(diǎn)DFT所需的次復(fù) 數(shù)乘加次數(shù)相比，兒乎減少了一半的計(jì)算量。假設(shè)N/2還是偶數(shù)，則X）和兀這兩個(gè)N/2點(diǎn)的DTF計(jì)算，乂分別可以通 過(guò)計(jì)算N/4點(diǎn)DFT和蝶形運(yùn)算

7、而得到。這時(shí)蝶形有兩組，每組N/4個(gè)，總數(shù)也是 N/2個(gè)，所以也需要N/2次復(fù)數(shù)乘法和N次復(fù)數(shù)加法。如果N=2W ,則分解過(guò)程一 直可以進(jìn)行下去，共分解M次，到2點(diǎn)DFT為止。以上所述就是FFT算法的核心思想?？梢钥闯觯@種罪基本形式的快速傅里 葉算法要求點(diǎn)數(shù)是2的幕次。當(dāng)序列長(zhǎng)度不具有尸的形式時(shí)，可以補(bǔ)上一段0,使 總長(zhǎng)度為2。13時(shí)間抽取過(guò)程的流圖表示現(xiàn)在就N=8的悄況對(duì)算法結(jié)構(gòu)，即計(jì)算流圖作詳細(xì)的說(shuō)明，山此可以直接得 到關(guān)于N=2”的情況下得一些一般性結(jié)論。圖3是8點(diǎn)DFT歐諾個(gè)兩個(gè)4點(diǎn)DFT和4 個(gè)蝶形運(yùn)算實(shí)現(xiàn)的圖示，圖4中進(jìn)一步將4點(diǎn)DFT作類(lèi)似的分解，最后將圖4中的2 點(diǎn)DFT也畫(huà)

8、成蝶形運(yùn)算流圖，如圖5所示。整個(gè)運(yùn)算有3輪蝶形運(yùn)算構(gòu)成，每一輪 有4個(gè)蝶形，但它們的旋轉(zhuǎn)因子和循環(huán)方式各輪有所不同。X-1A(0) 丄 雙4) x(6)x 班3) x(5) i(7)x(0)(4)1兀鞏6)14)1x(5) -x(3) -兀e-X x(7 x(X XXX(l)X(2)XXX(5)X(6)圖4 8點(diǎn)DFT時(shí)間抽選分解（續(xù)）X(0)X(l)X(2)XXX(5)X(6)X圖5反序輸入順序輸出時(shí)間抽選FFT算法流圖圖5那樣的結(jié)構(gòu)對(duì)任何N=2的序列都是適用的。整個(gè)N點(diǎn)DFT的訃算可用M=log2JV輪蝶形運(yùn)算構(gòu)成，每輪有N/2個(gè)蝶形，因此總的運(yùn)算就由這log川個(gè)蝶形 2所構(gòu)成。因此總的復(fù)

9、數(shù)乘法次數(shù)是|Nlog2A復(fù)數(shù)加法次數(shù)是Nlog2 N o21.4運(yùn)算量下面討論按時(shí)間抽取法的運(yùn)算量，對(duì)于任意一個(gè)N=2“的DFT運(yùn)算，都可以通 過(guò)M次分解，最后分解成2點(diǎn)DFT運(yùn)算的組合。這樣的M次分解，就構(gòu)成山x(”)到X(R) 的M級(jí)運(yùn)算過(guò)程，每一級(jí)運(yùn)算都有N/2個(gè)蝶式運(yùn)算完成。而每一蝶式運(yùn)算有一次 復(fù)數(shù)乘法和兩次復(fù)數(shù)加法。因此，每一級(jí)運(yùn)算都需要N/2次復(fù)數(shù)乘法和N次復(fù)數(shù) 加法。這樣M級(jí)運(yùn)算共需要復(fù)數(shù)乘法：叫=(N / 2)M = (N / 2) log, N復(fù)數(shù)加法： ap=N-M=N log2 N為了比較FFT算法與直接DFT運(yùn)算的運(yùn)算量，現(xiàn)將二者在不同N值時(shí)的乘法 次數(shù)列于表1中，

10、并將直接計(jì)算DFT與FFT算法的計(jì)算量之比N2 / (N / 2) log? N也 列于表中。從表中可以看出，當(dāng)N較大時(shí)，時(shí)間抽取法要比直接計(jì)算塊一、二個(gè)數(shù)量級(jí)。表1 FFT算法與直接計(jì)算DFT算法的比較NN2(AT/2)log2 AT/V2/(A/2)log2/V2斗14.0416斗4.0864125.41625632&03210248012.864409619221.41281638444836.6256262144102464.051210485762304113.81024419430411264204.815按時(shí)間抽取的FFT算法特點(diǎn)（1）原位運(yùn)算（同址運(yùn)算）山所述算法原理及圖5的N

11、2的信號(hào)流圖，F(xiàn)FT的每級(jí)（列）訂算都是山N個(gè) 復(fù)數(shù)據(jù)經(jīng)N/2個(gè)蝶形運(yùn)算變成了另外N個(gè)復(fù)數(shù)據(jù)。每一個(gè)蝶形運(yùn)算結(jié)構(gòu)完成下述 基本迭代運(yùn)算4.（0 = 4,-盅皿心-盅_心）嗆式中加表示第加列迭代，（、_/為數(shù)據(jù)所在的行數(shù)。上式的運(yùn)算如圖6所示，由一次復(fù)數(shù)乘法和一次復(fù)數(shù)加減法組成。任何兩個(gè)節(jié)點(diǎn)i和丿的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算后，其結(jié)果為下一列的i和丿兩點(diǎn)的節(jié)點(diǎn)變量，而和其他節(jié)點(diǎn)變量無(wú)關(guān)。4”（刀=九/）一心心）叭圖6時(shí)間抽取蝶形運(yùn)算結(jié)構(gòu)（2）輸入序列的序號(hào)及整序規(guī)律由圖5可見(jiàn),當(dāng)按原位進(jìn)行計(jì)算時(shí),FFT輸出端的X（的次序正好是順序排列的。但這時(shí)輸入沏）卻不能按自然順序存入存儲(chǔ)單元中，而是按x（0）,

12、x（4）, x（2）, a（6）,的順序存入存儲(chǔ)單元，因而是亂序的。這就使得運(yùn)算時(shí)取數(shù)據(jù)的地址編 排“混亂無(wú)序”。但實(shí)際上是山規(guī)律可尋的，我們稱(chēng)為倒位序。亂序的原因是山按時(shí)間抽取進(jìn)行的FFT運(yùn)算的原理造成的，是山于輸入測(cè)按 標(biāo)號(hào)n的奇偶不斷分組造成的。如果你用二進(jìn)制數(shù)表示為（”2耳心）2 （當(dāng)N=8=23時(shí), 用三位二進(jìn)制符號(hào)表示），第一次分組n為偶數(shù)，在上半部分，n為奇數(shù)在下半部 分，這可以觀察n的二進(jìn)制的最低位心，弘嚴(yán)0則序列值對(duì)應(yīng)于偶數(shù)抽樣，必=1 則序列值對(duì)應(yīng)于奇數(shù)抽樣。下一次則根據(jù)次最低位q的0, 1來(lái)分奇偶（而不管原 來(lái)的序列的子序列是偶序列還是奇序列）。在實(shí)際運(yùn)算中，直接將輸入數(shù)

13、據(jù).g）按原位運(yùn)算要求的“亂序”存放是很不 方便的。因此總是先按自然順序?qū)⑤斎胄蛄写嫒氪鎯?chǔ)單元，再通過(guò)變址運(yùn)算將自 然順序變換成按時(shí)間抽取的FFT算法要求的順序。變址的過(guò)程可以用程序安排加 以實(shí)現(xiàn)，稱(chēng)為“整序”或“重排”。整序規(guī)律：如果輸入序列的序號(hào)用二進(jìn)制數(shù)表示，如若其反序二進(jìn)制 気用來(lái)表示，則原來(lái)在自然順序時(shí)應(yīng)該放測(cè)的存儲(chǔ)單元，現(xiàn)在放著的是応）。表 2列出了 N=8時(shí)的自然順序二進(jìn)制數(shù)以及相應(yīng)的倒位序二進(jìn)制數(shù)。表2順序和倒位序二進(jìn)制數(shù)自然順序n二進(jìn)制數(shù)倒位序二進(jìn)制數(shù)倒位序順序J）00000000100110042010010230111106410000115101101561100113

14、71111117將按自然順序存放在存儲(chǔ)單元中的數(shù)據(jù)，調(diào)換成FFT原位運(yùn)算所要求的倒位序的變址功能，如圖8所示。當(dāng)n二“時(shí)，數(shù)據(jù)不需要調(diào)換，當(dāng)nH”時(shí)，必須將原 來(lái)存放數(shù)據(jù),W的存單元內(nèi)調(diào)入數(shù)據(jù)山），而將存放加）的單元內(nèi)調(diào)入勸）。為了避免再次考慮前面已調(diào)換過(guò)得數(shù)據(jù)，保證調(diào)換只進(jìn)行一次，只需檢查一下：是否比n小，假若,;比n小，則意味著在前面已經(jīng)與互相調(diào)換過(guò)。只有當(dāng),;比口大時(shí),才將原存放射）及存放的存儲(chǔ)單元的內(nèi)容互相調(diào)換。經(jīng)過(guò)上述變址運(yùn)算以后, 所得到的數(shù)據(jù)順序正是輸入所要求的順序。與圖4的輸入順序?qū)φ眨浑y看出上 述整序規(guī)律是正確的。存儲(chǔ)單元期0)A(l)AAAA(5)A（6）A(7)自然順序

15、輸入x(0)x(l)x(2)xx(4)x(5)1 一x(6)x(7)變址1二*亡1倒位序x(0)x(4)x(2)x(6)X(l)x(5)rx(7)圖8倒位序的變址處理（3）各類(lèi)蝶形運(yùn)算兩節(jié)點(diǎn)的“距離”及嗆的變化規(guī)律以8點(diǎn)FFT為例。第一列蝶形運(yùn)算只有一種類(lèi)型：系數(shù)為V=l,參加蝶形運(yùn) 算的兩個(gè)數(shù)據(jù)點(diǎn)的間距為1。笫二列有兩種類(lèi)型的蝶形運(yùn)算：系數(shù)分別為叨，w” 參加蝶形運(yùn)算的兩個(gè)數(shù)據(jù)節(jié)點(diǎn)間的間距為2。第三列有4種類(lèi)型的蝶形運(yùn)算：系數(shù) 分別為。，眈，吧，町，參加蝶形運(yùn)算的兩個(gè)數(shù)據(jù)節(jié)點(diǎn)間的間距為4?？梢?jiàn)， 每列的蝶形運(yùn)算類(lèi)型比前1列增加1倍，參加蝶形運(yùn)算的兩個(gè)數(shù)據(jù)節(jié)點(diǎn)間的間距也 增大1倍。由此尅類(lèi)推，

16、對(duì)于N=2w點(diǎn)FFT,當(dāng)輸入為倒位序，輸出為正常順序時(shí), 其第m級(jí)運(yùn)算，參與每個(gè)蝶形運(yùn)算的數(shù)據(jù)節(jié)點(diǎn)間距為2i。針對(duì)具體的第m級(jí)蝶形運(yùn)算，一個(gè)DIT蝶形運(yùn)算的來(lái)那個(gè)節(jié)點(diǎn)“距離”為2心， 因而式7和式8可寫(xiě)成A） = A,M）+ A”/ + 2心）（9）4貝+ 2心）=盅_衛(wèi)）-札/ + 2心）叭(10)叱;中的廠可用兩種方法確定O第1種方法：逆推法。山于蝶形運(yùn)算最后一級(jí)適用了N/2中系數(shù)，分寫(xiě)為WSW ，用尹，而前一列使用了它后面一列所用系數(shù)對(duì)應(yīng)的偶數(shù)序號(hào)的一半，依 次類(lèi)推，可以求出所有的需要適用的系數(shù)。第2種方法：直接求解。首先把式9或式10中兩節(jié)點(diǎn)中的第1個(gè)節(jié)點(diǎn)標(biāo)號(hào) 值，即i值，表示成L位（

17、N=2l）二進(jìn)制數(shù)。然后把此二進(jìn)制數(shù)乘上2,即 將此L位二進(jìn)制數(shù)左移厶-加位（注意m是第m級(jí)的運(yùn)算），把右邊空出的位置 補(bǔ)零，此數(shù)就是所求廠的二進(jìn)制數(shù)。1.6快速傅立葉變換的C語(yǔ)言實(shí)現(xiàn)IIIDIT基-2FFT算法原理及特點(diǎn)可知，F(xiàn)FTi十算程序主要包括變址和M級(jí)遞推計(jì) 算兩部分。（1）變址（倒序）運(yùn)算在實(shí)際運(yùn)算中，一般總是按自然順序?qū)⑤斎胄蛄写嫒氪鎯?chǔ)單元，因此為實(shí)現(xiàn) DIT基-2FFT首先必須將輸入序列巾）按二進(jìn)制倒序數(shù)重排。由表2知，自然順序每 增加1,相應(yīng)于二進(jìn)制順序數(shù)最低位加1,并向左進(jìn)位。而倒序數(shù)則是在M位二進(jìn) 制數(shù)最高位加1,并向右進(jìn)位，用這種反向進(jìn)位加法可以從當(dāng)前任意倒序數(shù)值求 得

18、下一個(gè)倒序數(shù)值，完成倒序重排。反向進(jìn)位加法實(shí)現(xiàn)變址運(yùn)算的訃算流程，如圖9所示，稱(chēng)為雷德（Rader）算 法。設(shè)A（I）表示存放原自然順序輸入數(shù)據(jù)的內(nèi)存單元，A（J）表示存放倒序順序數(shù) 據(jù)的內(nèi)存單元。IJ都從0開(kāi)始。若已知某倒序數(shù)J,要求下一個(gè)倒序數(shù)，先判斷J 的最高位是否為“0”，這可與K二N/2相比較，因?yàn)镹/2的二進(jìn)制描述中，最高位 為1,其余位為0,因?yàn)槿绻鸎 J,貝叮的最高位為“0”，只要該位變?yōu)椤?”（令J 與K二N/2相加即可），就得到了下一個(gè)倒序數(shù)。如果KWJ,則J的最高位為“1”， 可將其改為“0”（減去/2即可）。然后再判斷次高位（與K二N/4比較），若其為“0”, 則將它改

19、為“1”（與N/4相加），此外還需判斷再下一位（與K二N/8比較）。如此依 次進(jìn)行，直到遇到某位為“0”，這時(shí)把這個(gè)“0”改為“1”，就得到了下一個(gè)倒 序數(shù)。求得新的倒序數(shù)J以后，當(dāng)IM?圖10 S-2DIT-FFT算法流程圖2.用FFT計(jì)算相關(guān)函數(shù)相關(guān)概念很重要，互相關(guān)運(yùn)算廣泛應(yīng)用于信號(hào)分析與統(tǒng)計(jì)分析，如通過(guò)相關(guān) 函數(shù)峰值的檢測(cè)測(cè)量?jī)蓚€(gè)信號(hào)的時(shí)延差等。兩個(gè)長(zhǎng)為的實(shí)離散時(shí)間序列x (n)與y (n)的互相關(guān)函數(shù)定義為N-lN-1% CO =工戈0 -廠)丿(刃)=工戈0)歹5 + r)72=0M=0則可以證明，rX7( t )的離散付里葉變換為RS7 (k) =X*(k)Y(k)OWkWN-1

20、其中 X (k)二DFT x (n), Y (k)二DFT y (n) , R=. (k)二DFT g ( t )證：將x(n)、y(n)的逆離散付里葉變換代入互相關(guān)函數(shù)定義式吧牡1匕芒恥1阻的(旳切(刃=工戒捕血+勺=立亦丫總眇艸叱閃沁啊NgN口因x (n)是實(shí)序列，所以x (n) =x* (n),得K-l 1 AT-1-ihi 1 1V-1酥）=藝豈2X（砂5嘰1竺那fe =o TV zjN mi1N-1N-1_仗丄 ” d-SN fc=o f=o1 生 1 J 務(wù)（i） 1 1 0 N亓百 -N 那旳1 0 NNTi NT二工戈0”（刃 F 二Z=0Ng證畢。當(dāng)x (n) =y (n)時(shí)

21、，得到x (n)的自相關(guān)函數(shù)為：N-Li K-lj 蘭乂冷二x(n)x(n + i：) =石習(xí)乂仮)比11it-oN k.o上面的推導(dǎo)表明，互相關(guān)和自相關(guān)函數(shù)的計(jì)算可利用FFT實(shí)現(xiàn)。山于離散付 里葉變換隱含著周期性，所以用FFT計(jì)算離散相關(guān)函數(shù)也是對(duì)周期序列而言的。 直接做點(diǎn)FFT相當(dāng)于對(duì)兩個(gè)點(diǎn)序列x(n)、y(n)作周期延拓，作相關(guān)后再取 主值(類(lèi)似圓周卷積)。而實(shí)際一般要求的是兩個(gè)有限長(zhǎng)序列的線性相關(guān)，為避 免混淆，需采用與圓周卷積求線性卷積相類(lèi)似的方法，先將序列延長(zhǎng)補(bǔ)0后再用 上述方法。利用FFT求兩個(gè)有限長(zhǎng)序列線性相關(guān)的步驟：設(shè)x (n)長(zhǎng)為Nl, y (n)長(zhǎng)為N2,求線性相關(guān)。(1

22、)為了使兩個(gè)有限長(zhǎng)序列的線性相關(guān)可用其圓周相關(guān)代替而不產(chǎn)生混淆，選 擇周期NMNi+M ,N=2,以便使用FFT,將x (n) , y (n)補(bǔ)零至長(zhǎng)為N。即：n = o丄葛-1=pO)”=站，“十 1,.,預(yù)一 1 0(2)用 FFT 計(jì)算 X(k), Y(k) (k=O, 1,NT)R (k)二片(k) Y (k)對(duì) R (k)作 IFFT,得到 r (n) (n二0, b ，N-l)先創(chuàng)建兩個(gè)M格式的文件，分別為和y.txt,用來(lái)存放互相關(guān)變換的兩組數(shù)據(jù)；運(yùn)行 完產(chǎn)生result.txt和result2.txt和result.txt,用來(lái)保存結(jié)果。#includc includeincl

23、ude#includc #includc using namespace std;#dcfine PI 3.1415926535897932384626433832795028841971#dcfine N 128 typedef structdoublereal;/* 實(shí)部*/doubleimg;/* 虛部 *7 complex;complex xN, yN.RN, *W7* 輸岀序列的值*/ intsizc_x=0sizc_y=0cngth=0: 序列長(zhǎng)度bxomplexvoid add(complex a,complexc-real=a.real+b.real;c- img=a.i mg

24、+b .i mg;void sub(complex a,complex plexc-real=a.real-b.real;c-img=a.img-b.img;plexvoid mul(complex axomplexc-real=a.real*b.real a.img*b.img; c-img=a.real*b.img + a.img *b.rcal:void divi(complex plex b,complex *c)c-real=( a.real*b.real+a.img*b.img )/(b.real*b.rcal+b.img*bimg);c-img=( a.img*b.real-a.

25、real*b.img)/(b.real*b.real+b.img*b.img);void initW(int sizc_x)核運(yùn)算int i;W=(complex *)malloc(sizeof(complex) * sizc_x);for(i=0;i0)jink &1)；k=kl;if(ji)temp=xi;xi=xj;xj=temp;void fftx()快速傅里葉函數(shù)int i=OJ=O,k=O,l=O:complex up.dowikproduct;change();for(i=0;ilog(size_x)/log (2);i+) / 一級(jí)蝶形運(yùn)算 */1=1 i;for(j=0;js

26、ize_x;j+=2*l) /* 一組蝶形運(yùn)算 */for(k=0;kl;k+) /個(gè)蝶形運(yùn)算 */mul(x j+k+1 , Wsize_x *k/2/l product);add(xj+k,product,&up);sub(xj+duct.&down);xj+k=up;xj+k+l=down;void ffty()快速傅里葉函數(shù)int i=Oj=O,k=OJ=O;complexup.down. product;change();for(i=0;ilog (size_y)/log (2);i+) I*級(jí)蝶形運(yùn)算 */ 1=1 i;for(j=0;jsize_jf；j+=2*l) /

27、 一組蝶形運(yùn)算 */for(k=0;kl;k+) /個(gè)蝶形運(yùn)算 */mul(yj+k+l.Wsize_x*k/2/l,&product);add(yj+kbproduct.&up);sub(yj+duct.&down);yj+k=up:yj+k+l=down: void ifft() 傅里葉逆變換length=size_x+size_y-l;int i=OJ=O,k=OJ=length;complex up.down;for(i=0;i(int)(log(length)/log(2);i+) /*級(jí)蝶形運(yùn)算*/l/=2;for(j=0;jlength;j+=2*I) /*組蝶形運(yùn)算

28、*/ for(k=0;kl;k+) /*個(gè)蝶形運(yùn)算add(Rj+k.Rj+k+l.&up);up. real/=2:up.i mg/=2;sub(Rj+k,Rj+k+l,&down);down.reaI/=2;down.img/=2;divi(down,WIength*k/2/l,&down);Rj+k=up:Rj+k+l=down;change();void outputone() /*輸出 x 結(jié)果*/int i;printf(iix傅里葉變換結(jié)果n”)；for(i=0;i=0.0001)printf(H+%.4fj ,xi.img);else if(fabs(xi.img)0.000I)

29、printf(H+O.OOOOj ”)；elseprintf(N%.4fj H,xi.img);printf(HnH);void outputtwo()/*輸岀 x 結(jié)果*/int i;printf(ny傅里葉變換結(jié)果n”)；for(i=0;i=0.0001)printf(H+%.4fj ,yi.img);else if(fabs(yi.img)0.0001)printf(H+0.0000j ”)；elseprintf(N%.4fj jfiJ.img);printf(HnH);void outputthree() /*輸出 xy 互相關(guān)結(jié)果*/int i;printf(nxy互相關(guān)變換結(jié)果n”

30、)；for(i=0;i=0.0001)printf(N+%.4fj ：Riimg);else if(fabs(Ri.img)0.0001)printf(H+O.OOOOj N);elseprintf(N%.4fj ：Riimg);printf(Hn,r);void saveone()/保存x傅里葉變換結(jié)果int i;ofstream outfile(nresult l.txfjosout);if(!outfile) cerrHopen resultl.txt erro!Hendk exit( 1);outfileHx傅里葉變換結(jié)果:Mendl;for(i=0;i=0.0001)else if(

31、fabs(xi.img)=0.0001) ouifilcvv”+”vvyiimgvv 了 vv_;else if(fabs(yi.img)0.0001) outfilcvv*Nv”0000(Tvvyvv”H;elseoutfileyi.imgHj,H ”；printf(HnH);outfile.closeO;void savethree()/保存xy互相關(guān)變換結(jié)果int i;ofstream oiHfilc(rcsulUxF；ios:out); if(!outfile)cerrHopen result.txt erro!Mendl;exit(l);)outfileHxy互相關(guān)變換結(jié)果:,end

32、l: for(i=0:i=0.0001) outfile,+,Ri.img,jHH H;else if(fabs(Ri.img)0.0001)outfileH+,M0.0000,l,j,rH M; elseoutfileRi.img,j,t H;primfTE);outfile.close();)int main()/*對(duì)X進(jìn)行傅里葉變換*/int i;size_x=10;ifstream infile(Mx.txtHJos:in);if(!infile)cerrHopen erro!Hendl; exit(l);coutHx 序列的值:,rendl; for(i=0;isize_x;i+)i

33、nfilexi.real;xi.img=0;printf(H%.4f M,xi); initW(size_x);outputoneO;saveone(); infile.closeO;/*對(duì)y進(jìn)行傅里葉變換*/intj;size_y=10;ifstream infiles(Hy.txt,ios:in); if( Jinfiles)cerrHopen erro!Mendl; exit(l);coutHy 序列的值:,fendl;for(J=0：jvsizc_y;j+)infilesy|j.real;yj.img=O:printf(n%.4f ”,yj); initW(size_y); ffty();outputt