C++ FFT函数代码示例

// create DCT/DST plan
//  _nfft   :   FFT size
//  _x      :   input array [size: _nfft x 1]
//  _y      :   output array [size: _nfft x 1]
//  _type   :   type (e.g. LIQUID_FFT_REDFT00)
//  _method :   fft method
FFT(plan) FFT(_create_plan_r2r_1d)(unsigned int _nfft,
                                   T *          _x,
                                   T *          _y,
                                   int          _type,
                                   int          _flags)
{
    // allocate plan and initialize all internal arrays to NULL
    FFT(plan) q = (FFT(plan)) malloc(sizeof(struct FFT(plan_s)));

    q->nfft   = _nfft;
    q->xr     = _x;
    q->yr     = _y;
    q->type   = _type;
    q->flags  = _flags;

    // TODO : use separate 'method' for real-to-real types
    //q->method = LIQUID_FFT_METHOD_NONE;

    switch (q->type) {
    case LIQUID_FFT_REDFT00:  q->execute = &FFT(_execute_REDFT00);  break;  // DCT-I
    case LIQUID_FFT_REDFT10:  q->execute = &FFT(_execute_REDFT10);  break;  // DCT-II
    case LIQUID_FFT_REDFT01:  q->execute = &FFT(_execute_REDFT01);  break;  // DCT-III
    case LIQUID_FFT_REDFT11:  q->execute = &FFT(_execute_REDFT11);  break;  // DCT-IV

    case LIQUID_FFT_RODFT00:  q->execute = &FFT(_execute_RODFT00);  break;  // DST-I
    case LIQUID_FFT_RODFT10:  q->execute = &FFT(_execute_RODFT10);  break;  // DST-II
    case LIQUID_FFT_RODFT01:  q->execute = &FFT(_execute_RODFT01);  break;  // DST-III
    case LIQUID_FFT_RODFT11:  q->execute = &FFT(_execute_RODFT11);  break;  // DST-IV
    default:
        fprintf(stderr,"error: fft_create_plan_r2r_1d(), invalid type, %d\n", q->type);
        exit(1);
    }

    return q;
}

void Ocean::mainComputation() {
	
	// first FFT
	for(int n = 0 ; n < _nx ; n++) { 
		
        // puts heights in _hRf and _hIf
		getSineAmp(n, (double)glutGet(GLUT_ELAPSED_TIME)/1000, &_hRf[n], &_hIf[n]);
		
		_fft = FFT(_ny, _hRf[n], _hIf[n]);
		_fft.calculR();
		_fft.getResult(&_hRf[n], &_hIf[n]);
		
	}
	
	// second one, since it is a 2D FFT
	for(int y = 0 ; y < _ny ; y++) {
		
		int n;
		std::vector<double>::iterator it;
		
        // puts heights in _hRf and _hIf
		for(it = _hRt[y].begin(), n = 0 ; it != _hRt[y].end() ; it++, n++) *it = _hRf[n][y];
		for(it = _hIt[y].begin(), n = 0 ; it != _hIt[y].end() ; it++, n++) *it = _hIf[n][y];
		
		_fft = FFT(_nx, _hRt[y], _hIt[y]);
		_fft.calculR();
		_fft.getResult(&_hRt[y], &_hIt[y]);
		
	}
	
}

