cfft演算法實現

A. 求用C++實現的FFT演算法

很早以前的，如果管用別忘了給我加分呀
/*
This computes an in-place complex-to-complex FFT
x and y are the real and imaginary arrays of 2^m points.
dir = 1 gives forward transform
dir = -1 gives reverse transform
*/
short FFT(short int dir,long m,double *x,double *y)
{
long n,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;

/* Calculate the number of points */
n = 1;
for (i=0;i<m;i++)
n *= 2;

/* Do the bit reversal */
i2 = n >> 1;
j = 0;
for (i=0;i<n-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}

/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<n;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}

/* Scaling for forward transform */
if (dir == 1) {
for (i=0;i<n;i++) {
x[i] /= n;
y[i] /= n;
}
}

return(TRUE);
}

---------------------------------------------------------------------------------
/*****************fft programe*********************/
#include "typedef.h"
#include "math.h"

struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,nv2,nm1,i,k,j=1,l;
/*int f,m,nv2,nm1,i,k,j=N/2,l;*/
struct compx v,w,t;
nv2=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++){;}
nm1=N-1;

/*變址運算*/
for(i=1;i <=nm1;i++)
{
if(i <j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=nv2;
while(k <j){j=j-k;k=k/2;}
j=j+k;
}

{
int le,lei,ip;
float pi;
for(l=1;l <=m;l++)
{ le=pow(2,l);// 這里用的是L而不是1 !!!!
lei =le/2;
pi=3.14159;
v.real=1.0;
v.imag=0.0;
w.real=cos(pi/lei);
w.imag=-sin(pi/lei);
for(j=1;j <=lei;j++)
{
/*double p=pow(2,m-l)*j;
double ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);*/
for(i=j;i <=N;i=i+le)
{ /* w.real=cos(ps);
w.imag=-sin(ps);*/
ip=i+lei;
t=EE(xin[ip],v);
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
v=EE(v,w);
}
}
}
return;
}

/*****************main programe********************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "typedef.h"

float result[257];
struct compx s[257];
int Num=256;
const float pp=3.14159;

main()
{

int i=1;
for(;i <0x101;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0;
}

FFT(s,Num);

for(i=1;i <0x101;i++)
{
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
}

}

-----------------------------------------------------------------------------------
FFT變換 C源代碼

FFT C source code (Simple radix-2)

void fft_float (
unsigned NumSamples,
int InverseTransform,
float *RealIn,
float *ImagIn,
float *RealOut,
float *ImagOut )
{
unsigned NumBits; /* Number of bits needed to store indices */
unsigned i, j, k, n;
unsigned BlockSize, BlockEnd;
double angle_numerator = 2.0 * DDC_PI;
double tr, ti; /* temp real, temp imaginary */
if ( !IsPowerOfTwo(NumSamples) )
{
fprintf (
stderr,
"Error in fft(): NumSamples=%u is not power of two\n",
NumSamples );
exit(1);
}
if ( InverseTransform )
angle_numerator = -angle_numerator;
CHECKPOINTER ( RealIn );
CHECKPOINTER ( RealOut );
CHECKPOINTER ( ImagOut );
NumBits = NumberOfBitsNeeded ( NumSamples );
/*
** Do simultaneous data and bit-reversal ordering into outputs...
*/
for ( i=0; i < NumSamples; i++ )
{
j = ReverseBits ( i, NumBits );
RealOut[j] = RealIn;
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn;
}
/*
** Do the FFT itself...
*/
BlockEnd = 1;
for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
{
double delta_angle = angle_numerator / (double)BlockSize;
double sm2 = sin ( -2 * delta_angle );
double sm1 = sin ( -delta_angle );
double cm2 = cos ( -2 * delta_angle );
double cm1 = cos ( -delta_angle );
double w = 2 * cm1;
double ar[3], ai[3];
double temp;
for ( i=0; i < NumSamples; i += BlockSize )
{
ar[2] = cm2;
ar[1] = cm1;
ai[2] = sm2;
ai[1] = sm1;
for ( j=i, n=0; n < BlockEnd; j++, n++ )
{
ar[0] = w*ar[1] - ar[2];
ar[2] = ar[1];
ar[1] = ar[0];
ai[0] = w*ai[1] - ai[2];
ai[2] = ai[1];
ai[1] = ai[0];
k = j + BlockEnd;
tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];
ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];
RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;
RealOut[j] += tr;
ImagOut[j] += ti;
}
}
BlockEnd = BlockSize;
}
/*
** Need to normalize if inverse transform...
*/
if ( InverseTransform )
{
double denom = (double)NumSamples;
for ( i=0; i < NumSamples; i++ )
{
RealOut /= denom;
ImagOut /= denom;
}
}
}

