您的位置:

fftc语言源代码,fft的c语言代码

本文目录一览:

用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;ln2;l++)

{

if(lj)

{

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 (kj)

{

j=j-k;

k=k/2;

}

j=j+k;

}

for(i=1;i=s;i++)

{

u1=1;

u2=0;

m=(1i);

k=m1;

w1=fcos[i-1];

w2=-fsin[i-1];

for(j=1;j=k;j++)

{

for(l=j;lnn;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]); /* 幅值 */

}

}

(;si=2)

打字不易,如满意,望采纳。

matlab fft2的c代码

傅立叶变换的c语言源代码

128点DIT FFT函数:

/* 采样来的数据放在dataR[ ]数组中,运算前dataI[ ]数组初始化为0 */

void FFT(float dataR[],float dataI[])

{int x0,x1,x2,x3,x4,x5,x6;

int L,j,k,b,p;

float TR,TI,temp;

/********** following code invert sequence ************/

for(i=0;i128;i++)

{ x0=x1=x2=x3=x4=x5=x6=0;

x0=i0x01; x1=(i/2)0x01; x2=(i/4)0x01; x3=(i/8)0x01;x4=(i/16)0x01; x5=(i/32)0x01; x6=(i/64)0x01;

xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;

dataI[xx]=dataR[i];

}

for(i=0;i128;i++)

{ dataR[i]=dataI[i]; dataI[i]=0; }

/************** following code FFT *******************/

for(L=1;L=7;L++) { /* for(1) */

b=1; i=L-1;

while(i0)

{b=b*2; i--;} /* b= 2^(L-1) */

for(j=0;j=b-1;j++) /* for (2) */

{ p=1; i=7-L;

while(i0) /* p=pow(2,7-L)*j; */

{p=p*2; i--;}

p=p*j;

for(k=j;k128;k=k+2*b) /* for (3) */

{ TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];

dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];

dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];

dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];

dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];

} /* END for (3) */

} /* END for (2) */

} /* END for (1) */

for(i=0;i32;i++){ /* 只需要32次以下的谐波进行分析 */

w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);

w[i]=w[i]/64;}

w[0]=w[0]/2;

} /* END FFT */

基于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; iN; 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; qp; q++)

{

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

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

}

p *= 2;

}

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

for(i=0; iN; 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; iN; 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; jb; 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; kN; 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; igr_num; i++)

{

/*

* 对样本点进行循环, 注意上限和步进

*/

for(p=0; psample_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; iN; 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博客,转载请标明出处:

怎样用C语言实现FFT算法啊?

1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下:

先对各行逐一进行一维FFT,然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示:

for (int i=0; iM; i++)

FFT_1D(ROW[i],N);

for (int j=0; jN; j++)

FFT_1D(COL[j],M);

其中,ROW[i]表示矩阵的第i行。注意这只是一个简单的记法,并不能完全照抄。还需要通过一些语句来生成各行的数据。同理,COL[i]是对矩阵的第i列的一种简单表示方法。

所以,关键是一维FFT算法的实现。

2、例程:

#include stdio.h

#include math.h

#include stdlib.h

#define N 1000

/*定义复数类型*/

typedef struct{

double real;

double img;

}complex;

complex x[N], *W; /*输入序列,变换核*/

int size_x=0;      /*输入序列的大小,在本程序中仅限2的次幂*/

double PI;         /*圆周率*/

void fft();     /*快速傅里叶变换*/

void initW();   /*初始化变换核*/

void change(); /*变址*/

void add(complex ,complex ,complex *); /*复数加法*/

void mul(complex ,complex ,complex *); /*复数乘法*/

void sub(complex ,complex ,complex *); /*复数减法*/

void output();

int main(){

int i;                             /*输出结果*/

system("cls");

PI=atan(1)*4;

printf("Please input the size of x:\n");

scanf("%d",size_x);

printf("Please input the data in x[N]:\n");

for(i=0;isize_x;i++)

   scanf("%lf%lf",x[i].real,x[i].img);

initW();

fft();

output();

return 0;

}

/*快速傅里叶变换*/

void fft(){

int i=0,j=0,k=0,l=0;

complex up,down,product;

change();

for(i=0;i log(size_x)/log(2) ;i++){   /*一级蝶形运算*/

   l=1i;

   for(j=0;jsize_x;j+= 2*l ){             /*一组蝶形运算*/

    for(k=0;kl;k++){        /*一个蝶形运算*/

      mul(x[j+k+l],W[size_x*k/2/l],product);

      add(x[j+k],product,up);

      sub(x[j+k],product,down);

      x[j+k]=up;

      x[j+k+l]=down;

    }

   }

}

}

/*初始化变换核*/

void initW(){

int i;

W=(complex *)malloc(sizeof(complex) * size_x);

for(i=0;isize_x;i++){

   W[i].real=cos(2*PI/size_x*i);

   W[i].img=-1*sin(2*PI/size_x*i);

}

}

/*变址计算,将x(n)码位倒置*/

void change(){

complex temp;

unsigned short i=0,j=0,k=0;

double t;

for(i=0;isize_x;i++){

   k=i;j=0;

   t=(log(size_x)/log(2));

   while( (t--)0 ){

    j=j1;

    j|=(k  1);

    k=k1;

   }

   if(ji){

    temp=x[i];

    x[i]=x[j];

    x[j]=temp;

   }

}

}

/*输出傅里叶变换的结果*/

void output(){

int i;

printf("The result are as follows\n");

for(i=0;isize_x;i++){

   printf("%.4f",x[i].real);

   if(x[i].img=0.0001)printf("+%.4fj\n",x[i].img);

   else if(fabs(x[i].img)0.0001)printf("\n");

   else printf("%.4fj\n",x[i].img);

}

}

void add(complex a,complex b,complex *c){

c-real=a.real+b.real;

c-img=a.img+b.img;

}

void mul(complex a,complex b,complex *c){

c-real=a.real*b.real - a.img*b.img;

c-img=a.real*b.img + a.img*b.real;

}

void sub(complex a,complex b,complex *c){

c-real=a.real-b.real;

c-img=a.img-b.img;

}

求FFT的C语言程序……最好是1024点的……希望大家帮帮我!

float ar[1024],ai[1024];/* 实部,虚部 */

float a[2050]; /* 实际值 */

void fft()

{

int n1,n2,i,j,k,l,m,s=10,nn=1024,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,};

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;ln2;l++)

{

if(lj)

{

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 (kj)

{

j=j-k;

k=k/2;

}

j=j+k;

}

for(i=1;i=s;i++)

{

u1=1;

u2=0;

m=(1i);

k=m1;

w1=fcos[i-1];

w2=-fsin[i-1];

for(j=1;j=k;j++)

{

for(l=j;lnn;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]=a[2*i+2]/nn;

ai[i]=-a[2*i+3]/nn;

a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]);

}

}

C语言 1024点快速傅里叶变换(FFT)程序,最好经过优化,执行速度快

void fft()

{

int nn,n1,n2,i,j,k,l,m,s,l1;

float ar[1024],ai[1024]; // 实部 虚部

float a[2050];

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,};

nn=1024;

s=10;

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;ln2;l++)

{

if(lj)

{

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 (kj)

{

j=j-k;

k=k/2;

}

j=j+k;

}

for(i=1;i=s;i++)

{

u1=1;

u2=0;

m=(1i);

k=m1;

w1=fcos[i-1];

w2=-fsin[i-1];

for(j=1;j=k;j++)

{

for(l=j;lnn;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]=a[2*i+2]/nn;

ai[i]=-a[2*i+3]/nn;

a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); // 幅值

}

}