void stable_solve ( int n, fftw_real * u, fftw_real * v, fftw_real * u0, fftw_real * v0, fftw_real visc, fftw_real dt ) 
{
	fftw_real x, y, x0, y0, f, r, U[2], V[2], s, t;
	int i, j, i0, j0, i1, j1;

	for (i=0;i<n*n;i++) 
	{ u[i] += dt*u0[i]; u0[i] = u[i]; v[i] += dt*v0[i]; v0[i] = v[i]; }    

	for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n ) 
	   for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n ) 
	   {
	      x0 = n*(x-dt*u0[i+n*j])-0.5f; 
	      y0 = n*(y-dt*v0[i+n*j])-0.5f;
	      i0 = floor(x0); s = x0-i0;
	      i0 = (n+(i0%n))%n;
	      i1 = (i0+1)%n;
	      j0 = floor(y0); t = y0-j0;
	      j0 = (n+(j0%n))%n;
	      j1 = (j0+1)%n;
	      u[i+n*j] = (1-s)*((1-t)*u0[i0+n*j0]+t*u0[i0+n*j1])+                        
			  s *((1-t)*u0[i1+n*j0]+t*u0[i1+n*j1]);
	      v[i+n*j] = (1-s)*((1-t)*v0[i0+n*j0]+t*v0[i0+n*j1])+
			  s *((1-t)*v0[i1+n*j0]+t*v0[i1+n*j1]);
	   }     
	
	for(i=0; i<n; i++)
	  for(j=0; j<n; j++) 
	  {  u0[i+(n+2)*j] = u[i+n*j]; v0[i+(n+2)*j] = v[i+n*j]; }

	FFT(1,u0);
	FFT(1,v0);

	for (i=0;i<=n;i+=2) 
	{
	   x = 0.5f*i;
	   for (j=0;j<n;j++) 
	   {
	      y = j<=n/2 ? (fftw_real)j : (fftw_real)j-n;
	      r = x*x+y*y;
	      if ( r==0.0f ) continue;
	      f = (fftw_real)exp(-r*dt*visc);
	      U[0] = u0[i  +(n+2)*j]; V[0] = v0[i  +(n+2)*j];
	      U[1] = u0[i+1+(n+2)*j]; V[1] = v0[i+1+(n+2)*j];

	      u0[i  +(n+2)*j] = f*( (1-x*x/r)*U[0]     -x*y/r *V[0] );
	      u0[i+1+(n+2)*j] = f*( (1-x*x/r)*U[1]     -x*y/r *V[1] );
	      v0[i+  (n+2)*j] = f*(   -y*x/r *U[0] + (1-y*y/r)*V[0] );
	      v0[i+1+(n+2)*j] = f*(   -y*x/r *U[1] + (1-y*y/r)*V[1] );
	   }    
	}

	FFT(-1,u0); 
	FFT(-1,v0);

	f = 1.0/(n*n);
 	for (i=0;i<n;i++)
	   for (j=0;j<n;j++) 
	   { u[i+n*j] = f*u0[i+(n+2)*j]; v[i+n*j] = f*v0[i+(n+2)*j]; }
} 

static void fastConvolution(float* fastConvResult,float* ir,float* x_bucket,int ir_startindex,int ir_Size,int x_bucket_startindex,int x_Size,int ish0convolution)
{
	int x_fft_input_Size;
	int ir_fft_input_Size;

	//if convolving with h0, do this
	if(ish0convolution == TRUE)
	{
		//set sizes of the arrays that will be passed to FFT function. Both need to be set to 4*N
		x_fft_input_Size = 4*N;
		ir_fft_input_Size = 4*N;
	}
	else
	{
		//set sizes of the arrays that will be passed to FFT function. Both need to be set to 4*N
		x_fft_input_Size = 2*x_Size;
		ir_fft_input_Size = 2*ir_Size;	
	}
		
	//FFT input arrays declared, set to size as stipulated above
	float x_fft_input[x_fft_input_Size];
	float ir_fft_input[ir_fft_input_Size];
		
	//copy over input bucket data to FFT input array for the input signal
	int i;
	for(i=0;i<x_Size;i++)
	{
		x_fft_input[i] = x_bucket[x_bucket_startindex+i];
	}
		
	//copy over input bucket data to FFT input array for the input signal
	for(i=0;i<ir_Size;i++)
	{
		ir_fft_input[i] = ir[ir_startindex+i];
	}
		
	//set sizes of the arrays that will contain the result of the FFTs. These will be double of the input and ir arrays
	int X_Size = 2*x_fft_input_Size;
	int IR_Size = 2*ir_fft_input_Size;
		
	//FFT output arrays declared, set to double the size of the FFT input arrays
	float X[X_Size];
	float IR[IR_Size];
		
	//Call FFT function, write to X_fft_output and IR_fft_output
	FFT(x_fft_input,X,x_fft_input_Size);
	FFT(ir_fft_input,IR,ir_fft_input_Size);
		
	//Declare function that will contain the complex multiplication results
	float complexFastMultResult[X_Size];
		
	//Call complexFastMult function
	complexFastMult(X,IR,complexFastMultResult,X_Size);
		
	//Call IFFT to convert data in complexFastMultResult to the time domain and write it to fastConvResult
	IFFT(complexFastMultResult,fastConvResult,x_fft_input_Size);		
}