int IsPowerOfTwo ( unsigned x )
{
if ( x < 2 )
return FALSE;
if ( x & (x-1) ) // Thanks to 'byang' for this cute trick!
return FALSE;
return TRUE;
}

unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo )
{
unsigned i;
if ( PowerOfTwo < 2 )
{
fprintf (
stderr,
">>> Error in fftmisc.c: argument %d to NumberOfBitsNeeded is too small.\n",
PowerOfTwo );
exit(1);
}
for ( i=0; ; i++ )
{
if ( PowerOfTwo & (1 << i) )
return i;
}
}

unsigned ReverseBits ( unsigned index, unsigned NumBits )
{
unsigned i, rev;
for ( i=rev=0; i < NumBits; i++ )
{
rev = (rev << 1) | (index & 1);
index >>= 1;
}
return rev;
}

double Index_to_frequency ( unsigned NumSamples, unsigned Index )
{
if ( Index >= NumSamples )
return 0.0;
else if ( Index <= NumSamples/2 )
return (double)Index / (double)NumSamples;
return -(double)(NumSamples-Index) / (double)NumSamples;
}

B. 在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上實現的程序~

C. FFT的演算法

FFT是一種DFT的高效演算法，稱為快速傅立葉變換（fast Fourier transform），它根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的演算法進行改進獲得的。FFT演算法可分為按時間抽取演算法和按頻率抽取演算法，先簡要介紹FFT的基本原理。從DFT運算開始，說明FFT的基本原理。
DFT的運算為：

式中
由這種方法計算DFT對於 的每個K值，需要進行4N次實數相乘和（4N-2）次相加，對於N個k值，共需4N*4N次實數相乘和（4N-2)(4N-2）次實數相加。改進DFT演算法，減小它的運算量，利用DFT中 的周期性和對稱性，使整個DFT的計算變成一系列迭代運算，可大幅度提高運算過程和運算量，這就是FFT的基本思想。
FFT對傅氏變換的理論並沒有新的發現，但是對於在計算機系統或者說數字系統中應用離散傅立葉變換，可以說是進了一大步。
設x(n)為N項的復數序列，由DFT變換，任一X（m）的計算都需要N次復數乘法和N-1次復數加法，而一次復數乘法等於四次實數乘法和兩次實數加法，一次復數加法等於兩次實數加法，即使把一次復數乘法和一次復數加法定義成一次「運算」（四次實數乘法和四次實數加法），那麼求出N項復數序列的X（m）,即N點DFT變換大約就需要N^2次運算。當N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列（設N=2k,k為正整數），分為兩個N/2項的子序列，每個N/2點DFT變換需要（N/2）2次運算，再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以後，總的運算次數就變成N+2*（N/2)^2=N+（N^2）/2。繼續上面的例子，N=1024時，總的運算次數就變成了525312次，節省了大約50%的運算量。而如果我們將這種「一分為二」的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那麼N點的DFT變換就只需要Nlog2N次的運算，N在1024點時，運算量僅有10240次，是先前的直接演算法的1%，點數越多，運算量的節約就越大，這就是FFT的優越性。

D. 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變換。

E. FFT運算，在信號處理中是怎樣運用的啊

FFT演算法實現因為PIC16F877片內有高達368×8位（相當於184×16位）的數據存儲器（RAM），故用片內RAM最多可以完成64點FFT（16位實部和虛部數據）。現在僅實現16點FFT，主要是起拋磚引玉的作用。這里的FFT是按頻率抽取的。在調用FFT子程序前，輸入數據按正常次序輸入，而輸出數據是經FFT變換整序處理後輸出。原始數據被變換後的數據覆蓋存放在RAM中，這是通過分解序列實現的；然而分解序列將引起DFT的項序混亂，所以在變換結束，所有的數據需要進行「整序」,以恢復DFT的正常次序。某些應用可以不進行整序；因而整序程序編成子程序形式，當需要時隨時可以調用。輸入數據為32位，前為16位實部，後為16位虛部，中間結果為32位；輸出數據也是前為實部，後為虛部。這樣計算的結果具有相當高的精度。

