可用于嵌入式系統(tǒng)的傅立葉變換的C語(yǔ)言實(shí)現(xiàn)方法

傅立葉變換的重要性不用我說(shuō)，想必大家也很清楚，有了傅立葉變換，我們就可以從信號(hào)的頻域特征去分析信號(hào)。尤其在無(wú)線(xiàn)通信系統(tǒng)中，傅里葉變換的重要性就更加明顯了，無(wú)論是設(shè)計(jì)者還是測(cè)試工程師，在工作中都會(huì)和傅立葉變換打交道。在以下的文章中，我給出一種傅里葉變換的C語(yǔ)言實(shí)現(xiàn)方法（參考了C常用算法集），可以用于在嵌入式系統(tǒng)中實(shí)現(xiàn)傅立葉變換。

常規(guī)的傅立葉變換算法并不適用于嵌入式控制系統(tǒng)，原因是運(yùn)算量太大（涉及到復(fù)數(shù)運(yùn)算），比如離散的傅立葉變換等同于用序列Y(n×1列矢量)乘以n×n 矩陣Fn，需要n×n次乘法。若n=1024,則是104，8576次乘法運(yùn)算。哇，這么多呀！什么概念呢？如果你選用的CPU單周期指令為25ns,?單周期也可以完成一次乘法運(yùn)算，那么要計(jì)算1024點(diǎn)的傅立葉變換則需要26.2144ms,這還不包括加法或其它運(yùn)算，對(duì)于大多數(shù)實(shí)時(shí)系統(tǒng)，這個(gè)處理時(shí)間實(shí)在太長(zhǎng)。于是尋找一個(gè)快速的傅立葉變換算法是人們所期望的。

本來(lái)我想把FFT的整個(gè)數(shù)學(xué)推導(dǎo)過(guò)程列完出來(lái)，但當(dāng)自己硬著頭皮看完后，發(fā)現(xiàn)對(duì)我沒(méi)有任何用處，我又不是專(zhuān)門(mén)研究數(shù)學(xué)算法的，哪有那么多時(shí)間跟著書(shū)本的公式去慢慢推導(dǎo)。我想，這些推導(dǎo)問(wèn)題還是讓數(shù)學(xué)家想去吧。我需要的不過(guò)是理解它，然后學(xué)會(huì)應(yīng)用它就行。有興趣的讀者可以參考相關(guān)的資料，這方面的資料實(shí)在太多了。

雖然FFT大幅度地降低了常規(guī)傅立葉變換的運(yùn)算量，但對(duì)于一般的單片機(jī)而言，處理FFT運(yùn)算還是力不從心。主要原因是FFT計(jì)算過(guò)程中的蝶形運(yùn)算是復(fù)數(shù)運(yùn)算，要分開(kāi)實(shí)部和虛部分別計(jì)算，想想這是多么繁瑣的事情。可能會(huì)有些初學(xué)者認(rèn)為，有這么復(fù)雜嗎？我在PC上使用C++一樣可以對(duì)復(fù)數(shù)直接進(jìn)行加、減、乘、除運(yùn)算。你說(shuō)得不錯(cuò)，可以這么做，但那是C++封裝了對(duì)復(fù)數(shù)處理的類(lèi)，直接調(diào)用就行。在PC上運(yùn)算這種類(lèi)型的算法一般不考慮時(shí)間和空間，多一兩秒的運(yùn)行時(shí)間不會(huì)有什么災(zāi)難性的結(jié)果。

所以我們要衡量一個(gè)處理器有沒(méi)有足夠的能力來(lái)運(yùn)行FFT算法，根據(jù)以上的簡(jiǎn)單介紹可以得出以下兩點(diǎn)：

處理器要在一個(gè)指令周期能完成乘和累加的工作，因?yàn)閺?fù)數(shù)運(yùn)算要多次查表相乘才能實(shí)現(xiàn)。

間接尋址，可以實(shí)現(xiàn)增/減1個(gè)變址量，方便各種查表方法。FFT要對(duì)原始序列進(jìn)行反序排列，處理器要有反序間接尋址的能力。

所以，在數(shù)字信號(hào)的分析處理應(yīng)用中，DSP比其它的處理器有絕對(duì)的優(yōu)勢(shì)，因?yàn)镈SP完全具備以上條件。這就是單片機(jī)（51系列，AVR，PIC等等）或ARM處理器很少用來(lái)進(jìn)行數(shù)字信號(hào)分析的原因。

重點(diǎn)來(lái)了，下面的這段程序就是用C語(yǔ)言實(shí)現(xiàn)傅里葉變換

