#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define PI 3.141592654
#define TWOPI (2.0*PI)
struct COMPLEX
{
float re;
float im;
} cplx , * Hfield,*A;
void four1(int nn, int isign);
void initiate ();
int loop (int l);
//void flybh(double * y1,cplx * y2);
double *data;
int m,n;
void main()
{
int i,j;
int a,b,c;
m=256;n=256;
A = (struct COMPLEX *) malloc ((n*m+1)* sizeof (cplx));
for(i=0;i<=m-1;i++)
{
for(j=1;j<=n;j++)
{
A[i*n+j].re=(float)(1);
A[i*n+j].im=0;
// printf("<%.1f %.1f>",A[i*n+j].re,A[i*n+j].im);
}
// printf("\n");
}
initiate ();//判断m,n;对m,n填充
a=clock();
//对每一行进行fft
data = (double *) malloc ((n*2+1) * sizeof (double));
for(i=0;i<=m-1;i++)
{
for(j=1;j<=n;j++)
{
data[2*j-1]=A[i*n+j].re;
data[2*j]=A[i*n+j].im;
}
four1(n, 1);
for(j=1;j<=n;j++)
{
A[i*n+j].re=data[2*j-1];
A[i*n+j].im=data[2*j];
}
}
//对每一列进行fft
data = (double *) malloc ((m*2+1) * sizeof (double));
for(j=1;j<=n;j++)
{
for(i=0;i<=m-1;i++)
{
data[2*i+1]=A[i*n+j].re;
data[2*i+2]=A[i*n+j].im;
}
four1(m, 1);
for(i=0;i<=m-1;i++)
{
A[i*n+j].re=data[2*i+1];
A[i*n+j].im=data[2*i+2];
}
}
b=clock();
c=b-a;
printf("%d毫秒\n",c);
// a=clock();
// four1(n, 1);
// b=clock();
// c=b-a;
// printf("%d毫秒\n",c);
// for(i=1;i<=256*2;i+=2)
// {
// data[i]=data[i]/(256);
// data[i+1]=data[i+1]/(256);
//
// }
// four1(256, -1);
for(i=0;i<=m-1;i++)
{
for(j=1;j<=n;j++)
{
// printf("<%.2f %.2f>",A[i*n+j].re,A[i*n+j].im);
}
// printf("\n");
}
}
void four1(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 initiate ()
{
int ln,lm;
if ((ln = loop (n)) == -1)
{
printf (" 列数不是2的整数次幂 ");
exit (1);
}
if ((lm = loop (m)) == -1)
{
printf (" 行数不是2的整数次幂 ");
exit (1);
}
// Hfield = (struct COMPLEX *) malloc (n * m * sizeof (cplx));
// data = (double *) malloc ((n*2+1) * sizeof (double));
// for(i=0;i<=m-1;i++)
// {
// for(j=0;j<=n-1;j++)
// {
// Hfield[i*n+j].im=0;
// Hfield[i*n+j].re=i*n+j+1;
// }
// }
// for(i=1;i<=n*2-1;i=i+2)
// {
// data[i]=Hfield[(i+1)/2-1].re;
// data[i+1]=Hfield[(i+1)/2-1].im;
// }
}
int loop (int l)
{//检验输入数据是否为2的整数次幂，如果是返回用2进制表示时的位数
int i , m;
if (l != 0)
{
for (i = 1 ; i < 32 ; i++)
{
m = l >> i;
if (m == 0)
break;
}
if (l == (1 << (i - 1)))
return i - 1;
}
return -1;
}
评论1