int FFT2D(COMPLEX c[][32],int nx,int ny,int dir)
{
   int i,j;
   int m,twopm;
   

   /* Transform the rows */
   if (realx == 0) {
    realx = (double *)malloc(nx * sizeof(double));
    imagx = (double *)malloc(nx * sizeof(double));
    realy = (double *)malloc(ny * sizeof(double));
    imagy = (double *)malloc(ny * sizeof(double));
   }
   //if (real == NULL || imag == NULL)
     // return(FALSE);
   if (!Powerof2(nx,&m,&twopm) || twopm != nx)
      return(FALSE);
   for (j=0;j<ny;j++) {
      for (i=0;i<nx;i++) {
         realx[i] = c[i][j].real;
         imagx[i] = c[i][j].imag;
      }
      FFT(dir,m,realx,imagx);
      for (i=0;i<nx;i++) {
         c[i][j].real = realx[i];
         c[i][j].imag = imagx[i];
      }
   }
   //free(real);
   //free(imag);

   /* Transform the columns */
   //real = (double *)malloc(ny * sizeof(double));
   //imag = (double *)malloc(ny * sizeof(double));
   //if (real == NULL || imag == NULL)
//      return(FALSE);
   if (!Powerof2(ny,&m,&twopm) || twopm != ny)
      return(FALSE);
   for (i=0;i<nx;i++) {
      for (j=0;j<ny;j++) {
         realy[j] = c[i][j].real;
         imagy[j] = c[i][j].imag;
      }
      FFT(dir,m,realy,imagy);
      for (j=0;j<ny;j++) {
         c[i][j].real = realy[j];
         c[i][j].imag = imagy[j];
      }
   }
   //free(real);
   //free(imag);

   return(TRUE);
}

int FFT2D(COMPLEX *c,int nx,int ny,int dir)
{
   int i,j;
   int m,twopm;

   //double* real = new double[nx];
   double* real = malloc(nx*sizeof(double));
   //double* imag = new double[nx];
   double* imag = malloc(nx*sizeof(double));
   
   if (!Powerof2(nx,&m,&twopm) || twopm != nx)
      return(0);
   for (j=0;j<ny;j++) {
      for (i=0;i<nx;i++) {
         real[i] = c[j*nx+i].real;
         imag[i] = c[j*nx+i].imag;
      }
      FFT(real,imag,m,dir);
      for (i=0;i<nx;i++) {
         c[j*nx+i].real = real[i];
         c[j*nx+i].imag = imag[i];
      }
   }
   //delete[] real;
   free(real);
//delete[] imag;
	free(imag);
   
   //real = new double[ny];
  real = malloc(ny*sizeof(double));
   //imag = new double[ny];
   imag = malloc(ny*sizeof(double));
   
   if (!Powerof2(ny,&m,&twopm) || twopm != ny)
      return(0);
   for (i=0;i<nx;i++) {
      for (j=0;j<ny;j++) {
         real[j] = c[j*nx+i].real;
         imag[j] = c[j*nx+i].imag;
      }
      FFT(real,imag,m,dir);
      for (j=0;j<ny;j++) {
         c[j*nx+i].real = real[j];
         c[j*nx+i].imag = imag[j];
      }
   }
   //delete[] real;
   free(real);
//delete[] imag;
	free(imag);

   return(1);
}

