fft演算法的實現
Ⅰ 怎樣用C語言實現FFT演算法啊
1、二維FFT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下：
先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩陣的第i行。注意這只是一個簡單的記法，並不能完全照抄。還需要通過一些語句來生成各行的數據。同理，COL[i]是對矩陣的第i列的一種簡單表示方法。
所以，關鍵是一維FFT演算法的實現。
2、常式：
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#defineN1000
/*定義復數類型*/
typedefstruct{
doublereal;
doubleimg;
}complex;
complexx[N],*W;/*輸入序列,變換核*/
intsize_x=0;/*輸入序列的大小，在本程序中僅限2的次冪*/
doublePI;/*圓周率*/
voidfft();/*快速傅里葉變換*/
voidinitW();/*初始化變換核*/
voidchange();/*變址*/
voidadd(complex,complex,complex*);/*復數加法*/
voidmul(complex,complex,complex*);/*復數乘法*/
voidsub(complex,complex,complex*);/*復數減法*/
voidoutput();
intmain(){
inti;/*輸出結果*/
system("cls");
PI=atan(1)*4;
printf("Pleaseinputthesizeofx: ");
scanf("%d",&size_x);
printf("Pleaseinputthedatainx[N]: ");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x[i].real,&x[i].img);
initW();
fft();
output();
return0;
}
/*快速傅里葉變換*/
voidfft(){
inti=0,j=0,k=0,l=0;
complexup,down,proct;
change();
for(i=0;i<log(size_x)/log(2);i++){/*一級蝶形運算*/
l=1<<i;
for(j=0;j<size_x;j+=2*l){/*一組蝶形運算*/
for(k=0;k<l;k++){/*一個蝶形運算*/
mul(x[j+k+l],W[size_x*k/2/l],&proct);
add(x[j+k],proct,&up);
sub(x[j+k],proct,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化變換核*/
voidinitW(){
inti;
W=(complex*)malloc(sizeof(complex)*size_x);
for(i=0;i<size_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*變址計算，將x(n)碼位倒置*/
voidchange(){
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
for(i=0;i<size_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while((t--)>0){
j=j<<1;
j|=(k&1);
k=k>>1;
}
if(j>i){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*輸出傅里葉變換的結果*/
voidoutput(){
inti;
printf("Theresultareasfollows ");
for(i=0;i<size_x;i++){
printf("%.4f",x[i].real);
if(x[i].img>=0.0001)printf("+%.4fj ",x[i].img);
elseif(fabs(x[i].img)<0.0001)printf(" ");
elseprintf("%.4fj ",x[i].img);
}
}
voidadd(complexa,complexb,complex*c){
c->real=a.real+b.real;
c->img=a.img+b.img;
}
voidmul(complexa,complexb,complex*c){
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
voidsub(complexa,complexb,complex*c){
c->real=a.real-b.real;
c->img=a.img-b.img;
}
Ⅱ matlab中fft（）函數是什麼意思
FFT（快速傅里葉變換）是一種實現DFT（離散傅里葉變換）的快速演算法，是利用復數形式的離散傅里葉變換來計算實數形式的離散傅里葉變換，matlab中的fft()函數是實現該演算法的實現。
MATLAB它將數值分析、矩陣計算、科學數據可視化以及非線性動態系統的建模和模擬等諸多強大功能集成在一個易於使用的視窗環境中，為科學研究、工程設計以及必須進行有效數值計算的眾多科學領域提供了一種全面的解決方案，並在很大程度上擺脫了傳統非互動式程序設計語言（如C、Fortran）的編輯模式，代表了當今國際科學計算軟體的先進水平。
快速傅里葉變換, 即利用計算機計算離散傅里葉變換（DFT)的高效、快速計算方法的統稱，簡稱FFT。快速傅里葉變換是1965年由J.W.庫利和T.W.圖基提出的。採用這種演算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數N越多，FFT演算法計算量的節省就越顯著。
(2)fft演算法的實現擴展閱讀：
matlab優勢特點：
1、高效的數值計算及符號計算功能，能使用戶從繁雜的數學運算分析中解脫出來；
2、具有完備的圖形處理功能，實現計算結果和編程的可視化；
3、友好的用戶界面及接近數學表達式的自然化語言，使學者易於學習和掌握；
4、功能豐富的應用工具箱(如信號處理工具箱、通信工具箱等) ，為用戶提供了大量方便實用的處理工具。
參考資料來源：
網路-快速傅里葉變換
網路-MATLAB
Ⅲ 在DSP上實現FFT演算法
void FFT( COMPLEX *Y, int N) /* input sample array, number of points */
{
COMPLEX temp1,temp2; /*temporary storage variables */蔽型
int i,j,k; /*loop counter variables */
int upper_leg, lower_leg; /*index of upper/lower butterfly leg */
int leg_diff; /宏握猜*difference between upper/lower leg */
int num_stages=0; /*number of FFT stages, or iterations */
int index, step; /*index and step between twiddle factor*/
/* log(base 2) of # of points = # of stages */
i=1;
do
{
num_stages+=1;
i = i *2 ;
} while (i!=N);
/* starting difference between upper and lower butterfly legs*/
leg_diff = N/2;
/* step between values in twiddle factor array twiddle.h */
step = 512 / N;
/* For N-point FFT */
for ( i = 0 ; i < num_stages ; i++ )
{
index = 0;
for ( j = 0; j < leg_diff ; j++ )
{
for ( upper_leg = j; upper_leg < N ; upper_leg += (2*leg_diff) )
{
lower_leg = upper_leg + leg_diff;
temp1.real=(Y[upper_leg]).real + (Y[lower_leg]).real;
temp1.imag=(Y[upper_leg]).imag + (Y[lower_leg]).imag;
temp2.real=(Y[upper_leg]).real - (Y[lower_leg]).real;
temp2.imag=(Y[upper_leg]).imag - (Y[lower_leg]).imag;
(Y[lower_leg]).real = ((long)temp2.real * (w[index]).real)/8192;
(Y[lower_leg]).real -= ((long)temp2.imag * (w[index]).imag)/8192;
(Y[lower_leg]).imag = ((long)temp2.real * (w[index]).imag)/8192;
(Y[lower_leg]).imag += ((long)temp2.imag * (w[index]).real)/8192;
(Y[upper_leg]).real = temp1.real;
(Y[upper_leg]).imag = temp1.imag;
}
index+=step;
}
leg_diff = leg_diff / 2;
step *= 2;
}
/* bit reversal for resequencing data */
j=0;
for ( i=1 ; i < (N-1) ; i++ )
{
k = N / 2;
while ( k <= j)
{
j = j - k;
k >>皮旅= 1;
}
j = j + k;
if ( i < j )
{
temp1.real = (Y[j]).real;
temp1.imag = (Y[j]).imag;
(Y[j]).real = (Y[i]).real;
(Y[j]).imag = (Y[i]).imag;
(Y[i]).real = temp1.real;
(Y[i]).imag = temp1.imag;
}
}
return;
}
參考一下的吧，這個是TI官方的在5416上實現的程序~
Ⅳ 實序列的FFT演算法
在以上討論FFT演算法中，均假定序列x(l)為復的，但實際問題中的序列大多為實的。當然，我們可以把實序列處理成虛部為零的復序列。因此，就要引進許多零參加運算。這樣一來，在機器運算時間和存儲單元方面都將造成很大的浪費。在本段中，我們介紹對實序列x(l)應用FFT演算法的一個有效方法。
1.同時計算兩個實序列的FFT演算法
設有N=4的兩個實序列x1(l)與x2(l)。為了求得它們的譜X1(m)與X2(m)，我們用此二實序列構造成如下復序列
物探數字信號分析與處理技術
利用上一段的方法，可以求得復序列x(l)的譜X(m)。根據(7-3-1)得到
物探數字信號分析與處理技術
上式中的m用N-m代替，則得
物探數字信號分析與處理技術
將上式兩端取共軛，根據對稱性有
物探數字信號分析與處理技術
根據DFT的復共軛性質，對於實序列x1(l)與x2(l)，有
物探數字信號分析與處理技術
於是從(7-3-4)得到
物探數字信號分析與處理技術
聯立求解(7-3-2)和(7-3-6)便得到
物探數字信號分析與處理技術
例如設有兩個N=4點的實序列，
物探數字信號分析與處理技術
我們用它們構造一個N=4點的復序列
物探數字信號分析與處理技術
利用FFT演算法求X(m)，m=0，1，2，3(圖7-3-1)，
圖7-3-1 N=4點的FFT演算法流程圖
於是得到
物探數字信號分析與處理技術
因此從式(7-3-7)得到
物探數字信號分析與處理技術
物探數字信號分析與處理技術
2.實序列的FFT演算法
設有N點的實序列x(l)，l=0，1，2，…，N-1。按照點的奇偶編號，將它們分成N/2個點的兩個子序列
物探數字信號分析與處理技術
設x1(l)的譜與x2(l)的譜分別為X1(m)與X2(m)
物探數字信號分析與處理技術
其中
於是可以將實序列x(l)的譜X(m)，用兩個子序列x1(l)，x2(l)的譜X1(m)，X2(m)來表示
物探數字信號分析與處理技術
其中
物探數字信號分析與處理技術
注意，x1(l)，x2(l)與X1(m)，X2(m)均以N/2為周期，
利用x1(l)、x2(l)構成如下復序列
物探數字信號分析與處理技術
利用FFT演算法可以求得復序列 的譜 。根據(7-3-7)就求得兩個實子序列的譜X1(m)與X2(m)
物探數字信號分析與處理技術
有了X1(m)，X2(m)，根據(7-3-10)就可求得X(m)。以上就是用FFT演算法求實序列x(l)的譜X(m)的方法。必須指出，用公式(7-3-10)求X(m)時，第一，兩個實子序列的譜X1(m)，X2(m)及復序列x珓(l)的譜珘X(m)均是以N/2為周期的周期序列;第二，由於x
(l)是實序列，根據DFT的復共軛性質有X(m)=X*(N-m)，m=0，1，…，N/2，故只需求得前(N/2)+1個點的X(m)，就得到全部N個點的X(m)了
例如，有N=8點的實序列，
物探數字信號分析與處理技術
首先，按點的奇偶編號分成兩個實子序列，
物探數字信號分析與處理技術
其次用它們構造如下復序列，
物探數字信號分析與處理技術
用FFT演算法求此復序列的譜 (圖7-3-2)
圖7-3-2 N=4點的FFT演算法流程圖
於是得到:
根據周期性，有
物探數字信號分析與處理技術
根據(7-3-12)式，
物探數字信號分析與處理技術
根據周期性，有
物探數字信號分析與處理技術
故最終由(7-3-10)得到
物探數字信號分析與處理技術