寻求快速傅立叶算法的Ｃ语言实现

jack_yu

受教了。谢谢

yinhe

ji

skyzyf

mark

MeSleve

mark

longfeixue

学习

jy6715

学习 了

lxd_youxi

这个好 ，的好好研究研究！

fesu02

mark

xuhaikun

牛啊，记号记号

andriy

标记，估计会用上的

hecb999

FFT MARK

shuli163love

make

xuyqhd

MARK

caiyue3577

不错！学习！

hlswx

MARK，标记一下～

luvemcu

这个运算需要多少个机器周期？

416446891

mark

xph114

mark

century99

mark 各位大侠好厉害！

njust_zgx

Mark!

IGO_AVR

mark

wenming

还可以

qianbo0920

mark

guohua0710

感谢分享！

yfwuh

#include<conio.h>
     #include<graphics.h>
     #include<iostream.h>
     #include<math.h>
     #include<stdio.h>
     #include<stdlib.h>
     #include<string.h>
     #define    PI   3.1415926
//  #define    N    64
     #define    N1   64

struct rptr
{
double *r;
double *i;
}COMPLEX;

void fft(double *dr,double *di,int N);
struct rptr *dftn(double a,double h,int n,double *f,int nw,struct rptr *tran);
void draw(double *x,double *h,int lenx,int lenh);
void main()
{
  int fp=3;
  int i;
  struct rptr *tran;
  double *fw,*f;
  double time=64;//time of signal
  double h=0.25;
  double a=-time/2;
  int N=time/h;
  tran->i=(double *)calloc(N+2,sizeof(double));
  tran->r=(double *)calloc(N+2,sizeof(double));
  fw=(double *)calloc(N+2,sizeof(double));
  f=(double *)calloc(N+2,sizeof(double));

  cout<<"\nPlease select the input function:\n";
  cout<<"1.cos(0.125*PI*t)        2.sin(PI*t)/(PI*t)\n";
  cout<<"3.sin(0.25*PI*t)+2*cos(0.5*PI*t)      4.cos(0.125*PI*t)+2*cos(0.25*PI*t)\n";
  cout<<"5.sin(0.25*PI*t)\n<-----------";
  lp: cin>>fp;
       for(i=0;i<=N;i++)
       {
       double t=(a+h*i);
       switch(fp)
       {
         case 1:   f=tran->r=cos(0.125*PI*t);tran->i=0;
         case 2:   if(t==0)  f=tran->r=1;
                   else
  f=tran->r=sin(1.0*PI*t)/(1.0*PI*t);
             tran->i=0;
             break;
         case 3:
  f=tran->r=sin(0.25*PI*t)+2*cos(0.5*PI*t);tran->i=0;
             break;
         case 4:
  f=tran->r=cos(0.125*PI*t)+2*cos(0.25*PI*t);tran->i=0;
             break;
         case 5:   f=tran->r=sin(0.25*PI*t);tran->i=0;
             break;
         default:   cout<<"Wrong!\n";
                  goto lp;
         }
  }
  fft(tran->r,tran->i,N);
  cout<<"\n\n";
  for(int j=0;j<N;j++)
  {fw[j]=pow(tran->r[j],2.0)+pow(tran->i[j],2.0);
   printf("%-10.2f",fw[j]);
   }getch();
draw(f,fw,N,N);
}