int main(void)
{
  printf("If rows come in identical pairs, then everything works.\n");
  
  cpx a[8] = {0, 1, cpx(1,3), cpx(0,5), 1, 0, 2, 0};
  cpx b[8] = {1, cpx(0,-2), cpx(0,1), 3, -1, -3, 1, -2};
  cpx A[8];
  cpx B[8];
  FFT(a, A, 1, 8, 1);
  FFT(b, B, 1, 8, 1);
  
  for(int i = 0 ; i < 8 ; i++)
  {
    printf("%7.2lf%7.2lf", A[i].a, A[i].b);
  }
  printf("\n");
  for(int i = 0 ; i < 8 ; i++)
  {
    cpx Ai(0,0);
    for(int j = 0 ; j < 8 ; j++)
    {
      Ai = Ai + a[j] * EXP(j * i * two_pi / 8);
    }
    printf("%7.2lf%7.2lf", Ai.a, Ai.b);
  }
  printf("\n");
  
  cpx AB[8];
  for(int i = 0 ; i < 8 ; i++)
    AB[i] = A[i] * B[i];
  cpx aconvb[8];
  FFT(AB, aconvb, 1, 8, -1);
  for(int i = 0 ; i < 8 ; i++)
    aconvb[i] = aconvb[i] / 8;
  for(int i = 0 ; i < 8 ; i++)
  {
    printf("%7.2lf%7.2lf", aconvb[i].a, aconvb[i].b);
  }
  printf("\n");
  for(int i = 0 ; i < 8 ; i++)
  {
    cpx aconvbi(0,0);
    for(int j = 0 ; j < 8 ; j++)
    {
      aconvbi = aconvbi + a[j] * b[(8 + i - j) % 8];
    }
    printf("%7.2lf%7.2lf", aconvbi.a, aconvbi.b);
  }
  printf("\n");
  
  return 0;
}

/*-------------------------------------------------------------------------
   Perform a 2D FFT inplace given a complex 2D array
   The direction dir, 1 for forward, -1 for reverse
   The size of the array (nx,ny)
   Return false if there are memory problems or
      the dimensions are not powers of 2
*/
procStatus CFFTMachine::FFT2D( COMPLEX **c,int nx,int ny,int dir )
{
   int i,j;
   int m,twopm;
   double *real,*imag;

   /* Transform the rows */
   real = (double *)malloc(nx * sizeof(double));
   imag = (double *)malloc(nx * sizeof(double));
   if (real == NULL || imag == NULL)
      return(eMemAllocErr);
   if (!Powerof2(nx,&m,&twopm) || twopm != nx)
      return(eInvalideArg);
   for (j=0;j<ny;j++) {
      for (i=0;i<nx;i++) {
         real[i] = c[i][j].real;
         imag[i] = c[i][j].imag;
      }
      FFT(dir,m,real,imag);
      for (i=0;i<nx;i++) {
         c[i][j].real = real[i];
         c[i][j].imag = imag[i];
      }
   }
   free(real);
   free(imag);

   /* Transform the columns */
   real = (double *)malloc(ny * sizeof(double));
   imag = (double *)malloc(ny * sizeof(double));
   if (real == NULL || imag == NULL)
      return(eMemAllocErr);
   if (!Powerof2(ny,&m,&twopm) || twopm != ny)
      return(eInvalideArg);
   for (i=0;i<nx;i++) {
      for (j=0;j<ny;j++) {
         real[j] = c[i][j].real;
         imag[j] = c[i][j].imag;
      }
      FFT(dir,m,real,imag);
      for (j=0;j<ny;j++) {
         c[i][j].real = real[j];
         c[i][j].imag = imag[j];
      }
   }
   free(real);
   free(imag);

   return(eNormal);
}