F. 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)來表示

G. 實序列的FFT演算法

在以上討論FFT演算法中，均假定序列x(l)為復的，但實際問題中的序列大多為實的。當然，我們可以把實序列處理成虛部為零的復序列。因此，就要引進許多零參加運算。這樣一來，在機器運算時間和存儲單元方面都將造成很大的浪費。在本段中，我們介紹對實序列x(l)應用FFT演算法的一個有效方法。

1.同時計算兩個實序列的FFT演算法

設有N=4的兩個實序列x₁(l)與x₂(l)。為了求得它們的譜X₁(m)與X₂(m)，我們用此二實序列構造成如下復序列

物探數字信號分析與處理技術

利用上一段的方法，可以求得復序列x(l)的譜X(m)。根據(7-3-1)得到

物探數字信號分析與處理技術

上式中的m用N-m代替，則得

物探數字信號分析與處理技術

將上式兩端取共軛，根據對稱性有

物探數字信號分析與處理技術

根據DFT的復共軛性質，對於實序列x₁(l)與x₂(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個點的兩個子序列

物探數字信號分析與處理技術

設x₁(l)的譜與x₂(l)的譜分別為X₁(m)與X₂(m)

物探數字信號分析與處理技術

其中

於是可以將實序列x(l)的譜X(m)，用兩個子序列x₁(l)，x₂(l)的譜X₁(m)，X₂(m)來表示

物探數字信號分析與處理技術

其中

物探數字信號分析與處理技術

注意，x₁(l)，x₂(l)與X₁(m)，X₂(m)均以N/2為周期，

利用x₁(l)、x₂(l)構成如下復序列

物探數字信號分析與處理技術

利用FFT演算法可以求得復序列 的譜 。根據(7-3-7)就求得兩個實子序列的譜X₁(m)與X₂(m)

物探數字信號分析與處理技術

有了X₁(m)，X₂(m)，根據(7-3-10)就可求得X(m)。以上就是用FFT演算法求實序列x(l)的譜X(m)的方法。必須指出，用公式(7-3-10)求X(m)時，第一，兩個實子序列的譜X₁(m)，X₂(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)得到

物探數字信號分析與處理技術

H. 怎麼用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)

I. fft演算法matlab的實現代碼！完整版的！

function result = MyFFT(vector)
result = fft(vector);

J. 基2—fft演算法的軟體實現(MATLAB代碼)

參考網路： clc; clear all; close all; x=ones(1,128); %輸入的信號，自己可以改變 %整體運用原位計算 m=nextpow2(x);N=2^m; % 求x的長度對應的2的最低冪次m if length(x)<N x=[x,zeros(1,N-length(x))]; % 若x的長度不是2的冪，補零到2的整數冪 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m數列序號的倒序 y=x(nxd); % 將x倒序排列作為y的初始值 for mm=1:m % 將DFT作m次基2分解,從左到右，對每次分解作DFT運算,共做m級蝶形運算，每一級都有2^(mm-1)個蝶形結 Nz=2^mm;u=1; % 旋轉因子u初始化為WN^0=1 WN=exp(-i*2*pi/Nz); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nz) for j=1:Nz/2 % 本次跨越間隔內的各次蝶形運算,在進行第mm級運算時需要2^(mm-1)個 蝶形 for k=j:Nz:N % 本次蝶形運算的跨越間隔為Nz=2^mm kp=k+Nz/2; % 蝶形運算的兩個因子對應單元下標的關系 t=y(kp)*u; % 蝶形運算的乘積項 y(kp)=y(k)-t; % 蝶形運算 y(k)=y(k)+t; % 蝶形運算 end u=u*WN; % 修改旋轉因子,多乘一個基本DFT因子WN end end y y1=fft(x) %自己編的FFT跟直接調用的函數運算以後的結果進行對比