void fft(double *dr,double *di,int N)
{
     int lh,m;
     lh=N/2;
     m=int(log(N)/log(2)+0.9999);
     int j=lh;
     int k;
     for(int i=1;i<(N-1);i++)



{
     if(i<j)
     {
        double tr,ti;
        tr=dr;ti=di;
        dr=dr[j];di=di[j];
        dr[j]=tr;di[j]=ti;
     }
     k=lh;
     while(j>=k) {j=j-k;k=k/2;}
     j=j+k;
}
for(i=1;i<=m;i++)
{
     int b=int(pow(2.0,(i-1)));
     for(j=0;j<b;j++;)
     {
        double p;
        p=j*pow(2.0,(m-i))*2.0*PI/N;
        for(k=j;k<N;)
        {
           double tr,ti;
           tr=dr[k+b]*cos(p)+di[k+b]*sin(p);
           ti=di[k+b]*cos(p)-dr[k+b]*sin(p);
           dr[k+b]=dr[k]-tr;
           di[k+b]=di[k]-ti;
           dr[k]=dr[k]+tr;
           di[k]=di[k]+ti;
           k=int(k+pow(2.0,i));
           }
         }
       }

    }


      void draw(double *x,double *h,int lenx,int lenh)
      {
      int i,j,k;
      double mx,mh,my;
      int gdriver=DETECT,gmode,errorcode;
//    errorcode=registerbgidriver(EGAVGA_driver);

      initgraph(&gdriver,&gmode,"c:\\borlandc\\bgi");
      errorcode=graphresult();

      if(errorcode!=grOk)
      {
          cout<<"Graphics error: %s\n",grapherrormsg(errorcode);
          cout<<"Press any key to halt!\n";
          getch();
          exit(1);
      }
      setbkcolor(WHITE);
      setcolor(2);

      line(80,140,560,140);
      line(80,140,80,40);

      line(80,370,560,370);
      line(80,370,80,270);

//   line(300,250,600,250);
//   line(300,250,300,100);
      double heigh=80;
      long int a;
      int sign=1;
      char tmp[40],tmp1[40];
      mx=x[0];j=0;
      double m=x[0];k=0;
      for(i=1;i<lenx;i++)
      {  if(x<0)   {m=(x<m) ? x:m;if(m==x) k=i;}
         else         {mx=(x>mx) ? x:mx;if(mx==x) j=i;}
      }
      if((-1)*m>mx)
         {
         mx=(-1)*m;
         sign=-1;
         j=k;
         }
      a=mx;
      if(sign==-1)   strcpy(tmp,"-");
      else           strcpy(tmp,"");
      ltoa(a,tmp1,10);
      strcat(tmp,tmp1);
      strcat(tmp,".");
      a=(mx-a)*1000;
      itoa(a,tmp1,10);
      strcat(tmp,tmp1);
      if(sign==-1)   outtextxy(220.0*j/lenx,140+heigh,tmp);

      else           outtextxy(220.0*j/lenx,130-heigh,tmp);
      setcolor(1);
      mx=heigh/mx;
      moveto(80,140-mx*x[0]);
      for(i=0;i<lenx;i++)
         lineto(80+450.0*i/lenx,140-mx*x);


      mh=h[0];j=0;
      m=h[0];k=0;
      sign=1;
      for(i=1;i<lenh;i++)
      {
         if(h<0)   {m=(h<m ? h:m;if(m==h) k=i;}
         else         {mh=(h>mh) ? h:mh;if(mh==h) j=i;}
      }
      if((-1)*m>mh)
         {
         mh=(-1)*m;
         sign=-1;
         j=k;
         }
      a=mh;
      if(sign==-1)   strcpy(tmp,"-");
      else           strcpy(tmp,"");
      ltoa(a,tmp1,10);
      strcat(tmp,tmp1);
      strcat(tmp,".");
      a=(mh-a)*1000;
      itoa(a,tmp1,10);
      strcat(tmp,tmp1);
      setcolor(2);
      if(sign==-1)      outtextxy(220.0*j/lenh,370+heigh,tmp);
      else              outtextxy(220.0*j/lenh,358-heigh,tmp);
      mh=heigh/mh;
      setcolor(1);
      moveto(80,370-mh*h[0]);
      for(i=0;i<lenh;i++)
      {
         moveto(80+450.0*i/lenh,370);
         lineto(80+450.0*i/lenh,370-mh*h);
      }
       settextstyle(DEFAULT_FONT,HORIZ_DIR,2);
       setcolor(7);
       outtextxy(75,15,"Input x(t)");
       outtextxy(70,230,"Response h(t)");
//    outtextxy(320,68,"Onput x(t)*h(t)");
       settextstyle(DEFAULT_FONT,HORIZ_DIR,1);
       setcolor(2);
       outtextxy(44,38,"x(t)");
       outtextxy(44,266,"h(t)");
//    outtextxy(302,96,"x(t)*h(t)");
       outtextxy(272,140,"t");
       outtextxy(272,380,"t");
//    outtextxy(602,252,"t");
       outtextxy(39,150,"0");
       outtextxy(39,374,"0");
//    outtextxy(300,254,"0");
       getch();
       closegraph();
       }

yfwuh

//分 裂 基 FFT 算 法 程 序

/*Free_Copy*/
/*主程序:64点分裂基FFT算法*/
/*输入: 64点任意序列*/
/*输出: 序列的FFT变换*/

#include "conio.h";
#include"math.h"
#include"stdio.h"
#define PI 3.1415926
#define N 128

void main()  
{
   float x[N],y[N],xt;
   float cc1,cc3,ss1,ss3;
   float r1,r2,r3,s1,s2,a,a3,e,m1;
   int n,n1,m,j,k,i;
   int is,id,i0,i1,i2,i3,n2,n4;

   printf("\nThis program is about  FFT by SPEFT way.  ");
   printf("\nplease enter  n : ");
   scanf("%d",&n1);
   n=n1;
   m1=log(n1)/log(2);
   m=log(n1)/log(2);
   if (m!=m1)  n=pow(2,m+1);
   for(i=0;i<=N;i++)
      {
       x=y=0.0;
      }
   printf("\n");
   for(i=1;i<=n1;i++)
      {
       printf("\nplease enter  data(%d)_[Re]: ",i);
       scanf("%f",&x);
       printf("\nplease enter  data(%d)_[Im]: ",i);
       scanf("%f",&y);
      }
   j=1;
   for (i=1;i<=n-1;i++)
      {
       if (i<j)
         {
          xt=x[j];
          x[j]=x;
          x=xt;
          xt=y[j];
          y[j]=y;
          y=xt;
         }
       k=n/2;
       while (k<j)
         {
           j=j-k;
           k=k/2;
         }
        j=j+k;
      }
   is=1;
   id=4;
   while (is<n)
     {
       for (i0=is;i0<=n;i0+=id)
           {
             i1=i0+1;
             r1=x[i0];
             x[i0]=r1+x[i1];
             x[i1]=r1-x[i1];
             r1=y[i0];
             y[i0]=r1+y[i1];
             y[i1]=r1-y[i1];
           }
        is=2*id-1;
        id=4*id;
     }
   n2=2;
   for (k=2;k<=m;k++)
     {
      n2=n2*2;
      n4=n2/4;
      e=2.0*PI/n2;
      a=0.0;
      for (j=1;j<=n4;j++)
          {
           a3=3.0*a;
           cc1=cos(a);
           ss1=sin(a);
           cc3=cos(a3);
           ss3=sin(a3);  
           a=j*e;
           is=j;
           id=2*n2;
           while (is<n)
              {
                for (i0=is;i0<=n-1;i0+=id)
         {
                    i1=i0+n4;
           i2=i1+n4;
           i3=i2+n4;
           r1=x[i2]*cc1+y[i2]*ss1;
           s1=y[i2]*cc1-x[i2]*ss1;
           r2=x[i3]*cc3+y[i3]*ss3;
           s2=y  [i3]*cc3-x[i3]*ss3;
           r3=r1+r2;
           r2=r1-r2;
           r1=s1+s2;
           s2=s1-s2;
           x[i2]=x[i0]-r3;
           x[i0]=x[i0]+r3;
           x[i3]=x[i1]-s2;
           x[i1]=x[i1]+s2;
           y[i2]=y[i0]-r1;
           y[i0]=y[i0]+r1;
           y[i3]=y[i1]+r2;
              y[i1]=y[i1]-r2;
          }
                is=2*id-n2+j;
                id=4*id;
              }
           }
      }
   printf("\n分裂基fft结果是： \n ");
   for (i=1;i<=n;i++)
     {
       printf("\n %7.3f, %7.3fj",x,y);
       y=-y;
     }
   getch();
   printf("\n\n");
}

yfwuh

基- 2 FFT ( 频率抽取DIF 法) 算法程序

/*Free_Copy*/
/*  C语言编写的频率抽取FFT算法(最大计算64点）*/  
/*  输入: 序列点数、序列值  * /
/*  输出: 序列FFT变换后的数值及反变换(应与原序列相同 )  */

#include "conio.h"
#include "math.h"
#include "stdio.h"

#define N 64
#define PI 3.1415926
#define w0 (0.125*PI)
#define Cmul(a,b,c) a.x=b.x*c.x-b.y*c.y;a.y=b.x*c.y+b.y*c.x;
#define   Cequal(a,b) a.x=b.x;a.y=b.y;
#define Cadd(a,b,c) a.x=b.x+c.x;a.y=b.y+c.y;
#define Csub(a,b,c) a.x=b.x-c.x;a.y=b.y-c.y;
#define Wn(w,r) w.x=cos(2.0*PI*r/n);w.y=-sin(2.0*PI*r/n);
struct comp
{
   float x;
   float y;
};
void main()
{
   int i,j,nu2,nm1,n,m,le,le1,k,ip,z;
   int flag,f,n1;
   struct comp a[N],t,t1,w,d;
   float a_ipx,m1;
   printf("\nThis program is about  FFT by DIF way.  ");
   printf("\nplease enter  N : ");
   scanf("%d",&n1);
   n=n1;
   m1=log(n1)/log(2);
   m=log(n1)/log(2);
   if (m!=m1)  n=pow(2,m+1);
   for(i=0;i<n;i++) {a.x=a.y=0.0;}
   printf("\n");
   for(i=0;i<n1;i++)
    {
     printf("\nplease enter  data(%d)_[Re]: ",i);
     scanf("%f",&a.x);
     printf("\nplease enter  data(%d)_[Im]: ",i);
     scanf("%f",&a.y);
    }
   for(z=0;z<=1;z++)
   {
    flag=-1;
    for (m=(log(n)/log(2));m>=1;m--)
      {
      le=pow(2,m);
      flag++;
      le1=le/2;
      for( j=0;j<le1;j++)
         {
          for (i=j;i<=(n-1);i+=le)
             {  
               ip=i+le1;
               Cequal(t,a);
               Cequal(t1,a[ip]);
               f=(int) (i*pow(2,flag))%n;
               Wn(w,f);  
               Cadd(a,t,t1);
               Csub(a[ip],t,t1);
               a_ipx=a[ip].x;
               if (z==1)
                 {  
                  w.y*=-1;
                 }
              a[ip].x=a[ip].x*w.x-a[ip].y*w.y;
              a[ip].y=a_ipx*w.y+a[ip].y*w.x;
             }
          }
      }
    nu2=n/2;
    nm1=n-2;
    j=0;i=0;  
   while(i<=nm1)
      {  
       if (i<j)
         {
          Cequal(d,a[j]);
          Cequal(a[j],a);
          Cequal(a,d);
         }
       k=nu2;
       while(k<=j)
         {
          j=j-k;k=k/2;
         }
        j=j+k;
        i=i+1;
       }
    if(z==0)
       {
        printf("\n序列的fft是:\n\n");
       }
    else
        printf("\n用ifft计算出的原序列是:\n\n"  ) ;
    for(i=0;i<n;i++)
       if(z==0)
          {
           printf("    %7.3f",a.x);
           if (a.y>=0)
           printf("  + %7.3f j \n",a.y);
           else
           printf("  - %7.3f j \n",fabs(a.y));
           a.y= -a.y;
           }
       else
          {   
           printf("   %7.3f",a.x/n);
           a.y=-a.y/n;
           if (a.y>=0)
           printf("  + %7.3f j \n",a.y);
           else
           printf("  - %7.3f j \n",fabs(a.y));
          }
    }
  printf("\n");
}

yfwuh

点击此处下载 ourdev_645194TKANLP.pdf(文件大小:3.41M) (原文件名:快速傅里叶变换及其C程序.pdf)

pujing

mark FFT

leang521

mark

jiabin1024

nb

wenming

支持一下~

huaxin

学习了， 留在后来用的在查，

471176877

mark

starsun

请问，这个算法有什么具体应用吗？

vcrm

mark

safan

mark!留着用的时候参考……

safan

mark!留着用的时候参考……

qq329889612

mark...

wcm_e

mark

ledatou

学习了~~~~

createc

mark

wlzxzp

有空看看

dszx

mark ,有空好好看看

zwc58

mark.

STC89C55

mark

gaomcu

mark

xiongxie007

mark

1988_coolboy

mark

longcai1988

FFT

xqxyh2015

MARK

yangbo_zd

mark fft

ShawnLinson

Mark，FFT最佳例子

wdh2152

mark

pcbddd

mark

CHENXIAOTIAN

mark 学习

eastbest

mark

bai2011

mark

FpvCamera

mark

leiwen_ourdev

mark

lanjingjing

hao

slzxriver

mark

jlhgold

mark

Roshan

mark

semonpic

mark

djl310

收藏~~~~~~~~~~

amote

mark

ss_soldier

强帖留名。
睡前看看。

Ca_guo

不错，FFT应用还是挺多的

gloryzkl

mark

luvemcu

Fast Fourier Transform?

stalker2

想做个ARM对声音信号处理的东西来着 这个帮大忙啦

luckseason

顶顶顶顶。。

myhonour

快速傅里叶，标记一个^_^

guowanling8061

应该是成形的东西了！！！

小孔523

可以先将浮点数放大一定倍数以后变为定点数，再除最后实际的放大倍数，这样可以大大提高运算速度~~

一棵葱

好贴，必须顶！！！

w625419875

标记fft，回头看看

机器人天空

mark.........

知行合一

mark                                       

chinamaken

必须收藏了，做工控的很实用哦

211LIRUISHUO

好东西,收藏!