//int FFT2D(COMPLEX **c,int nx,int ny,int dir)
int FFT2D(Complex	 c[WAVE_SIZE][WAVE_SIZE], int nx, int ny, int dir)
{
	int i, j;
	int m, twopm;
	double *real, *imag;

	/* Transform the rows */
	real = (double *)malloc(nx * sizeof(double));
	imag = (double *)malloc(nx * sizeof(double));
	if (real == NULL || imag == NULL)
		return  0;
	if (!Powerof2(nx, &m, &twopm) || twopm != nx)
		return  0;
	for (j = 0; j<ny; j++) {
		for (i = 0; i<nx; i++) {
			real[i] = c[i][j].real;
			imag[i] = c[i][j].imag;
		}
		FFT(dir, m, real, imag);
		for (i = 0; i<nx; i++) {
			c[i][j].real = real[i];
			c[i][j].imag = imag[i];
		}
	}
	free(real);
	free(imag);

	/* Transform the columns */
	real = (double *)malloc(ny * sizeof(double));
	imag = (double *)malloc(ny * sizeof(double));
	if (real == NULL || imag == NULL)
		return 0;
	if (!Powerof2(ny, &m, &twopm) || twopm != ny)
		return  0;
	for (i = 0; i<nx; i++) {
		for (j = 0; j<ny; j++) {
			real[j] = c[i][j].real;
			imag[j] = c[i][j].imag;
		}
		FFT(dir, m, real, imag);
		for (j = 0; j<ny; j++) {
			c[i][j].real = real[j];
			c[i][j].imag = imag[j];
		}
	}
	free(real);
	free(imag);

	return		1;
}

void FFT(cpx *in, cpx *out, int step, int size, int dir) {
    if (size < 1) return;
    if (size == 1) {
        out[0] = in[0];
        return;
    }
    FFT(in, out, step * 2, size / 2, dir);
    FFT(in + step, out + size / 2, step * 2, size / 2, dir);
    for (int i = 0; i < size / 2; i++) {
        cpx even = out[i];
        cpx odd = out[i + size / 2];
        out[i] = even + EXP(dir * two_pi * i / size) * odd;
        out[i + size / 2] = even + EXP(dir * two_pi * (i + size / 2) / size) * odd;
    }
}

void FFT(double *x_r,double *x_i,double *y_r,double *y_i,int N)
{
     if(N == 1)
     {
          y_r[0] = x_r[0];
          y_i[0] = x_i[0];
          return;
     }
     int k;
     double *u_r,*u_i,*v_r,*v_i,w_r,w_i;
     
     u_r = (double *) malloc (N*sizeof(double));
     u_i = (double *) malloc (N*sizeof(double));
     
     v_r = (double *) malloc (N*sizeof(double));
     v_i = (double *) malloc (N*sizeof(double));
     
     for(k = 0 ; k < N/2 ; k++)
     {
             u_r[k] = x_r[2*k];
             u_i[k] = x_i[2*k];
             u_r[k+N/2] = x_r[2*k+1];
             u_i[k+N/2] = x_i[2*k+1];
     }
     
     
     FFT(u_r,u_i,v_r,v_i,N/2);
     FFT(u_r+N/2,u_i+N/2,v_r+N/2,v_i+N/2,N/2);
     
     for(k = 0 ; k < N/2 ; k++)
     {
           w_r = cos(-k*2*M_PI/N);
           w_i = sin(-k*2*M_PI/N);
           
           //y[k] = v[k] + w*v[k+N/2] & y[k+N/2] = v[k] - w*v[k+N/2]
           y_r[k] = v_r[k] + w_r*v_r[k+N/2] - w_i*v_i[k+N/2];
           y_i[k] = v_i[k] + w_r*v_i[k+N/2] + w_i*v_r[k+N/2];
           y_r[k+N/2] = v_r[k] - (w_r*v_r[k+N/2] - w_i*v_i[k+N/2]);
           y_i[k+N/2] = v_i[k] - (w_r*v_i[k+N/2] + w_i*v_r[k+N/2]);
     }
     
     free(u_r);
     free(u_i);
     free(v_r);
     free(v_i);
     
     return;
}

