fft演算法語言

1. 怎樣用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;
}

2. FFT的最優演算法是什麼以及其代碼（C語言），謝謝！

應該是庫利-圖基演算法和桑德-圖基演算法吧。這兩種演算法的時間復雜度是一樣的，需要（N/2)log2N次的復數乘法和Nlog2N的復數加法。
當然你要是用基-4的FFT會更快，需要3/8Nlog2N次的復數乘法和Nlog2N次的加法。但這樣做的一個很麻煩的事是在做快速傅立葉變換時需要將原數據補足到2或4的整數次方。因此如果數據量合適的話基-4要快，如果數據不合適還是用基-2好。至於C語言代碼暫時沒有。還有為什麼要編C啊？用Matlab不是更好嗎？連循環都不用寫，甚至還有已經寫好的函數fft（），直接看這個函數演算法就好了

3. 基於FFT的演算法優化 要C語言完整程序（利用旋轉因子的性質），有的請留言，答謝！！！(有核心代碼，望指教

實現(C描述)

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

//#include "complex.h"

// --------------------------------------------------------------------------

#define N 8 //64

#define M 3 //6 //2^m=N

#define PI 3.1415926

// --------------------------------------------------------------------------

float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};

float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};

float x_i[N]; //N=8

/*

float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,

0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,

0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,

-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64

float x_r[N]={1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,};

float x_i[N];

*/

FILE *fp;

// ----------------------------------- func -----------------------------------

/**

* 初始化輸出虛部

*/

static void fft_init( void )

{

int i;

for(i=0; i<N; i++) x_i[i] = 0.0;

}

/**

* 反轉演算法.將時域信號重新排序.

* 這個演算法有改進的空間

*/

static void bitrev( void )

{

int p=1, q, i;

int bit_rev[ N ]; //

float xx_r[ N ]; //

bit_rev[ 0 ] = 0;

while( p < N )

{

for(q=0; q<p; q++)

{

bit_rev[ q ] = bit_rev[ q ] * 2;

bit_rev[ q + p ] = bit_rev[ q ] + 1;

}

p *= 2;

}

for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];

for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];

}

/* ------------ add by sshc625 ------------ */

static void bitrev2( void )

{

return ;

}

/* */

void display( void )

{

printf("\n\n");

int i;

for(i=0; i<N; i++)

printf("%f\t%f\n", x_r[i], x_i[i]);

}

/**

*

*/

void fft1( void )

{ fp = fopen("log1.txt", "a+");

int L, i, b, j, p, k, tx1, tx2;

float TR, TI, temp; // 臨時變數

float tw1, tw2;

/* 深M. 對層進行循環. L為當前層, 總層數為M. */

for(L=1; L<=M; L++)

{

fprintf(fp,"----------Layer=%d----------\n", L);

/* b的意義非常重大,b表示當前層的顆粒具有的輸入樣本點數 */

b = 1;

i = L - 1;

while(i > 0)

{

b *= 2;

i--;

}

// -------------- 是否外層對顆粒循環, 內層對樣本點循環邏輯性更強一些呢! --------------

/*

* outter對參與DFT的樣本點進行循環

* L=1, 循環了1次(4個顆粒, 每個顆粒2個樣本點)

* L=2, 循環了2次(2個顆粒, 每個顆粒4個樣本點)

* L=3, 循環了4次(1個顆粒, 每個顆粒8個樣本點)

*/

for(j=0; j<b; j++)

{

/* 求旋轉因子tw1 */

p = 1;

i = M - L; // M是為總層數, L為當前層.

while(i > 0)

{

p = p*2;

i--;

}

p = p * j;

tx1 = p % N;

tx2 = tx1 + 3*N/4;

tx2 = tx2 % N;

// tw1是cos部分, 實部; tw2是sin部分, 虛數部分.

tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];

tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];

/*

* inner對顆粒進行循環

* L=1, 循環了4次(4個顆粒, 每個顆粒2個輸入)

* L=2, 循環了2次(2個顆粒, 每個顆粒4個輸入)

* L=3, 循環了1次(1個顆粒, 每個顆粒8個輸入)

*/

for(k=j; k<N; k=k+2*b)

{

TR = x_r[k]; // TR就是A, x_r[k+b]就是B.

TI = x_i[k];

temp = x_r[k+b];

/*

* 如果復習一下 (a+j*b)(c+j*d)兩個復數相乘後的實部虛部分別是什麼

* 就能理解為什麼會如下運算了, 只有在L=1時候輸入才是實數, 之後層的

* 輸入都是復數, 為了讓所有的層的輸入都是復數, 我們只好讓L=1時候的

* 輸入虛部為0

* x_i[k+b]*tw2是兩個虛數相乘

*/

fprintf(fp, "tw1=%f, tw2=%f\n", tw1, tw2);

x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;

x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;

x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;

x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;

fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k, x_r[k], x_i[k]);

fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k+b, x_r[k+b], x_i[k+b]);

} //

} //

} //

}

