FFT程序說明

小弟近日在网上看见一个FFT（快速傅立叶变换）的C源代码，正好需要，可是看完代码后还是不是很明白。希望各大侠给与帮祝，能否给该代码加上注解？谢谢！

PS：这个程序的输入数据是什么,输出数据是什么？


代码如下：

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include <stdio.h>

void four1(double data[65], 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;
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];
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 twofft(double data1[], double data2[], double fft1[], double fft2[], int& n)
{
int j,n2,j2;
double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i;
c1r = 0.5;
c1i = 0.0;
c2r = 0.0;
c2i = -0.5;
for (j = 1; j<=n; j++)
{
fft1[2 * j - 1] = data1[j];
fft1[2 * j] = data2[j];
}
four1(fft1, n, 1);
fft2[1] = fft1[2];
fft2[2] = 0.0;
fft1[2] = 0.0;
n2 = 2 * (n + 2);
for (j = 2; j<=n / 2 + 1; j++)
{
j2 = 2 * j;
conjr = fft1[n2 - j2 - 1];
conji = -fft1[n2 - j2];
h1r=c1r*(fft1[j2-1]+conjr)-c1i*(fft1[j2]+conji);
h1i=c1i*(fft1[j2-1]+conjr)+c1r*(fft1[j2]+conji);
h2r=c2r*(fft1[j2-1]-conjr)-c2i*(fft1[j2]-conji);
h2i=c2i*(fft1[j2-1]-conjr)+c2r*(fft1[j2]-conji);
fft1[j2 - 1] = h1r;
fft1[j2] = h1i;
fft1[n2 - j2 - 1] = h1r;
fft1[n2 - j2] = -h1i;
fft2[j2 - 1] = h2r;
fft2[j2] = h2i;
fft2[n2 - j2 - 1] = h2r;
fft2[n2 - j2] = -h2i;
}
}


int cint(double x)
{
int temp;
double iprt;
if (x>0)
{
x=modf(x,&iprt);
if(fabs(x)<0.5)
temp=int(iprt);
else
temp=int(iprt+1);
}
else if(x==0)
temp=0;
else
{
x=modf(x,&iprt);
if(fabs(x)<0.5)
temp=int(iprt);
else
temp=int(iprt)-1;
}
return temp;
}

void prntft(double data[], double nn2)
{
int n,m,mm;
cout<<"n real(n) imag.(n) real(N-n) imag.(N-n)"<<endl;
cout<<setw(1)<<"0";
cout<<setw(14)<<data[1];
cout<<setw(14)<<data[2];
cout<<setw(14)<<data[1];
cout<<setw(14)<<data[2]<<endl;
for (n = 3; n<=(nn2 / 2) + 1; n=n+2)
{
m = (n - 1) / 2;
mm = nn2 + 2 - n;
cout<<setiosflags(ios::fixed);
cout<<setprecision(0)<<setw(1)<<m;
cout<<setprecision(6)<<setw(14)<<data[n];
cout<<setprecision(6)<<setw(14)<<data[n+1];
cout<<setprecision(6)<<setw(14)<<data[mm];
cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl;
}
}

void main()
{
//program d12r2
//driver for routine twofft
int n,i,n2,isign;
double data1[33], data2[33], fft1[65], fft2[65],per,x;
n = 32;
n2 = 2 * n;
per = 8.0;
const double pi = 3.1415926;
for(i = 1; i<=n; i++)
{
x = 2.0 * pi * i / per;
data1[i] = cint(cos(x));
data2[i] = cint(sin(x));
}
twofft(data1, data2, fft1, fft2, n);

cout<<setw(1)<<"Fourier transform of first function:"<<endl;
prntft(fft1, n2);
cout<<setw(1)<<"Fourier transform of second function:"<<endl;
prntft(fft2, n2);
//invert transform
isign = -1;
four1(fft1, n, isign);
cout<<setw(1)<<"Inverted transform = first function:"<<endl;
prntft(fft1, n2);
four1(fft2, n, isign);
cout<<setw(1)<<"Inverted transform = second function:"<<endl;
prntft(fft2, n2);
}