图片显示

开发平台：
Visual C++

fft.cpp：源码内容
							//#include "stdafx.h"
//#include "complex"
#include <math.h>
//using namespace std;
#define PI 3.1415926535
void fft(double y[][2], double * pCos, double *pSin, int * l, int N) //N is the power of 2 (as 8,64, 1024)
{
	double t1[2], t2[2], t[2];
	int u,mu,i,j,k,/*l,*/s;
	mu=0; u=N-1;
	while (u>0) {mu++; u/=2;}
	for (i=0; i<N; i++) { //sort input number array
	  if (l[i]>i) {
			t[0] = y[i][0]; t[1] = y[i][1];
			y[i][0] = y[l[i]][0]; y[i][1] = y[l[i]][1];
			y[l[i]][0] = t[0]; y[l[i]][1] = t[1];
	  }
	}
	s=1;
	for (i=0; i<mu; i++) { //fft by binary dividing
	  for (j=0; j<N; j+=s*2) {
		for (k=0; k<s; k++) {
			  t[0] = y[j+k+s][0]*pCos[k*N/s/2] - y[j+k+s][1]*pSin[k*N/s/2];
			  t[1] = y[j+k+s][0]*pSin[k*N/s/2] + y[j+k+s][1]*pCos[k*N/s/2];
			  t1[0] = y[j+k][0] + t[0]; t1[1] = y[j+k][1] + t[1];
			  t2[0] = y[j+k][0] - t[0]; t2[1] = y[j+k][1] - t[1];
			  y[j+k][0]=t1[0]; y[j+k][1]=t1[1];
			  y[j+k+s][0]=t2[0]; y[j+k+s][1]=t2[1];
		}
	  }
	  s*=2;
	}
}
void calW(double * pCos, double * pSin, int * l, int N)
{
  int i,j,k;
  for(i=0; i<N/2; i++)
  {
	  double theta = 2*PI/N*i;
	  pCos[i] = cos(theta);
	  pSin[i] = sin(theta);
  }
  for (i=0; i<N; i++) 
  { //sort input number array
    k=i; l[i]=0;
    for (j=2; j<=N; j=j*2) 
	{
      l[i] += (N/j)*(k%2); k=k/2;
    }
  }
}
void imposition(double y[][2], long width)
{
	long i,j,pos1,pos2;
	double temp[2];
	long height = width;
	for(j = 0; j < height; j ++)
	{
		for(i = 0; i < j; i ++)
		{
			pos1 = j*height + i;
			pos2 = i*width + j;
			temp[0] = y[pos1][0];
			temp[1] = y[pos1][1];
			y[pos1][0] = y[pos2][0];
			y[pos1][1] = y[pos2][1];
			y[pos2][0] = temp[0];
			y[pos2][1] = temp[1];
		}
	}
	return;
}
void imposition(double y[], long width)
{
	long i,j,pos1,pos2;
	double temp;
	long height = width;
	for(j = 0; j < height; j ++)
	{
		for(i = 0; i < j; i ++)
		{
			pos1 = j*height + i;
			pos2 = i*width + j;
			temp = y[pos1];
			y[pos1] = y[pos2];
			y[pos2] = temp;
		}
	}
	return;
}
void dct(double x[], double * pCos, double * pSin, int * l, int N) //N is the power of 2 (as 8,64, 1024)
{
	int i;
	double (* y)[2] = new double[N][2];
	double c, s, a;
	for(i = 0; i < N/2; i ++)
	{
		y[i][0] = x[2*i];
		y[i][1] = 0.0;
		y[N-1-i][0] = x[2*i+1];
		y[N-1-i][1] = 0.0;
	}
	fft(y, pCos, pSin, l, N);
	a = sqrt(2.0/(double)N);
	for(i = 0; i < N; i ++)
	{
		c = cos(i*PI/(double)(2*N));
		s = sin(i*PI/(double)(2*N));
		x[i] = (c*y[i][0] - s*y[i][1]) * a;
	}
	x[0] = x[0]/sqrt(2.0);
	if(y) delete[] y;
}
void fft(double pr[][2],int n,int k,
		        double f[][2],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;
		}
        f[it][0]=pr[is][0]; f[it][0]=pr[is][0];
    }
    pr[0][0]=1.0; pr[0][1]=0.0;
    p=6.283185306/(1.0*n);
    pr[1][0]=cos(p); pr[1][1]=-sin(p);
    if (l!=0)
	{
		pr[1][0]=-pr[1][0];
	}
    for (i=2; i<=n-1; i++)
    { 
		p=pr[i-1][0]*pr[1][0]; q=pr[i-1][1]*pr[1][1];
        s=(pr[i-1][0]+pr[i-1][1])*(pr[1][0]+pr[1][1]);
        pr[i][0]=p-q; pr[i][1]=s-p-q;
    }
    for (it=0; it<=n-2; it=it+2)
    { 
		vr=f[it][0]; vi=f[it][1];
        f[it][0]=vr+f[it+1][0]; f[it][1]=vi+f[it+1][1];
        f[it+1][0]=vr-f[it+1][0]; f[it+1][1]=vi-f[it+1][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][0]*f[it+j+nv/2][0];
              q=pr[m*j][1]*f[it+j+nv/2][1];
              s=pr[m*j][0]+pr[m*j][1];
              s=s*(f[it+j+nv/2][0]+f[it+j+nv/2][1]);
              poddr=p-q; poddi=s-p-q;
              f[it+j+nv/2][0]=f[it+j][0]-poddr;
              f[it+j+nv/2][1]=f[it+j][1]-poddi;
              f[it+j][0]=f[it+j][0]+poddr;
              f[it+j][1]=f[it+j][1]+poddi;
            }
		}
     }
    if (l!=0)
    {  for (i=0; i<=n-1; i++)
        { f[i][0]=f[i][0]/(1.0*n);
          f[i][1]=f[i][1]/(1.0*n);
        }
	}
    if (il!=0)
	{
		for (i=0; i<=n-1; i++)
        { pr[i][0]=sqrt(f[i][0]*f[i][0]+f[i][1]*f[i][1]);
          if (fabs(f[i][0])<0.000001*fabs(f[i][1]))
          { 
			  if ((f[i][1]*f[i][0])>0) 
				  pr[i][1]=90.0;
              else
				  pr[i][1]=-90.0;
          }
          else
            pr[i][1]=atan(f[i][1]/f[i][0])*360.0/6.283185306;
        }
	}
}

						