//**********************************************************
??????????????????????? 　　// 函數(shù)名: 快速傅立葉變換（來(lái)源《C常用算法集》）
??????????????????????? 　　// 本函數(shù)測(cè)試OK,可以在TC2.0,VC++6.0,Keil C51測(cè)試通過(guò)。
??????????????????????? 　　// 如果你的MCS51系統(tǒng)有足夠的RAM時(shí),可以驗(yàn)證一下用單片機(jī)處理FFT有多么的慢。
??????????????????????? 　　//
??????????????????????? 　　// 入口參數(shù)：
??????????????????????? 　　// l: l = 0, 傅立葉變換; l = 1, 逆傅立葉變換
??????????????????????? 　　// il: il = 0,不計(jì)算傅立葉變換或逆變換模和幅角；il = 1,計(jì)算模和幅角
??????????????????????? 　　// n: 輸入的點(diǎn)數(shù)，為偶數(shù)，一般為32，64，128，…,1024等
??????????????????????? 　　// k: 滿(mǎn)足n=2^k(k>0),實(shí)質(zhì)上k是n個(gè)采樣數(shù)據(jù)可以分解為偶次冪和奇次冪的次數(shù)
??????????????????????? 　　// pr[]: l=0時(shí)，存放N點(diǎn)采樣數(shù)據(jù)的實(shí)部
??????????????????????? 　　// l=1時(shí), 存放傅立葉變換的N個(gè)實(shí)部
??????????????????????? 　　// pi[]: l=0時(shí)，存放N點(diǎn)采樣數(shù)據(jù)的虛部
??????????????????????? 　　// l=1時(shí), 存放傅立葉變換的N個(gè)虛部
??????????????????????? 　　//
??????????????????????? 　　// 出口參數(shù)：
??????????????????????? 　　// fr[]: l=0, 返回傅立葉變換的實(shí)部
??????????????????????? 　　// l=1, 返回逆傅立葉變換的實(shí)部
??????????????????????? 　　// fi[]: l=0, 返回傅立葉變換的虛部
??????????????????????? 　　// l=1, 返回逆傅立葉變換的虛部
??????????????????????? 　　// pr[]: il = 1,i = 0 時(shí)，返回傅立葉變換的模
??????????????????????? 　　// il = 1,i = 1 時(shí)，返回逆傅立葉變換的模
??????????????????????? 　　// pi[]: il = 1,i = 0 時(shí)，返回傅立葉變換的輻角
??????????????????????? 　　// il = 1,i = 1 時(shí)，返回逆傅立葉變換的輻角
??????????????????????? 　　// data: 2005.8.15,Mend Xin Dong
??????????????????????? 　　void kkfft(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;
??????????????????????? 　　 }
??????????????????????? 　　 fr[it] = pr[is];
??????????????????????? 　　 fi[it] = pi[is];
??????????????????????? 　　 }
??????????????????????? 　　//—————————-
??????????????????????? 　　 pr[0] = 1.0;
??????????????????????? 　　 pi[0] = 0.0;
??????????????????????? 　　 p = 6.283185306/(1.0*n);
??????????????????????? 　　 pr[1] = cos(p);
??????????????????????? 　　 pi[1] = -sin(p);
???????????????????????
??????????????????????? 　　 if (l!=0)
??????????????????????? 　　 pi[1]=-pi[1];
???????????????????????
??????????????????????? 　　 for (i=2; i<=n-1; i++)
??????????????????????? 　　 {
??????????????????????? 　　 p = pr[i-1]*pr[1];
??????????????????????? 　　 q = pi[i-1]*pi[1];
??????????????????????? 　　 s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
??????????????????????? 　　 pr[i] = p-q;
??????????????????????? 　　 pi[i] = s-p-q;
??????????????????????? 　　 }
???????????????????????
??????????????????????? 　　 for (it=0; it<=n-2; it=it+2)
??????????????????????? 　　 {
??????????????????????? 　　 vr = fr[it];
??????????????????????? 　　 vi = fi[it];
??????????????????????? 　　 fr[it] = vr+fr[it+1];
??????????????????????? 　　 fi[it] = vi+fi[it+1];
??????????????????????? 　　 fr[it+1] = vr-fr[it+1];
??????????????????????? 　　 fi[it+1] = vi-fi[it+1];
??????????????????????? 　　 }
??????????????????????? 　　 m = n/2;
??????????????????????? 　　 nv = 2;
???????????????????????
??????????????????????? 　　 for (l0=k-2; l0>=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 = pr[m*j]*fr[it+j+nv/2];
??????????????????????? 　　 q = pi[m*j]*fi[it+j+nv/2];
??????????????????????? 　　 s = pr[m*j]+pi[m*j];
??????????????????????? 　　 s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
??????????????????????? 　　 poddr = p-q;
??????????????????????? 　　 poddi = s-p-q;
??????????????????????? 　　 fr[it+j+nv/2] = fr[it+j]-poddr;
??????????????????????? 　　 fi[it+j+nv/2] = fi[it+j]-poddi;
??????????????????????? 　　 fr[it+j] = fr[it+j]+poddr;
??????????????????????? 　　 fi[it+j] = fi[it+j]+poddi;
??????????????????????? 　　 }
??????????????????????? 　　 }
???????????????????????
??????????????????????? 　　 if(l!=0)
??????????????????????? 　　 for(i=0; i<=n-1; i++)
??????????????????????? 　　 {
??????????????????????? 　　 fr[i] = fr[i]/(1.0*n);
??????????????????????? 　　 fi[i] = fi[i]/(1.0*n);
??????????????????????? 　　 }
???????????????????????
??????????????????????? 　　 if(il!=0)
??????????????????????? 　　 for(i=0; i<=n-1; i++)
??????????????????????? 　　 {
??????????????????????? 　　 pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
??????????????????????? 　　 if(fabs(fr[i])<0.000001*fabs(fi[i]))
??????????????????????? 　　 {
??????????????????????? 　　 if ((fi[i]*fr[i])>0)
??????????????????????? 　　 pi[i] = 90.0;
??????????????????????? 　　 else
??????????????????????? 　　 pi[i] = -90.0;
??????????????????????? 　　 }
??????????????????????? 　　 else
??????????????????????? 　　 pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;
??????????????????????? 　　 }
??????????????????????? 　　return;
??????????????????????? 　　}