cfft算法实现

⑴ 实序列的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)得到

物探数字信号分析与处理技术

⑵ 怎样用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;
}

⑶ 怎么用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)

⑷ 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变换。

⑸ 在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算法实现因为PIC16F877片内有高达368×8位（相当于184×16位）的数据存储器（RAM），故用片内RAM最多可以完成64点FFT（16位实部和虚部数据）。现在仅实现16点FFT，主要神塌是起抛砖引玉的谨瞎嫌作用。这里的FFT是按频率抽取的。在调用FFT子程序前，输入数据按正常次序输入，而输出数据是经FFT变换整序处理后输出。原始数据被变换后的数据覆盖存放在RAM中，这是通过分解序列实现的；然而分解序列将引起DFT的项序混乱，所以在变换结束，所有的数据需要进行“整序”,以恢复DFT的正常次序。某些应用可以不进行整序；因而整序程序编成子程序形式，当需要时随时可以调用。输入数据为32位祥手，前为16位实部，后为16位虚部，中间结果为32位；输出数据也是前为实部，后为虚部。这样计算的结果具有相当高的精度。

⑺ 16点DFT的FFT算法

FFT（快速傅里叶变换）是DFT的一种特殊情况，就是当运算点的个数是2的整数次幂的时候进行的运算（不够用0补齐）。

FFT计算原理及流程图：

原理：FFT的计算要求点数必须为2的整数次幂，如果点数不够用0补齐。例如计算{2，3，5，8，4}的16点FFT，需要补11个0后进行计算。FFT计算运用蝶形运算，在蝶形运算中变化规律由W(N, p)推导，其中N为FFT计算点数，J为下角标的值。

L = 1时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0;

L = 2时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1;

L = 3时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3;

所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1

又因为2^L = 2^M*2^(L-M) = N*2^(L-M)，这里N为2的整数次幂，即N=2^M，

W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))

所以，p = J*2^(M-L)，此处J = 0, 1, ..., 2^(L-1)-1，当J遍历结束但计算点数不够N时，J=J+2^L，后继续遍历，直到计算点数为N时不再循环。

流程图：

/*======================================================================
*方法名：fft
*方法功能：计算数组的FFT，运用蝶形运算
*
*变量名称：
*yVector-原始数据
*length-原始数据长度
*N-FFT计算点数
*fftYreal-FFT后的实部
*fftYImg-FFT后的虚部
*
*返回值：是否成功的标志，若成功返回true，否则返回false
*=====================================================================*/

+(BOOL)fft:(floatfloat*)yVectorandOriginalLength:(NSInteger)lengthandFFTCount:(NSInteger)NandFFTReal:(floatfloat*)fftYRealandFFTYImg:(floatfloat*)fftYImg
{
//确保计算时时2的整数幂点数计算
NSIntegerN1=[selfnextNumOfPow2:N];

//定义FFT运算是否成功的标志
BOOLisFFTOK=false;

//判断计算点数是否为2的整数次幂
if(N!=N1)
{
//不是2的整数次幂，直接计算DFT
isFFTOK=[selfdft:yVectorandOriginalLength:lengthandFFTCount:NandFFTReal:fftYRealandFFTYImg:fftYImg];
轿誉
闭肢段//返回成功标志
returnisFFTOK;
}


//如果计算点数位2的整数次幂，用FFT计算，如下
//定义变量
floatyVectorN[N1];//N点运算的原始数据
NSIntegerpowOfN=log2(N1);//N=2^powOfN，用于标记最大运算级数(公式中表示为：M)
NSIntegerlevel=1;//运算级数（第几次运算），最大为powOfN，初值为第一级饥模运算(公式中表示为：L)
NSIntegerlineNum;//行号，倒序排列后的蝶形运算行号(公式中表示为：k)
floatinverseOrderY[N1];//yVector倒序x
NSIntegerdistanceLine=1;//行间距，第level级运算每个蝶形的两个节点距离为distanceLine=2^(L-1)(公式中表示为：B)
NSIntegerp;//旋转因子的阶数，旋转因子表示为W(N,p)，p=J*2^(M-L)
NSIntegerJ;//旋转因子的阶数，旋转因子表示为W(2^L,J)，J=0,1,2,...,2^(L-1)-1=distanceLine-1
floatrealTemp,imgTemp,twiddleReal,twiddleImg,twiddleTheta,twiddleTemp=PI_x_2/N1;
NSIntegerN_4=N1/4;

//判断点数是否够FFT运算点数
if(length<=N1)
{
//如果N至少为length，先把yVector全部赋值
for(NSIntegeri=0;i<length;i++)
{
yVectorN[i]=yVector[i];
}

if(length<N1)
{
//如果N>length后面补零
for(NSIntegeri=length;i<N1;i++)
{
yVectorN[i]=0.0;
}
}
}
else
{
//如果N<length截取相应长度的数据进行运算
for(NSIntegeri=0;i<N1;i++)
{
yVectorN[i]=yVector[i];
}
}

//调用倒序方法
[selfinverseOrder:yVectorNandN:N1andInverseOrderVector:inverseOrderY];

//初始值
for(NSIntegeri=0;i<N1;i++)
{
fftYReal[i]=inverseOrderY[i];
fftYImg[i]=0.0;
}

//三层循环
//第三层（最里）：完成相同旋转因子的蝶形运算
//第二层（中间）：完成旋转因子的变化，步进为2^level
//第一层（最外）：完成M次迭代过程，即计算出x(k)=A0(k),A1(k),...,Am(k)=X(k)

//第一层循环
while(level<=powOfN)
{
distanceLine=powf(2,level-1);//初始条件distanceLine=2^(level-1)
J=0;
NSIntegerpow2_Level=distanceLine*2;//2^level
NSIntegerpow2_NSubL=N1/pow2_Level;//2^(M-L)

//第二层循环
while(J<distanceLine)
{
p=J*pow2_NSubL;
lineNum=J;
NSIntegerstepCount=0;//J运算的步进计数

//求旋转因子
if(p==0)
{
twiddleReal=1.0;
twiddleImg=0.0;
}
elseif(p==N_4)
{
twiddleReal=0.0;
twiddleImg=-1.0;
}
else
{
//计算尤拉公式中的θ
twiddleTheta=twiddleTemp*p;

//计算复数的实部与虚部
twiddleReal=cos(twiddleTheta);
twiddleImg=-11*sin(twiddleTheta);
}

//第三层循环
while(lineNum<N1)
{
//计算下角标
NSIntegerfootNum=lineNum+distanceLine;

/*---------------------------------------
*用复数运算计算每级中各行的蝶形运算结果
*X(k)=X(k)+X(k+B)*W(N,p)
*X(k+B)=X(k)-X(k+B)*W(N,p)
*---------------------------------------*/
realTemp=fftYReal[footNum]*twiddleReal-fftYImg[footNum]*twiddleImg;
imgTemp=fftYReal[footNum]*twiddleImg+fftYImg[footNum]*twiddleReal;

//将计算后的实部和虚部分别存放在返回数组中
fftYReal[footNum]=fftYReal[lineNum]-realTemp;
fftYImg[footNum]=fftYImg[lineNum]-imgTemp;
fftYReal[lineNum]=fftYReal[lineNum]+realTemp;
fftYImg[lineNum]=fftYImg[lineNum]+imgTemp;

stepCount+=pow2_Level;

//行号改变
lineNum=J+stepCount;
}

//旋转因子的阶数变换，达到旋转因子改变的效果
J++;
}

//运算级数加一
level++;
}

isFFTOK=true;
returnisFFTOK;
}

