1 Reply Latest reply on Dec 2, 2013 2:04 PM by bragadeesh

    Bad results with a simple FFT 1D OpenCL example

    chris10

      I try to do a simple FFT 1D with clMath OpenCL libray. For this, I am using a cosinus signal with a frequency equal to 10 Hz and a sampling frequency of sizex*frequency_signal with sizex the number of sampling points.

       

      Into my results, I can't get to have the dirac impulsion at f=10 Hz. For the first frequencies (k*f_sampling/sizex with k =0 to sizex), I have the following outputs (column 1: k*f_sampling/size x and column 2: X_k)

       

       

          0.000000 -5.169492e-06

          10.000000 -5.169449e-06

          20.000000 2.498745e-02

          30.000000 2.499508e-02

          40.000000 2.832322e-06

          50.000000 9.587315e-07

          60.000000 4.648561e-07

          70.000000 9.241408e-08

          80.000000 1.400218e-07

          90.000000 3.699876e-07

          100.000000 2.663564e-07

          110.000000 2.523218e-07

          120.000000 1.631196e-07

          130.000000 7.783362e-08

          140.000000 7.932793e-08

          150.000000 9.296434e-08

          160.000000 1.039581e-07

          170.000000 7.396823e-08

          180.000000 6.584698e-08

          190.000000 1.468647e-08

          200.000000 1.694581e-08

          210.000000 3.367836e-08

          220.000000 3.538253e-08

          ...

       

      As you can see,there seems to be a dirac (I say that because of the "e-02" value compared to the others values but I'm not sure) between f=20 Hz and f=30 Hz and not for 10 Hz as expected.

      Here's the code :

       

      #include "clFFT.h"
      #include <stdio.h>
      #include <stdlib.h>
      #include <string.h>
      #include <math.h>
      
      /////////////////////////////////////
      // OpenCL FFT 1D function ///////////
      /////////////////////////////////////
      
      int FftOpenCL(float *tab, const char* direction, int Ng)
      {
       cl_int err;
          cl_platform_id platform = 0;
          cl_device_id device = 0;
          cl_context_properties props[3] = { CL_CONTEXT_PLATFORM, 0, 0 };
          cl_context ctx = 0;
          cl_command_queue queue = 0;
          cl_mem bufX;
          int ret = 0;
          size_t N = Ng/2;
      
          /* FFT library realted declarations */
          clfftPlanHandle planHandle;
          clfftDim dim = CLFFT_1D;
          size_t clLengths[1] = {N};
      
          /* Setup OpenCL environment. */
          err = clGetPlatformIDs( 1, &platform, NULL );
          err = clGetDeviceIDs( platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL );
      
          props[1] = (cl_context_properties)platform;
          ctx = clCreateContext( props, 1, &device, NULL, NULL, &err );
          queue = clCreateCommandQueue( ctx, device, 0, &err );
      
          /* Setup clFFT. */
          clfftSetupData fftSetup;
          err = clfftInitSetupData(&fftSetup);
          err = clfftSetup(&fftSetup);
      
          /* Prepare OpenCL memory objects and place data inside them. */
          bufX = clCreateBuffer( ctx, CL_MEM_READ_WRITE, N * 2 * sizeof(float), NULL, &err );
      
          err = clEnqueueWriteBuffer( queue, bufX, CL_TRUE, 0,
          N * 2 * sizeof(float), tab, 0, NULL, NULL );
      
          /* Create a default plan for a complex FFT. */
          err = clfftCreateDefaultPlan(&planHandle, ctx, dim, clLengths);
      
          /* Set plan parameters. */
          err = clfftSetPlanPrecision(planHandle, CLFFT_SINGLE);
          err = clfftSetLayout(planHandle, CLFFT_COMPLEX_INTERLEAVED, CLFFT_COMPLEX_INTERLEAVED);
          err = clfftSetResultLocation(planHandle, CLFFT_INPLACE);
      
          /* Bake the plan. */
          err = clfftBakePlan(planHandle, 1, &queue, NULL, NULL);
      
      
          if(strcmp(direction,"forward") == 0) 
          {
          /* Execute the plan. */
          err = clfftEnqueueTransform(planHandle, CLFFT_FORWARD, 1, &queue, 0, NULL, NULL, &bufX, NULL, NULL);
          }
          
          /* Wait for calculations to be finished. */
          err = clFinish(queue);
      
          /* Fetch results of calculations. */
          err = clEnqueueReadBuffer( queue, bufX, CL_TRUE, 0, N * 2 * sizeof(float), tab, 0, NULL, NULL );
      
          /* Release OpenCL memory objects. */
          clReleaseMemObject( bufX );
      
          /* Release the plan. */
          err = clfftDestroyPlan( &planHandle );
      
          /* Release clFFT library. */
          clfftTeardown( );
      
          /* Release OpenCL working objects. */
          clReleaseCommandQueue( queue );
          clReleaseContext( ctx );
      
          return ret;
      
      }
      
      
      int main(void) {
      
          int i;
          
          // temporal array
          float *x_temp;
          
          // signal array and FFT output array
          float *Array;
          
          // number of sampling points
          int sizex = 10000;
      
          float h = 0;
          
          // signal frequency
          float frequency_signal = 10;
          
          // size x points between 0 and T_signal
          float frequency_sampling = sizex*frequency_signal;
          
          // step = T_sampling
          float step = 1.0/frequency_sampling;
      
          FILE *Fft1D;
      
          // Allocation of Array
      
          Array = (float*) malloc(sizex*sizeof(float));
          x_temp = (float*) malloc(sizex*sizeof(float));
      
          // Initialization of 1D ArrayInput
      
          for(i=0; i<sizex; i++)
          {
              Array[i] = cos(2*M_PI*frequency_signal*h);
              x_temp[i] = h; 
              h = h + step;
          }
      
          if (FftOpenCL(Array,"forward", sizex) == 0)
              printf("Fft passed !\n");;
      
          // Printf ArrayOutput - Multiply FFT Array result by period T_sampling
      
          Fft1D = fopen("Fft1D.dat","w");
      
          for(i=0; i<sizex; i++)
              fprintf(Fft1D,"%f %e\n", i*frequency_sampling/sizex, 1.0/frequency_sampling*Array[i]);
      
          fclose(Fft1D);
      
          return 0;
      
      }
      
      

       

      Anyone could see what's wrong ?

      Thanks

        • Re: Bad results with a simple FFT 1D OpenCL example
          bragadeesh

          Well, one immediate thing I see is that you need to have the "imaginary" component in your "Array" buffer since you are asking the library to compute a complex FFT. So your memory allocation should be twice as big and you should write values in the RIRIRI.... format where R is real and I is imaginary.

          Also, you should sample at the rate of 1/sizex and not 1/(sizex*frequency_signal) because the sampling frequency for your 'sizex' transform is 1/sizex. Sampling frequency is different from the test signal frequency you are trying to create. Assuming your time period is 1 second and number of samples is 10000, then your sampling frequency is 10 kHz. Your signal frequency can be no greater than 5 kHz for the Fourier transform to detect it according to the Nyquist theorem. Th test case of 10 Hz is valid.