/**

* ------------ add by sshc625 ------------

* 該實現的流程為

* for( Layer )

* for( Granule )

* for( Sample )

*

*

*

*

*/

void fft2( void )

{ fp = fopen("log2.txt", "a+");

int cur_layer, gr_num, i, k, p;

float tmp_real, tmp_imag, temp; // 臨時變數, 記錄實部

float tw1, tw2;// 旋轉因子,tw1為旋轉因子的實部cos部分, tw2為旋轉因子的虛部sin部分.

int step; // 步進

int sample_num; // 顆粒的樣本總數(各層不同, 因為各層顆粒的輸入不同)

/* 對層循環 */

for(cur_layer=1; cur_layer<=M; cur_layer++)

{

/* 求當前層擁有多少個顆粒(gr_num) */

gr_num = 1;

i = M - cur_layer;

while(i > 0)

{

i--;

gr_num *= 2;

}

/* 每個顆粒的輸入樣本數N' */

sample_num = (int)pow(2, cur_layer);

/* 步進. 步進是N'/2 */

step = sample_num/2;

/* */

k = 0;

/* 對顆粒進行循環 */

for(i=0; i<gr_num; i++)

{

/*

* 對樣本點進行循環, 注意上限和步進

*/

for(p=0; p<sample_num/2; p++)

{

// 旋轉因子, 需要優化...

tw1 = cos(2*PI*p/pow(2, cur_layer));

tw2 = -sin(2*PI*p/pow(2, cur_layer));

tmp_real = x_r[k+p];

tmp_imag = x_i[k+p];

temp = x_r[k+p+step];

/*(tw1+jtw2)(x_r[k]+jx_i[k])

*

* real : tw1*x_r[k] - tw2*x_i[k]

* imag : tw1*x_i[k] + tw2*x_r[k]

* 我想不抽象出一個

* typedef struct {

* double real; // 實部

* double imag; // 虛部

* } complex; 以及針對complex的操作

* 來簡化復數運算是否是因為效率上的考慮!

*/

/* 蝶形演算法 */

x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );

x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );

/* X[k] = A(k)+WB(k)

* X[k+N/2] = A(k)-WB(k) 的性質可以優化這里*/

// 旋轉因子, 需要優化...

tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));

tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));

x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );

x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );

printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);

printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);

}

/* 開跳!:) */

k += 2*step;

}

}

}

/*

* 後記:

* 究竟是顆粒在外層循環還是樣本輸入在外層, 好象也差不多, 復雜度完全一樣.

* 但以我資質愚鈍花費了不少時間才弄明白這數十行代碼.

* 從中我發現一個於我非常有幫助的教訓, 很久以前我寫過一部分演算法, 其中絕大多數都是遞歸.

* 將數據量減少, 減少再減少, 用歸納的方式來找出數據量加大代碼的規律

* 比如FFT

* 1. 先寫死LayerI的代碼; 然後再把LayerI的輸出作為LayerII的輸入, 又寫死代碼; ......

* 大約3層就可以統計出規律來. 這和遞歸也是一樣, 先寫死一兩層, 自然就出來了!

* 2. 有的功能可以寫偽代碼, 不急於求出結果, 降低復雜性, 把邏輯結果定出來後再添加.

* 比如旋轉因子就可以寫死, 就寫1.0. 流程出來後再寫旋轉因子.

* 寥寥數語, 我可真是流了不少汗! Happy!

*/

void dft( void )

{

int i, n, k, tx1, tx2;

float tw1,tw2;

float xx_r[N],xx_i[N];

/*

* clear any data in Real and Imaginary result arrays prior to DFT

*/

for(k=0; k<=N-1; k++)

xx_r[k] = xx_i[k] = x_i[k] = 0.0;

// caculate the DFT

for(k=0; k<=(N-1); k++)

{

for(n=0; n<=(N-1); n++)

{

tx1 = (n*k);

tx2 = tx1+(3*N)/4;

tx1 = tx1%(N);

tx2 = tx2%(N);

if(tx1 >= (N/2))

tw1 = -twiddle[tx1-(N/2)];

else

tw1 = twiddle[tx1];

if(tx2 >= (N/2))

tw2 = -twiddle[tx2-(N/2)];

else

tw2 = twiddle[tx2];

xx_r[k] = xx_r[k]+x_r[n]*tw1;

xx_i[k] = xx_i[k]+x_r[n]*tw2;

}

xx_i[k] = -xx_i[k];

}

// display

for(i=0; i<N; i++)

printf("%f\t%f\n", xx_r[i], xx_i[i]);

}

