求个快速傅里叶变换的C语言程序

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

     }

}

具体实例如下：

1.对于时间序列，可以展开成傅立叶级数，进行频谱分析。对于时间序列xt其傅立叶级数展开式为展开成傅立叶级数：

2.傅立叶分析工具应用 *** 作步骤：

（1）输入数据并中心化：时间、时间序号t、观测值xt、中心化(减x平均值）、求频率fi(=i/N).

（2）由傅立叶分析工具求中心化数据序列的傅立叶变换。

（3）IMREAL和IMAGINARY提取实部和虚部，按公式5计算频率强度（或由IMCONJUGATE求得共轭复数，再由IMPRODUCT求得两共轭复数乘积，得频率强度。

（4）以频率为横坐标、频率强度为纵坐标，绘制频率强度图。

（5）分析周期性。由频率强度最大的所对应的配运频率倒数即得周期。

3.由图可见，序列显现周期性变化，在整个时期范围内，周期为4.下面利用傅立叶分析工具进行频谱分析。

（1）在B18单元格输入“=AVERAGE(B2:B17)”求得观测值的平均值；在C2单元格输入“=B2/B$18”，将观测值中心化（均值为0,并仍保持原序列的方差），并复制到C3:C17。

（2）从“数据”选项卡选择“数培大梁据分析”｜选择“傅利叶分析”d出对话框并设置如

（3）单击“确定”生成傅立叶变换序列（图 20‑2 D列）。

（4）在E2单元格输入“=IMCONJUGATE(D2)”求得傅利叶变换值的共轭复数，并复制到E3:E17；在F3至F17输入1至15,列出周期序列；在G3单元格输入“=F3/16”求得仿行频率，并复制到G4:G17；在H3单元格输入“=IMPRODUCT(D3:E3)*8”（即根据公式5）求得频率强度，并复制到H4:H17。

（5）以G3:H17为源数据，插入散点图，得图 20‑4所示频率强度频谱图。

由图可见，图形完全对称，通常只取左半部分。频率强度最大的所对应的频率为0.25,其倒数为4,即周期为4.

---