求个快速傅里叶变换的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])  // 幅值

     }

}

这是我写的1024点的快速傅里叶变换程序，下面有验证，你把数组

double

A[2049]={0}

double

B[1100]={0}

double

powerA[1025]={0}

改成

A[256]={0}

B[130]={0}

power[129]={0}就行了，

void

FFT(double

data[],

int

nn,

int

isign)

的程序可以针对任何点数，只要是2的n次方

具体程序如下：

#include

<iostream.h>

#include

"math.h"

#include<stdio.h>

#include<string.h>

#include

<stdlib.h>

#include

<fstream.h>

#include

<afx.h>

void

FFT(double

data[],

int

nn,

int

isign)

{

//复数的快速傅里叶变换

int

n,j,i,m,mmax,istep

double

tempr,tempi,theta,wpr,wpi,wr,wi,wtemp

n

=

2

*

nn

j

=

1

for

(i

=

1

i<=n

i=i+2)

//这个循环进行的是码位倒置。

{

if(

j

>

i)

{

tempr

=

data[j]

tempi

=

data[j

+

1]

data[j]

=

data[i]

data[j

+

1]

=

data[i

+

1]

data[i]

=

tempr

data[i

+

1]

=

tempi

}

m

=

n

/

2

while

(m

>=

2

&&

j

>

m)

{

j

=

j

-

m

m

=

m

/

2

}

j

=

j

+

m

}

mmax

=

2

while(

n

>

mmax

)

{

istep

=

2

*

mmax

//这里表示一次的数字的变化。也体现了级数，若第一级时，也就是书是的第0级，其为两个虚数，所以对应数组应该增加4，这样就可以进入下一组运算

theta

=

-6.28318530717959

/

(isign

*

mmax)

wpr

=

-2.0

*

sin(0.5

*

theta)*sin(0.5

*

theta)

wpi

=

sin(theta)

wr

=

1.0

wi

=

0.0

for(

m

=

1

m<=mmax

m=m+2)

{

for

(i

=

m

i<=n

i=i+istep)

{

j

=

i

+

mmax

tempr=double(wr)*data[j]-double(wi)*data[j+1]//这两句表示蝶形因子的下一个数乘以W因子所得的实部和虚部。

tempi=double(wr)*data[j+1]+double(wi)*data[j]

data[j]

=

data[i]

-

tempr

//蝶形单元计算后下面单元的实部，下面为虚部，注意其变换之后的数组序号与书上蝶形单元是一致的

data[j

+

1]

=

data[i

+

1]

-

tempi

data[i]

=

data[i]

+

tempr

data[i

+

1]

=

data[i

+

1]

+

tempi

}

wtemp

=

wr

wr

=

wr

*

wpr

-

wi

*

wpi

+

wr

wi

=

wi

*

wpr

+

wtemp

*

wpi

+

wi

}

mmax

=

istep

}

}

void

main()

{

//本程序已经和MATLAB运算结果对比，准确无误，需要注意的的是，计算告扒中数组都是从1开始取得，丢弃了A[0]等数据

double

A[2049]={0}

double

B[1100]={0}

double

powerA[1025]={0}

char

line[50]

char

dataA[20],

dataB[20]

int

ij

char

ch1[3]="巧带\t"

char

ch2[3]="\n"

int

strl1,strl2

CString

str1,str2

ij=1

//********************************读入文件data1024.txt中的数据,

其中的数据格式见该文件

FILE

*fp

=

fopen("data1024.txt","r")

if(!fp)

{

cout<<"Open

file

is

failing!"<<endl

return

}

while(!feof(fp))

//feof(fp)有两个返回值:如果遇到文件结束，函数feof（fp）的值为1，否则为0。

{

memset(line,0,50)

//清空为0

memset(dataA,0,20)

memset(dataB,0,20)

fgets(line,50,fp)

//函数的功能是从fp所指文件中读入n-1个字符放入line为起始地址的空间内

sscanf(line,

"%s%s",

dataA,

dataB)

//我同时读入孝友芦了两列值，但你要求1024个，那么我就只用了第一列的1024个值

//dataA读入第一列，dataB读入第二列

B[ij]=atof(dataA)

//将字符型的dataA值转化为float型

ij++

}

for

(int

mm=1mm<1025mm++)//A[2*mm-1]是实部，A[2*mm]是虚部，当只要输入实数时，那么保证虚部A[mm*2]为零即可

{

A[2*mm-1]=B[mm]

A[2*mm]=0

}

//*******************************************正式计算FFT

FFT(A,1024,1)

//********************************************写入数据到workout.txt文件中

for

(int

k=1k<2049k=k+2)

{

powerA[(k+1)/2]=sqrt(pow(A[k],2.0)+pow(A[k+1],2.0))//求功率谱

FILE

*pFile=fopen("workout.txt","a+")

//?a+只能在文件最后补充，光标在结尾。没有则创建

memset(ch1,0,15)

str1.Format("%.4f",powerA[(k+1)/2])

if

(A[k+1]>=0)

str2.Format("%d\t%6.4f%s%6.4f

%s",(k+1)/2,A[k],"+",A[k+1],"i")//保存fft计算的频谱，是复数频谱

else

str2.Format("%d\t%6.4f%6.4f

%s",(k+1)/2,A[k],A[k+1],"i")

strl1=strlen(str1)

strl2=strlen(str2)

//

用

法:fwrite(buffer,size,count,fp)

//

buffer：是一个指针，对fwrite来说，是要输出数据的地址。

//

size：要写入的字节数；

//

count:要进行写入size字节的数据项的个数；

//

fp:目标文件指针。

fwrite(str2,1,strl2,pFile)

fwrite(ch1,1,3,pFile)

fwrite(ch1,1,3,pFile)

fwrite(str1,1,strl1,pFile)

fwrite(ch2,1,3,pFile)

fclose(pFile)

}

cout<<"计算完毕，到fft_test\workout.txt查看结果"<<endl

}

---