電子發燒友App
傅立葉變換的重要性不用我說，想必大家也很清楚，有了傅立葉變換，我們就可以從信號的頻域特征去分析信號。尤其在無線通信系統中，傅里葉變換的重要性就更加明顯了，無論是設計者還是測試工程師，在工作中都會和傅立葉變換打交道。在以下的文章中，我給出一種傅里葉變換的C語言實現方法（參考了C常用算法集），可以用于在嵌入式系統中實現傅立葉變換。
常規的傅立葉變換算法并不適用于嵌入式控制系統，原因是運算量太大（涉及到復數運算），比如離散的傅立葉變換等同于用序列Y(n×1列矢量)乘以n×n 矩陣Fn，需要n×n次乘法。若n=1024,則是104，8576次乘法運算。哇，這么多呀！什么概念呢？如果你選用的CPU單周期指令為25ns,?單周期也可以完成一次乘法運算，那么要計算1024點的傅立葉變換則需要26.2144ms,這還不包括加法或其它運算，對于大多數實時系統，這個處理時間實在太長。于是尋找一個快速的傅立葉變換算法是人們所期望的。
本來我想把FFT的整個數學推導過程列完出來，但當自己硬著頭皮看完后，發現對我沒有任何用處，我又不是專門研究數學算法的，哪有那么多時間跟著書本的公式去慢慢推導。我想，這些推導問題還是讓數學家想去吧。我需要的不過是理解它，然后學會應用它就行。有興趣的讀者可以參考相關的資料，這方面的資料實在太多了。
雖然FFT大幅度地降低了常規傅立葉變換的運算量，但對于一般的單片機而言，處理FFT運算還是力不從心。主要原因是FFT計算過程中的蝶形運算是復數運算，要分開實部和虛部分別計算，想想這是多么繁瑣的事情?？赡軙行┏鯇W者認為，有這么復雜嗎？我在PC上使用C++一樣可以對復數直接進行加、減、乘、除運算。你說得不錯，可以這么做，但那是C++封裝了對復數處理的類，直接調用就行。在PC上運算這種類型的算法一般不考慮時間和空間，多一兩秒的運行時間不會有什么災難性的結果。
所以我們要衡量一個處理器有沒有足夠的能力來運行FFT算法，根據以上的簡單介紹可以得出以下兩點：
處理器要在一個指令周期能完成乘和累加的工作，因為復數運算要多次查表相乘才能實現。
間接尋址，可以實現增/減1個變址量，方便各種查表方法。FFT要對原始序列進行反序排列，處理器要有反序間接尋址的能力。
所以，在數字信號的分析處理應用中，DSP比其它的處理器有絕對的優勢，因為DSP完全具備以上條件。這就是單片機（51系列，AVR，PIC等等）或ARM處理器很少用來進行數字信號分析的原因。
重點來了，下面的這段程序就是用C語言實現傅里葉變換
//**********************************************************
??????????????????????? // 函數名: 快速傅立葉變換（來源《C常用算法集》）
??????????????????????? // 本函數測試OK,可以在TC2.0,VC++6.0,Keil C51測試通過。
??????????????????????? // 如果你的MCS51系統有足夠的RAM時,可以驗證一下用單片機處理FFT有多么的慢。
??????????????????????? //
??????????????????????? // 入口參數：
??????????????????????? // l: l = 0, 傅立葉變換; l = 1, 逆傅立葉變換
??????????????????????? // il: il = 0,不計算傅立葉變換或逆變換模和幅角；il = 1,計算模和幅角
??????????????????????? // n: 輸入的點數，為偶數，一般為32，64，128，…,1024等
??????????????????????? // k: 滿足n=2^k(k>0),實質上k是n個采樣數據可以分解為偶次冪和奇次冪的次數
??????????????????????? // pr[]: l=0時，存放N點采樣數據的實部
??????????????????????? // l=1時, 存放傅立葉變換的N個實部
??????????????????????? // pi[]: l=0時，存放N點采樣數據的虛部
??????????????????????? // l=1時, 存放傅立葉變換的N個虛部
??????????????????????? //
??????????????????????? // 出口參數：
??????????????????????? // fr[]: l=0, 返回傅立葉變換的實部
??????????????????????? // l=1, 返回逆傅立葉變換的實部
??????????????????????? // fi[]: l=0, 返回傅立葉變換的虛部
??????????????????????? // l=1, 返回逆傅立葉變換的虛部
??????????????????????? // pr[]: il = 1,i = 0 時，返回傅立葉變換的模
??????????????????????? // il = 1,i = 1 時，返回逆傅立葉變換的模
??????????????????????? // pi[]: il = 1,i = 0 時，返回傅立葉變換的輻角
??????????????????????? // il = 1,i = 1 時，返回逆傅立葉變換的輻角
??????????????????????? // 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;
??????????????????????? }
工商網監
評論