傅里叶变换用C语言程序怎么实现？

#include <math.h>

#include <stdio.h>

#define N 8

void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)

void main()

{

double xr[N],xi[N],Yr[N],Yi[N],l=0,il=0

int i,j,n=N,k=3

for(i=0i<Ni++)

{

xr[i]=i

xi[i]=0

}

printf("------FFT------\n")

l=0

kkfft(xr,xi,n,k,Yr,Yi,l,il)

for(i=0i<Ni++)

{

printf("%-11lf + j* %-11lf\n",Yr[i],Yi[i])

}

printf("-----DFFT-------\n")

l=1

kkfft(Yr,Yi,n,k,xr,xi,l,il)

for(i=0i<Ni++)

{

printf("%-11lf + j* %-11lf\n",xr[i],xi[i])

}

getch()

}

void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)

{

int it,m,is,i,j,nv,l0

double p,q,s,vr,vi,poddr,poddi

for (it=0it<=n-1it++)

{

m = it

is = 0

for(i=0i<=k-1i++)

{

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=2i<=n-1i++)

{

p = pr[i-1]*pr[1]

q = pi[i-1]*pi[1]

s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1])

pr[i] = p-q

pi[i] = s-p-q

}

for (it=0it<=n-2it=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-2l0>=0l0--)

{

m = m/2

nv = 2*nv

for(it=0it<=(m-1)*nvit=it+nv)

for (j=0j<=(nv/2)-1j++)

{

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=0i<=n-1i++)

{

fr[i] = fr[i]/(1.0*n)

fi[i] = fi[i]/(1.0*n)

}

}

/*是否计算模和相角*/

if(il!=0)

{

for(i=0i<=n-1i++)

{

pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i])

if(fabs(fr[i])<0.000001*fabs(fi[i]))

{

if ((fi[i]*fr[i])>0)

pi[i] = 90.0

else

pi[i] = -90.0

}

else

pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306

}

}

return

}

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=1i<=nni++)

      {

        a[2*i]=ar[i-1]

        a[2*i+1]=ai[i-1]

      }

      for(l=1l<n2l++)

      {

       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=1i<=si++)

     {

    u1=1

    u2=0

    m=(1<<i)

    k=m>>1

    w1=fcos[i-1]

    w2=-fsin[i-1]

    for(j=1j<=kj++)

    {

     for(l=jl<nnl=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=1i<=nn/2i++)

     {

    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])  // 幅值

     }

}

---