求傅里叶逆变换的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=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][/it]it] = pr[is]

       fi[[it][/it]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[i] = p-q 

      pi[i] = s-p-q

    }

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

    { 

      vr = fr[[it][/it]it] 

       vi = fi[[it][/it]it]

       fr[[it][/it]it] = vr+fr[it+1] 

       fi[[it][/it]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[i] = fr[i]/(1.0*n)

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

       }

    } 

    

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

    if(il!=0)

    {

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

       { 

        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

}

1.FFT：

// data为输入和输出的数据，N为长度

bool CFFT::Forward(complex *const Data, const unsigned int N)

{

   if (!Data || N < 1 || N & (N - 1))

      return false

   //   排序

   Rearrange(Data, N)

   //   FFT计改迟算：const bool Inverse = false

   Perform(Data, N)

   return true

}

 2.IFFT：

// Scale 为是否缩放

bool CFFT::Inverse(complex *const Data, const unsigned int N,

   const bool Scale /* = true */)

{

   if (!Data || N < 1 || N & (N - 1))

      return false

   //   排序

   Rearrange(Data, N)

   //   FFT计算，ture表示是逆运算

   Perform(Data, N, true)

   //   对结果进行缩放

   if (Scale)

      CFFT::Scale(Data, N)

   return true

}

3.排序：

void CFFT::Rearrange(complex *const Data, const unsigned int N)

{

   //   Swap position

   unsigned int Target = 0

   //   Process all positions of input signal

   for (unsigned int Position = 0 Position < N ++Position)

   {

      //   晌扰Only for not yet swapped entries

      if (Target > Position)

      {

         //   Swap entries

         const complex Temp(Data[Target])

         Data[Target] = Data[Position]

         宴歼旦Data[Position] = Temp

      }

      //   Bit mask

      unsigned int Mask = N

      //   While bit is set

      while (Target & (Mask >>= 1))

         //   Drop bit

         Target &= ~Mask

      //   The current bit is 0 - set it

      Target |= Mask

   }

}

 4.FFT计算：

void CFFT::Perform(complex *const Data, const unsigned int N,

   const bool Inverse /* = false */)

{

   const double pi = Inverse ? 3.14159265358979323846 : -3.14159265358979323846

   //   Iteration through dyads, quadruples, octads and so on...

   for (unsigned int Step = 1 Step < N Step <<= 1)

   {

      //   Jump to the next entry of the same transform factor

      const unsigned int Jump = Step << 1

      //   Angle increment

      const double delta = pi / double(Step)

      //   Auxiliary sin(delta / 2)

      const double Sine = sin(delta * .5)

      //   Multiplier for trigonometric recurrence

      const complex Multiplier(-2. * Sine * Sine, sin(delta))

      //   Start value for transform factor, fi = 0

      complex Factor(1.)

      //   Iteration through groups of different transform factor

      for (unsigned int Group = 0 Group < Step ++Group)

      {

         //   Iteration within group 

         for (unsigned int Pair = Group Pair < N Pair += Jump)

         {

            //   Match position

            const unsigned int Match = Pair + Step

            //   Second term of two-point transform

            const complex Product(Factor * Data[Match])

            //   Transform for fi + pi

            Data[Match] = Data[Pair] - Product

            //   Transform for fi

            Data[Pair] += Product

         }

         //   Successive transform factor via trigonometric recurrence

         Factor = Multiplier * Factor + Factor

      }

   }

}

 5.缩放：

void CFFT::Scale(complex *const Data, const unsigned int N)

{

   const double Factor = 1. / double(N)

   //   Scale all data entries

   for (unsigned int Position = 0 Position < N ++Position)

      Data[Position] *= Factor

}

---