/* EXPORT-> Realft: apply fft to real s */
void Realft (Vector s)
{
   int n, n2, i, i1, i2, i3, i4;
   double xr1, xi1, xr2, xi2, wrs, wis;
   double yr, yi, yr2, yi2, yr0, theta, x;

   n=VectorSize(s) / 2; n2 = n/2;
   theta = PI / n;
   FFT(s, FALSE);
   x = sin(0.5 * theta);
   yr2 = -2.0 * x * x;
   yi2 = sin(theta); yr = 1.0 + yr2; yi = yi2;
   for (i=2; i<=n2; i++) {
      i1 = i + i - 1;      i2 = i1 + 1;
      i3 = n + n + 3 - i2; i4 = i3 + 1;
      wrs = yr; wis = yi;
      xr1 = (s[i1] + s[i3])/2.0; xi1 = (s[i2] - s[i4])/2.0;
      xr2 = (s[i2] + s[i4])/2.0; xi2 = (s[i3] - s[i1])/2.0;
      s[i1] = xr1 + wrs * xr2 - wis * xi2;
      s[i2] = xi1 + wrs * xi2 + wis * xr2;
      s[i3] = xr1 - wrs * xr2 + wis * xi2;
      s[i4] = -xi1 + wrs * xi2 + wis * xr2;
      yr0 = yr;
      yr = yr * yr2 - yi  * yi2 + yr;
      yi = yi * yr2 + yr0 * yi2 + yi;
   }
   xr1 = s[1];
   s[1] = xr1 + s[2];
   s[2] = 0.0;
}

/****************************************************
	IFFT()

	参数：

		FD为频域值
		TD为时域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数利用快速傅立叶变换实现傅立叶反变换
****************************************************/
void IFFT(COMPLEX * FD, COMPLEX * TD, int power)
{
	int i, count;
	COMPLEX *x;

	/*计算傅立叶反变换点数*/
	count=1<<power;

	/*分配运算所需存储器*/
	x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);

	/*将频域点写入存储器*/
	memcpy(x,FD,sizeof(COMPLEX)*count);
	
	/*求频域点的共轭*/
	for(i=0;i<count;i++)
		x[i].im = -x[i].im;

	/*调用FFT*/
	FFT(x, TD, power);

	/*求时域点的共轭*/
	for(i=0;i<count;i++)
	{
		TD[i].re /= count;
		TD[i].im = -TD[i].im / count;
	}

	/*释放存储器*/
	free(x);
}

void doit(int iter, struct problem *p)
{
     int i;
     if (p->kind == PROBLEM_COMPLEX) {
	  bench_complex *in = (bench_complex *) p->in;
	  int two_n = 2 * p->n[0];

	  if (p->sign > 0)
	       for (i = 0; i < iter; ++i)
		    FFT(in, &two_n);
	  else
	       for (i = 0; i < iter; ++i)
		    IFFT(in, &two_n);
     }
     else /* PROBLEM_REAL */ {
	  bench_real *in = (bench_real *) p->in;
	  int n = p->n[0];

	  if (p->sign < 0)
	       for (i = 0; i < iter; ++i)
		    REAL_FFT(in, &n);
	  else
	       for (i = 0; i < iter; ++i)
		    REAL_IFFT(in, &n);
     }
}

void RealFFTI(const ColumnVector& A, const ColumnVector& B, ColumnVector& U)
{
   // inverse of a Fourier transform of a real series
   Tracer trace("RealFFTI");
   REPORT
   const int n21 = A.Nrows();                     // length of arrays
   if (n21 != B.Nrows() || n21 == 0)
      Throw(ProgramException("Vector lengths unequal or zero", A, B));
   const int n2 = n21 - 1;  const int n = 2 * n2;  int i = n2 - 1;

   ColumnVector X(n2), Y(n2);
   Real* a = A.Store(); Real* b = B.Store();  // first els of A and B
   Real* an = a + n2;   Real* bn = b + n2;    // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + i;    Real* yn = y + i;     // last els of X and Y

   Real hn = 0.5 / n2;
   *x++  = hn * (*a + *an);  *y++  = - hn * (*a - *an);
   a++; an--; b++; bn--;
   int j = -1;  i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am - c * bp; Real sbpcam = s * bp + c * am;
      *x++  =  hn * ( ap + samcbp); *y++  =  - hn * ( bm + sbpcam);
      *xn-- =  hn * ( ap - samcbp); *yn-- =  - hn * (-bm + sbpcam);
   }
   FFT(X,Y,X,Y);             // have done inverting elsewhere
   U.ReSize(n); i = n2;
   x = X.Store(); y = Y.Store(); Real* u = U.Store();
   while (i--) { *u++ = *x++; *u++ = - *y++; }
}