⑻ 基-2fft算法的软件实现 matlab代码

% 基于Matlab的时间抽取基2FFT算法
function y=myditfft(x)
%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次
%------------------------------------
% Leo's fft program(改编网上的一个程序)
%------------------------------------
m=log2(2^nextpow2(length(x))); %求的x长度对应的2的最低幂次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若长度不是2的幂，补0到2的整数幂
end
x;
%--------------------------------------------------------------------------
%对输入序列进行倒序
%如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示
%则其倒位序J对应的二进制数就是（n0n1n2），这样，在原来自然顺序时应该放x(I)的
%单元，现在倒位序后应放x(J)。
%--------------------------------------------------------------------------
%以下程序相当于以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m数列的倒序
%y=x(nxd); %将倒序排列作为初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;

while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
%--------------------------------------------------------------------------
%以下程序解释:
%第一级从x(0)开始,跨接一阶蝶形,再取每条对称
%第二级从x(0)开始,跨接两阶蝶形,再取每条对称
%第m级从x(0)开始,跨接2^(m-1)阶蝶形,再取每条对称....
%--------------------------------------------------------------------------
for mm=1:m %将DFT做m次基2分解，从左到右，对每次分解作DFT运算
Nmr=2^mm;
u=1; %旋转因子u初始化
WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算
for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mm
kp=k+Nmr/2; %确定碟形运算的对应单元下标(对称性)
t=x(kp)*u; %碟形运算的乘积项
x(kp)=x(k)-t; %碟形运算的加法项
x(k)=x(k)+t;
end
u=u*WN; %修改旋转因子，多乘一个基本DFT因子WN
end
end
y=x; %输出

⑼ 求用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;
}

⑽ fft窄带高分辨率算法

fft算法是什么
FFT算法（fast Fourier transform），即快速傅里叶变换，是指利用计算机计算离散傅里叶变换（DFT）的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显着。



概念
有限长睁卖序列可以通过离散傅里叶变换（DFT）将其频域也离散化成有限长序列。但其计算量太大，很难实时地处理问题，因此引出了快速傅里叶变换（FFT）。 1965年，Cooley和Tukey提出了计算离散傅里叶变换（DFT）的快速算法，将DFT的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基2DIT和基2DIF。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的优越性。

如何提高fft算法分辨率
FFT程序，输入是一组复数，输出也是一组复数，想问一下输入到底应该输入什么，输出的复数的含义是什么。

给定一组序列的抽样值，如何用FFT确定它的频率。

首先，fft函数出来的应该是个复数，每一个点分实部虚部两部分。

假设采用1024点fft，悉败逗采样频率是fs，那么第一个点对应0频率点，第512点对应的就是fs/2的频率点。然后从头开始找模值最大的那个点，其所对应的频率值应该就是要的基波频率了。

FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。

有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。

这就是很多信号分析采用FFT变换的原因。

另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。

虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什么意思、如何决定要使用多少点来做FFT。

一个模拟信号，经过ADC采样之后，就变成了数字信号。

采样定理告诉，采样频率要大于信号频率的两倍，这些就不在此罗嗦了。

采样得到的数字信号，就可以做FFT变换了。

N个采样点，经过FFT之后，就可以得到N个点的FFT结果。

为了方便进行FFT运算，通常N取2的整数次方。