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

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

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

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

4、葉變換。最初傅里葉分析是作為熱過程的解析分析的工具被提出的（參見：林家翹、西格爾著自然科學中確定性問題的應用數(shù)學，科學出版社，北京。原版書名為 c. c. lin & l. a. segel, mathematics applied to deterministic problems in the natural sciences, macmillan inc., new york, 1974）。 傅里葉變換屬于諧波分析。 傅里葉變換的逆變換容易求出,而且形式與正變換非常類似; 正弦基函數(shù)是微分運算的本征函數(shù),從而使得線性微分方程的求解可以轉化為常系數(shù)的代數(shù)方程的求解.在線性時不變的物理系統(tǒng)內

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

6、: 滿足n=2k(k0),實質上k是n個采樣數(shù)據(jù)可以分解為偶次冪和奇次冪的次數(shù)/ pr: l=0時，存放n點采樣數(shù)據(jù)的實部/ l=1時, 存放傅立葉變換的n個實部/ pi: l=0時，存放n點采樣數(shù)據(jù)的虛部 / l=1時, 存放傅立葉變換的n個虛部/ 出口參數(shù)：/ fr: l=0, 返回傅立葉變換的實部/ l=1, 返回逆傅立葉變換的實部/ fi: l=0, 返回傅立葉變換的虛部/ l=1, 返回逆傅立葉變換的虛部/ pr: il = 1,l = 0 時，返回傅立葉變換的模/ il = 1,l = 1 時，返回逆傅立葉變換的模/ pi: il = 1,l = 0 時，返回傅立葉變換的輻角/ i

7、l = 1,l = 1 時，返回逆傅立葉變換的輻角void kbfft(double *pr,double *pi,int n,int k,double *fr,double *fi,int l,int il) int it,m,is,i,j,nv,l0; double 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; /蝶形運算 pr0=1.0; pi0=0.0; p=6

8、.283185306/(1.0*n); pr1=cos(p); pi1=-sin(p); if (l!=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

9、=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: 輸入的點數(shù)，為偶數(shù)，一般為32，64，128，.,1024等/ k: 滿足n=2k(k0),實質上k是

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

12、立葉變換的輻角void kbfft(double *pr,double *pi,int n,int k,double *fr,double *fi,int l,int il) int it,m,is,i,j,nv,l0; double 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; /蝶形運算 pr0=1.0; pi0=0.0; p=6.283185306/(1.0*n)

13、; pr1=cos(p); pi1=-sin(p); if (l!=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*n)