// ---------------------------------------------------------------------------

int main( void )

{

fft_init( );

bitrev( );

// bitrev2( );

//fft1( );

fft2( );

display( );

system( "pause" );

// dft();

return 1;

}

本文來自CSDN博客，轉載請標明出處：http://blog.csdn.net/sshcx/archive/2007/06/14/1651616.aspx

4. FFT的公式是什麼和演算法是怎樣實現

二維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演算法的實現。下面討論一維FFT的演算法原理。

【1D-FFT的演算法實現】
設序列h(n)長度為N，將其按下標的奇偶性分成兩組，即he和ho序列，它們的長度都是N/2。這樣，可以將h(n)的FFT計算公式改寫如下 ：

(A)
由於

所以，(A)式可以改寫成下面的形式：

按照FFT的定義，上面的式子實際上是：

其中，k的取值范圍是 0~N-1。
我們注意到He(k)和Ho(k)是N/2點的DFT，其周期是N/2。因此，H(k)DFT的前N/2點和後N/2點都可以用He(k)和Ho(k)來表示

5. matlab中fft（）函數是什麼意思

FFT（快速傅里葉變換）是一種實現DFT（離散傅里葉變換）的快速演算法，是利用復數形式的離散傅里葉變換來計算實數形式的離散傅里葉變換，matlab中的fft()函數是實現該演算法的實現。

MATLAB它將數值分析、矩陣計算、科學數據可視化以及非線性動態系統的建模和模擬等諸多強大功能集成在一個易於使用的視窗環境中，為科學研究、工程設計以及必須進行有效數值計算的眾多科學領域提供了一種全面的解決方案，並在很大程度上擺脫了傳統非互動式程序設計語言（如C、Fortran）的編輯模式，代表了當今國際科學計算軟體的先進水平。

快速傅里葉變換, 即利用計算機計算離散傅里葉變換（DFT)的高效、快速計算方法的統稱，簡稱FFT。快速傅里葉變換是1965年由J.W.庫利和T.W.圖基提出的。採用這種演算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數N越多，FFT演算法計算量的節省就越顯著。

(5)fft演算法語言擴展閱讀：

matlab優勢特點：

1、高效的數值計算及符號計算功能，能使用戶從繁雜的數學運算分析中解脫出來；

2、具有完備的圖形處理功能，實現計算結果和編程的可視化；

3、友好的用戶界面及接近數學表達式的自然化語言，使學者易於學習和掌握；

4、功能豐富的應用工具箱(如信號處理工具箱、通信工具箱等) ，為用戶提供了大量方便實用的處理工具。

參考資料來源：

網路-快速傅里葉變換

網路-MATLAB

6. FFT 演算法

FFT即Fast Fourier Transform，中文翻譯：快速傅立葉演算法。下面是網上找到的演算法實現。留以備用。

7. FFT演算法分幾種

FFT演算法分析FFT演算法的基本原理是把長序列的DFT逐次分解為較短序列的DFT。按照抽取方式的不同可分為DIT-FFT（按時間抽取）和DIF-FFT（按頻率抽取）演算法。按照蝶形運算的構成不同可分為基2、基4、基8以及任意因子（2n,n為大於1的整數），基2、基4演算法較為常用。 網上有幫助文檔: http://www.5doc.com/doc/123035(右上角有點擊下載)

8. 什麼是FFT

計算離散傅里葉變換的一種快速演算法,簡稱FFT(Fast Fourier Transform)。快速傅里葉變換是1965年由J.W.庫利和T.W.圖基提出的。採用這種演算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少,特別是被變換的抽樣點數N越多,FFT演算法計算量的節省就越顯著。
FFT 的出現，使信號分析從時域分析向頻域分析成為可能，極大地推動了信號分析在各領域的實際應用。
http://ke..com/view/1006229.htm

9. FFT原理的FFT基本原理

FFT是一種DFT的高效演算法，稱為快速傅立葉變換（fast Fourier transform）。FFT演算法可分為按時間抽取演算法和按頻率抽取演算法，先簡要介紹FFT的基本原理。從DFT運算開始，說明FFT的基本原理。
DFT的運算為：

式中

由這種方法計算DFT對於X（K）的每個K值，需要進行4N次實數相乘和（4N-2）次相加，對於N個k值，共需N*N乘和N（4N-2）次實數相加。改進DFT演算法，減小它的運算量，利用DFT中

