自学教程：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 methodFFT(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++; }}

void RealFFT(const ColumnVector& U, ColumnVector& X, ColumnVector& Y){   // Fourier transform of a real series   Tracer trace("RealFFT");   REPORT   const int n = U.Nrows();                     // length of arrays   const int n2 = n / 2;   if (n != 2 * n2)      Throw(ProgramException("Vector length not multiple of 2", U));   ColumnVector A(n2), B(n2);   Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;   while (i--) { *a++ = *u++; *b++ = *u++; }   FFT(A,B,A,B);   int n21 = n2 + 1;   X.ReSize(n21); Y.ReSize(n21);   i = n2 - 1;   a = A.Store(); b = B.Store();              // first els of A and B   Real* an = a + i; Real* bn = b + i;        // last els of A and B   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y   Real* xn = x + n2; Real* yn = y + n2;      // last els of X and Y   *x++ = *a + *b; *y++ = 0.0;                // first complex element   *xn-- = *a++ - *b++; *yn-- = 0.0;          // last complex element   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++  =  0.5 * ( ap + samcbp); *y++  =  0.5 * ( bm + sbpcam);      *xn-- =  0.5 * ( ap - samcbp); *yn-- =  0.5 * (-bm + sbpcam);   }}

int main(void){  MONO_PCM pcm0;  int n, k, N;  double *x_real, *x_imag;    mono_wave_read(&pcm0, "ex2_1.wav"); /* WAVEファイルからモノラルの音デ
C++ FF_ARRAY_ELEMS函数代码示例
C++ FFMIN函数代码示例