网站首页
这些小活动你都参加了吗？快来围观一下吧！>>
电子产品世界 » 论坛首页 » 活动中心 » 合作大赛 » 傅立叶算法的Ｃ语言实现

共1条 1/1 1 跳转至页

傅立叶算法的Ｃ语言实现

diansaizhuanti
高工
2015-06-11 11:09:18     打赏
/***** 作者  _yuming           *************/ #include "math.h" 
  void kfft(pr,pi,n,k,fr,fi,l,il) 
  int n,k,l,il; 
  double pr[],pi[],fr[],fi[]; 
  { int it,m,is,i,j,nv,l0; 
    double p,q,s,vr,vi,poddr,poddi; 
    for (it=0; it<=n-1; it++) 
      { m=it; is=0; 
        for (i=0; i<=k-1; i++) 
          { j=m/2; is=2*is+(m-2*j); m=j;} 
        fr[it]=pr[is]; fi[it]=pi[is]; 
      } 
    pr[0]=1.0; pi[0]=0.0; 
    p=6.283185306/(1.0*n); 
    pr[1]=cos(p); pi[1]=-sin(p); 
    if (l!=0) pi[1]=-pi[1]; 
    for (i=2; i<=n-1; i++) 
      { p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1]; 
        s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]); 
        pr=p-q; pi=s-p-q; 
      } 
    for (it=0; it<=n-2; it=it+2) 
      { vr=fr[it]; vi=fi[it]; 
        fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1]; 
        fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1]; 
      } 
    m=n/2; nv=2; 
    for (l0=k-2; l0>=0; l0--) 
      { m=m/2; nv=2*nv; 
        for (it=0; it<=(m-1)*nv; it=it+nv) 
          for (j=0; j<=(nv/2)-1; j++) 
            { p=pr[m*j]*fr[it+j+nv/2]; 
              q=pi[m*j]*fi[it+j+nv/2]; 
              s=pr[m*j]+pi[m*j]; 
              s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]); 
              poddr=p-q; poddi=s-p-q; 
              fr[it+j+nv/2]=fr[it+j]-poddr; 
              fi[it+j+nv/2]=fi[it+j]-poddi; 
              fr[it+j]=fr[it+j]+poddr; 
              fi[it+j]=fi[it+j]+poddi; 
            } 
      } 
    if (l!=0) 
      for (i=0; i<=n-1; i++) 
        { fr=fr/(1.0*n); 
          fi=fi/(1.0*n); 
        } 
    if (il!=0) 
      for (i=0; i<=n-1; i++) 
        { pr=sqrt(fr*fr+fi*fi); 
          if (fabs(fr)<0.000001*fabs(fi)) 
            { if ((fi*fr)>0) pi=90.0; 
              else pi=-90.0; 
            } 
          else 
            pi=atan(fi/fr)*360.0/6.283185306; 
        } 
    return; 
  }
回复
收藏
0


共1条 1/1 1 跳转至页

回复

匿名不能发帖！请先 [ 登陆 注册 ]