的周期性和對稱性，使整個DFT的計算變成一系列迭代運算，可大幅度提高運算過程和運算量，這就是FFT的基本思想。
FFT基本上可分為兩類，時間抽取法和頻率抽取法，而一般的時間抽取法和頻率抽取法只能處理長度N=2^M的情況，另外還有組合數基四FFT來處理一般長度的FFT 設N點序列x(n),,將x(n)按奇偶分組，公式如下圖

改寫為：

一個N點DFT分解為兩個 N/2點的DFT，繼續分解，迭代下去，其運算量約為

其演算法有如下規律
兩個4點組成的8點DFT

四個2點組成的8點DFT
按時間抽取的8點DFT
原位計算
當數據輸入到存儲器中以後，每一級運算的結果仍然儲存在同一組存儲器中，直到最後輸出，中間無需其它存儲器
序數重排
對按時間抽取FFT的原位運算結構，當運算完畢時，這種結構存儲單元A(1)、A(2)，…，A(8)中正好順序存放著X(0)，X(1)，X(2)，…，X(7)，因此可直接按順序輸出，但這種原位運算的輸入x(n)卻不能按這種自然順序存入存儲單元中，而是按X(0)，X(4)，X(2)，X(6)，…，X(7)的順序存入存儲單元，這種順序看起來相當雜亂，然而它也是有規律的。當用二進製表示這個順序時，它正好是「碼位倒置」的順序。
蝶形類型隨迭代次數成倍增加
每次迭代的蝶形類型比上一次蝶代增加一倍，數據點間隔也增大一倍 頻率抽取2FFT演算法是按頻率進行抽取的演算法。
設N=2^M，將x(n)按前後兩部分進行分解，
按K的奇偶分為兩組，即
得到兩個N/2 點的DFT運算。如此分解，並迭代，總的計算量和時間抽取（DIT）基2FFT演算法相同。
演算法規律如下：
蝶形結構和時間抽取不一樣但是蝶形個數一樣，同樣具有原位計算規律，其迭代次數成倍減小 時，可採取補零使其成為
，或者先分解為兩個p,q的序列，其中p*q=N，然後進行計算。 前面介紹，採用FFT演算法可以很快算出全部N點DFT值，即z變換X（z）在z平面單位圓上的全部等間隔取樣值。實際中也許①不需要計算整個單位圓上z變換的取樣，如對於窄帶信號，只需要對信號所在的一段頻帶進行分析，這時希望頻譜的采樣集中在這一頻帶內，以獲得較高的解析度，而頻帶以外的部分可不考慮，②或者對其它圍線上的z變換取樣感興趣，例如語音信號處理中，需要知道z變換的極點所在頻率，如極點位置離單位圓較遠，則其單位圓上的頻譜就很平滑，這時很難從中識別出極點所在的頻率，如果采樣不是沿單位圓而是沿一條接近這些極點的弧線進行，則在極點所在頻率上的頻譜將出現明顯的尖峰，由此可較准確地測定極點頻率。③或者要求能有效地計算當N是素數時序列的DFT，因此提高DFT計算的靈活性非常有意義。
螺旋線采樣是一種適合於這種需要的變換，且可以採用FFT來快速計算，這種變換也稱作Chirp-z變換。

10. 怎麼用C語言實現FFT演算法 呀

float ar[1024],ai[1024];/* 原始數據實部，虛部 */
float a[2050];

void fft(int nn) /* nn數據長度 */
{
int n1,n2,i,j,k,l,m,s,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,};

switch(nn)
{
case 1024: s=10; break;
case 512: s=9; break;
case 256: s=8; break;
}

n1=nn/2; n2=nn-1;
j=1;
for(i=1;i<=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;l<n2;l++)
{
if(l<j)
{
t1=a[2*j];
t2=a[2*j+1];
a[2*j]=a[2*l];
a[2*j+1]=a[2*l+1];
a[2*l]=t1;
a[2*l+1]=t2;
}
k=n1;
while (k<j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i<=s;i++)
{
u1=1;
u2=0;
m=(1<<i);
k=m>>1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j<=k;j++)
{
for(l=j;l<nn;l=l+m)
{
l1=l+k;
t1=a[2*l1]*u1-a[2*l1+1]*u2;
t2=a[2*l1]*u2+a[2*l1+1]*u1;
a[2*l1]=a[2*l]-t1;
a[2*l1+1]=a[2*l+1]-t2;
a[2*l]=a[2*l]+t1;
a[2*l+1]=a[2*l+1]+t2;
}
z=u1*w1-u2*w2;
u2=u1*w2+u2*w1;
u1=z;
}
}
for(i=1;i<=nn/2;i++)
{
ar[i]=4*a[2*i+2]/nn; /* 實部 */
ai[i]=-4*a[2*i+3]/nn; /* 虛部 */
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); /* 幅值 */
}
}

(http://..com/question/284943905.html?an=0&si=2)