傅里葉變換及C語言實現(xiàn)

1、快速傅氏變換，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對 離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機(jī) 系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進(jìn)了一大步。設(shè)x（n）為N項的復(fù)數(shù)序列，由DFT變換，任一 X（ m）的計算都需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù) 加法，而一次復(fù)數(shù)乘法等于四次實數(shù)乘法和兩次實數(shù)加法，一次復(fù)數(shù)加法等于兩次實數(shù)加法，即使 把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次運(yùn)算”（四次實數(shù)乘法和四次實數(shù)加法），那么求出 N項復(fù)數(shù)序列的X（m）,即N點(diǎn)DFT變換大約就需要N2次運(yùn)算。當(dāng)N=1024點(diǎn)甚至更多的時候

2、， 需要N2=1048576次運(yùn)算，在FFT中，利用WN的周期性和對稱性，把一個 N項序列（設(shè)N=2k,k 為正整數(shù)），分為兩個 N/2項的子序列，每個 N/2點(diǎn)DFT變換需要（N/2）2次運(yùn)算，再用N次運(yùn) 算把兩個N/2點(diǎn)的DFT變換組合成一個N點(diǎn)的DFT變換。這樣變換以后，總的運(yùn)算次數(shù)就變成 N+2（ N/2）2二N+N2/2。繼續(xù)上面的例子，N=1024時，總的運(yùn)算次數(shù)就變成了 525312次，節(jié)省 了大約50%的運(yùn)算量。而如果我們將這種一分為二”的思想不斷進(jìn)行下去，直到分成兩兩一組的 DFT 運(yùn)算單元，那么N點(diǎn)的DFT變換就只需要Nlog2N次的運(yùn)算，N在1024點(diǎn)時，運(yùn)算量僅有102

3、40 次，是先前的直接算法的1%,點(diǎn)數(shù)越多，運(yùn)算量的節(jié)約就越大，這就是 FFT的優(yōu)越性.傅里葉變換（Transform ee de Fourier ）是一種積分變換。因其基本思想首先由法國學(xué)者傅里葉系統(tǒng) 地提出，所以以其名字來命名以示紀(jì)念。應(yīng)用傅里葉變換在物理學(xué)、數(shù)論、組合數(shù)學(xué)、信號處理、概率論、統(tǒng)計學(xué)、密碼學(xué)、聲學(xué)、光學(xué)、海洋 學(xué)、結(jié)構(gòu)動力學(xué)等領(lǐng)域都有著廣泛的應(yīng)用（例如在信號處理中，傅里葉變換的典型用途是將信號分 解成幅值分量和頻率分量）。概要介紹傅里葉變換能將滿足一定條件的某個函數(shù)表示成三角函數(shù)（正弦和/或余弦函數(shù)）或者它們的積分的線性組合。在不同的研究領(lǐng)域，傅里葉變換具有多種不同的變體形

4、式，如連續(xù)傅里葉變換和離散傅 里葉變換。最初傅里葉分析是作為熱過程的解析分析的工具被提出的（參見：林家翹、西格爾著自然科學(xué)中確定性問題的應(yīng)用數(shù)學(xué)，科學(xué)出版社，北京。原版書名為C. C. Lin & L. A. Segel,Mathematics Applied to Deterministic Problems in the Natural Sciences, Macmillan Inc., New York, 1974 ）。傅里葉變換屬于諧波分析。傅里葉變換的逆變換容易求出，而且形式與正變換非常類似；正弦基函數(shù)是微分運(yùn)算的本征函數(shù)，從而使得線性微分方程的求解可以轉(zhuǎn)化為常系數(shù)的代數(shù)方程的 求解

5、在線性時不變的物理系統(tǒng)內(nèi)，頻率是個不變的性質(zhì)，從而系統(tǒng)對于復(fù)雜激勵的響應(yīng)可以通過組合 其對不同頻率正弦信號的響應(yīng)來獲??；卷積定理指出：傅里葉變換可以化復(fù)雜的卷積運(yùn)算為簡單的乘積運(yùn)算，從而提供了計算卷積的一種簡 單手段；離散形式的傅里葉變換可以利用數(shù)字計算機(jī)快速的算出（其算法稱為快速傅里葉變換算法（FFT）.快速傅立葉變換FFT的C語言程序2008年 08 月 29 日 星期五 22:59/ 入口參數(shù)：/ l: l = 0, 傅立葉變換 ; l = 1, 逆傅立葉變換/ il: il = 0, 不計算傅立葉變換或逆變換模和幅角； il = 1, 計算模和幅角/ n: 輸入的點(diǎn)數(shù)，為偶數(shù)，一般為3

6、2，64，128，.,1024等/ k: 滿足n=2Ak(k0),實質(zhì)上k是n個采樣數(shù)據(jù)可以分解為偶次幕和奇次幕的次數(shù)/ pr: 1=0 時，存放N點(diǎn)采樣數(shù)據(jù)的實部/ 1=1 時，存放傅立葉變換的N個實部/ pi： 1=0 時，存放N點(diǎn)采樣數(shù)據(jù)的虛部/ 1=1 時，存放傅立葉變換的N個虛部/ / 出口參數(shù)：/ fr: 1=0, 返回傅立葉變換的實部/ 1=1, 返回逆傅立葉變換的實部/ fi: 1=0,返回傅立葉變換的虛部/ 1=1, 返回逆傅立葉變換的虛部/ pr: i1 = 1,1 = 0/ i1 = 1,1 = 1/ pi: i1 = 1,1 = 0/ i1 = 1,1 = 1時，返回傅

7、立葉變換的模時，返回逆傅立葉變換的模時，返回傅立葉變換的輻角時，返回逆傅立葉變換的輻角void kbfft(doub1e *pr,doub1e *pi,int n,int k,doub1e *fr,doub1e *fi,int 1,int i1) int it,m,is,i,j,nv,10;doub1e p,q,s,vr,vi,poddr,poddi;/ 排序for (it=0; it=n-1; it+) m=it; is=0;for (i=0; i=k-1; i+) j=m/2; is=2*is+(m-2*j); m=j;frit=pris; fiit=piis;/ 蝶形運(yùn)算pr0=1.0;

8、 pi0=0.0; p=6.283185306/(1.0*n);pr1=cos(p); pi1=-sin(p);if (1!=0) pi1=-pi1;for (i=2; i=n-1; i+) p=pri-1*pr1; q=pii-1*pi1;s=(pri-1+pii-1)*(pr1+pi1); pri=p-q; pii=s-p-q;for (it=0; it=0; l0-) m=m/2; nv=2*nv;for (it=0; it=(m-1)*nv; it=it+nv)for (j=0; j=(nv/2)-1; j+) p=prm*j*frit+j+nv/2;q=pim*j*fiit+j+nv

9、/2;s=prm*j+pim*j; s=s*(frit+j+nv/2+fiit+j+nv/2); poddr=p-q; poddi=s-p-q; frit+j+nv/2=frit+j-poddr; fiit+j+nv/2=fiit+j-poddi; frit+j=frit+j+poddr; fiit+j=fiit+j+poddi;if (l!=0)for (i=0; i=n-1; i+) fri=fri/(1.0*n);fii=fii/(1.0*n);if (il!=0)for (i=0; i=n-1; i+) pri=sqrt(fri*fri+fii*fii);pri=(pri/(n/2);

10、 / 各次諧波幅值，其中 pr1 為基波幅值if (fabs(fri)0.00001 & f0) pii=90.0;else pii=-90.0;else pii=atan(fii/fri)*360.0/6.283185306;return;快速傅立葉變換FFT的C語言程序2008年 08月 29日 星期五 22:59/ 入口參數(shù)：/ l: l = 0, 傅立葉變換 ; l = 1, 逆傅立葉變換/ il: il = 0, 不計算傅立葉變換或逆變換模和幅角； il = 1, 計算模和幅角/ n: 輸入的點(diǎn)數(shù)，為偶數(shù)，一般為32，64，128，.,1024等/ k: 滿足n=2Ak(k0),實質(zhì)

11、上k是n個采樣數(shù)據(jù)可以分解為偶次幕和奇次幕的次數(shù)/ pr: 1=0 時，存放N點(diǎn)采樣數(shù)據(jù)的實部/ 1=1 時，存放傅立葉變換的N個實部/ pi： 1=0 時，存放N點(diǎn)采樣數(shù)據(jù)的虛部/ 1=1 時，存放傅立葉變換的N個虛部/ / 出口參數(shù)：/ fr: 1=0,返回傅立葉變換的實部/ 1=1, 返回逆傅立葉變換的實部/ fi: 1=0,返回傅立葉變換的虛部/ 1=1, 返回逆傅立葉變換的虛部/ pr: i1 = 1,1 = 0/ i1 = 1,1 = 1/ pi: i1 = 1,1 = 0/ i1 = 1,1 = 1時，返回傅立葉變換的模時，返回逆傅立葉變換的模時，返回傅立葉變換的輻角 時，返回逆

12、傅立葉變換的輻角void kbfft(doub1e *pr,doub1e *pi,int n,int k,doub1e *fr,doub1e *fi,int 1,int i1) int it,m,is,i,j,nv,10;doub1e p,q,s,vr,vi,poddr,poddi;/ 排序for (it=0; it=n-1; it+) m=it; is=0;for (i=0; i=k-1; i+) j=m/2; is=2*is+(m-2*j); m=j;frit=pris; fiit=piis;/ 蝶形運(yùn)算pr0=1.0; pi0=0.0; p=6.283185306/(1.0*n);pr1

13、=cos(p); pi1=-sin(p);if (1!=0) pi1=-pi1;for (i=2; i=n-1; i+) p=pri-1*pr1; q=pii-1*pi1; s=(pri-1+pii-1)*(pr1+pi1);pri=p-q; pii=s-p-q;for (it=0; it=0; l0-) m=m/2; nv=2*nv;for (it=0; it=(m-1)*nv; it=it+nv)for (j=0; j=(nv/2)-1; j+) p=prm*j*frit+j+nv/2; q=pim*j*fiit+j+nv/2; s=prm*j+pim*j; s=s*(frit+j+nv/2+fiit+j+nv/2); poddr=p-q; poddi=s-p-q; frit+j+nv/2=frit+j-poddr; fiit+j+nv/2=fiit+j-poddi; frit+j=frit+j+poddr; fiit+j=fiit+j+poddi;if (l!=0)for (i=0; i=n-1; i+) fri=fri/(1.0