about frequency spectrum computed from cuFFT

Hi,

For confirming frequency spectrum of cuFFT with that of FFT of Matlab, I did cuFFT for a sin function. The result of the spectrum of cuFFT is different from that of FFT of Matlab. My codes are below:

#include “cuda_runtime.h”
#include “device_launch_parameters.h”
#include “cufft.h”

#include
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define GNUPLOT_NAME “C:/gnuplot/bin/gnuplot -persist”

#define NX 256 // number of fft
typedef float2 Complex;

int main()
{
#ifdef WIN32
FILE *pipe = _popen(GNUPLOT_NAME, “w”);
#else
FILE *pipe = popen(GNUPLOT_NAME, “w”);
#endif
//Complex *x;
float arg3, pi = 3.1415926, dt=0.1;
int size = 128;
float t;
int fc = 35;
cufftComplex *d_data, *x;

// define *h_data with cufftComplex or float2 has same result
float2 *h_data=0;


// alloc memory
x = (cufftComplex *)malloc(sizeof(cufftComplex) * size);


h_data = (float2 *)malloc(sizeof(float2) * size);

for (int i = 0; i < size; i++) {		
	printf("%d %lf \n", i,t);		
	x[i].x = sin((float)(i+1) / (float)5);
	x[i].y = 0.0f;
}

// alloc memory for device variables and copy variables from host to device
cudaMalloc((void **)&d_data, sizeof(cufftComplex)*size);
cudaMemcpy(d_data, x, sizeof(cufftComplex)*size, cudaMemcpyHostToDevice);

// set up fft plan
cufftHandle plan;	
cufftPlan1d(&plan, sizeof(cufftComplex)*size, CUFFT_C2C, 1);

// do fft
cufftExecC2C(plan, d_data, d_data, CUFFT_FORWARD);


// copy varibales from device to host for plotting
cudaMemcpy(h_data, d_data, sizeof(cufftComplex)*size, cudaMemcpyDeviceToHost);

cudaDeviceSynchronize();

    // plot sin wave and frequency spectrum by using gnuplot
if (pipe != NULL)
{
	
	fprintf(pipe, "set multiplot layout 1,2 \n");         // set multiple plots		
	fprintf(pipe, "plot '-' with linespoints \n");        // plot type
											
	for (int i = 0; i < size; i++)             // loop over the data [0,...,9]
		fprintf(pipe, "%lf \n", x[i].x);           // data terminated with 
	   
	fprintf(pipe, "%s\n", "e");             // termination character
	fflush(pipe);                           // flush the pipe
	

	fprintf(pipe, "plot '-' with linespoints\n"); // plot type
											
	for (int i = 0; i < size; i++)             // loop over the data [0,...,9]
		fprintf(pipe, "%lf \n", cuCabsf(h_data[i]));      // cuCabsf		
	
	fprintf(pipe, "%s\n", "e");             // termination character
	fflush(pipe);                           // flush the pipe
	
	fprintf(pipe, "unset multiplot \n");      
	
	// wait for key press
	std::cin.clear();
	std::cin.ignore(std::cin.rdbuf()->in_avail());
	std::cin.get();

#ifdef WIN32
_pclose(pipe);
#else
pclose(pipe);
#endif
}
else
std::cout << “Could not open pipe” << std::endl;

cufftDestroy(plan);
free(x);
free(h_data);
cudaFree(d_data);

return 0;

}


So, please help me find out what causes the frequency spectrum of cuFFT is different from that of FFT by Matlab (spectrum in Matlab is correct).


 Thanks very much in advanced.

Dawn

You say “spectrum in Matlab is correct” but you don’t show any matlab code.

The FFT (or spectrum) of a sine wave should just be single spike at the frequency of the sinewave.

In what way is the cuFFT data not correct?

Do you get a single frequency spike?
Is it at the correct frequency?

You may simply need to normalize the results in cuFFT - divide by the size of the transform.

Thank you for your help.

Yes, by Matlab the single frequency spike is located at 5Hz while that is located in 30Hz with cuFFT. I know spike should be in 5Hz. So, can you see what is wrong with my cuda codes?

Dawn

What makes you think this gives you a frequency of 5Hz:

x[i].x = sin((float)(i+1) / (float)5);

?

What is your intention for your time step/sample interval? I don’t see it defined anywhere in your program, and the variable:

float t;

isn’t set to anything (although you are printing it out).

You also have a float dt variable, set to 0.1, but not sure what you mean by that, it is not used in your program. a 0.1 second time step implied?

??

After carefully cleaning my cuda codes, I repost codes as below:


#include “cuda_runtime.h”
#include “device_launch_parameters.h”
#include “cufft.h”

#include
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define GNUPLOT_NAME “C:/gnuplot/bin/gnuplot -persist”

typedef float2 Complex;

int main()
{
#ifdef WIN32
FILE *pipe = _popen(GNUPLOT_NAME, “w”);
#else
FILE *pipe = popen(GNUPLOT_NAME, “w”);
#endif

float arg3, pi = 3.1415926, dt=0.1;
int size = 128;

cufftComplex *d_data, *x;

float2 *h_data=0;

// alloc memory
x = (cufftComplex *)malloc(sizeof(cufftComplex) * size);

h_data = (float2 *)malloc(sizeof(float2) * size);

for (int i = 0; i < size; i++) {
printf("%d %lf \n", i,t);
x[i].x = sin((float)(i+1) / (float)5);
x[i].y = 0.0f;
}

// alloc memory for device variables and copy variables from host to device
cudaMalloc((void **)&d_data, sizeof(cufftComplex)*size);
cudaMemcpy(d_data, x, sizeof(cufftComplex)*size, cudaMemcpyHostToDevice);

// set up fft plan
cufftHandle plan;
cufftPlan1d(&plan, sizeof(cufftComplex)*size, CUFFT_C2C, 1);

// do fft
cufftExecC2C(plan, d_data, d_data, CUFFT_FORWARD);

// copy varibales from device to host for plotting
cudaMemcpy(h_data, d_data, sizeof(cufftComplex)*size, cudaMemcpyDeviceToHost);

cudaDeviceSynchronize();

// plot sin wave and frequency spectrum by using gnuplot
if (pipe != NULL)
{

fprintf(pipe, “set multiplot layout 1,2 \n”); // set multiple plots
fprintf(pipe, “plot ‘-’ with linespoints \n”); // plot type

for (int i = 0; i < size; i++) // loop over the data [0,…,size-1]
fprintf(pipe, “%lf \n”, x[i].x);

fprintf(pipe, “%s\n”, “e”); // termination character
fflush(pipe); // flush the pipe

fprintf(pipe, “plot ‘-’ with linespoints\n”); // plot type

for (int i = 0; i < size; i++) // loop over the data [0,…,size-1]
fprintf(pipe, “%lf \n”, cuCabsf(h_data[i])); // cuCabsf

fprintf(pipe, “%s\n”, “e”); // termination character
fflush(pipe); // flush the pipe

fprintf(pipe, “unset multiplot \n”);

// wait for key press
std::cin.clear();
std::cin.ignore(std::cin.rdbuf()->in_avail());
std::cin.get();

#ifdef WIN32
_pclose(pipe);
#else
pclose(pipe);
#endif
}
else
std::cout << “Could not open pipe” << std::endl;

cufftDestroy(plan);
free(x);
free(h_data);
cudaFree(d_data);

return 0;
}


My Matlab codes are:


t=1:128;

x=sin(t./5);

subplot(121)
plot(t,x);

xfft=fft(x);

subplot(122)
plot(abs(xfft));


The frequency spectrum of sin wave by Matlab has two peak frequencies at k=5 and k=125, respectively, while that of cuFFT has one at k=30. Each k number is 1/128 wave numbers. So, I think if you know Fourier coefficients of cuFFT are correct or not. Or, there is some proportional scale I need to add in.

Thanks very much in advanced.

Dawn

now your posted CUDA/CUFFT code does not compile.

In any event, this is incorrect in your code:

cufftPlan1d(&plan, sizeof(cufftComplex)*size, CUFFT_C2C, 1);

the second parameter should be just size, not sizeof(cufftComplex)*size. It is supposed to be the size of the transform in elements, not bytes.

When I make that fix, I get spikes at offset 4 and 124 in the cufft output, which I think matches what you are seeing with Matlab. The spike at 124 is just the symmetric version of the spike at 4. A C2C transform of purely real data produces symmetric data.

In the future, I would recommend:

  1. Use proper CUDA error checking
  2. Use proper CUFFT error checking - check the return value of each API call.
  3. As a final test, run your code with cuda-memcheck and make sure there are no reported errors.

Thank you very much for your instruction. It is y mistake - careless, size is an integer rather than a cuFFTcomplex. So, my problem has been solved.

Thanks a lot one more time!

Dawn

Hi, I used cuda error check funtion to check cufftPlan1d as below:

  checkCudaErrors(cufftPlan1d(&plan, sizeof(cufftComplex)*size, CUFFT_C2C, 1));


 It passed the compilation without error message and is executed successfully. But, the frequency spectrum is incorrect.  Is there any way that can check this kind of syntax error - wrong data type and wrong data dimension? 

  In c/c++, they can be automatically checked for syntax error.


  Thanks in advance!

Dawn

First of all, this can easily be made into ordinary c/c++ code. CUFFT activity does not need to be in a .cu file. You can put it into a c/c++ file and compile it with an ordinary c/c++ compiler. If you need examples, look at the CUFFT sample codes.

Having said that, then, I don’t really believe your claim that “In c/c++, they can be automatically checked for syntax error.” If you believe that, please provide an example.

My suggestion is that this type of error cannot be easily caught at compile-time, but can be caught at run-time.

The call itself will not return a runtime error, because of the asynchronous nature of the underlying activity. But subsequent calls will return an error of some sort, if you do rigorous error checking on all the calls.

And, you can also run the (broken) code with cuda-memcheck, and it will spew out errors.

Hi,

 Thank you for your instruction and kindly guidance. CUFFT works well for me. Right now, I meet a FILE I/O for 2D array problem. I tried to read in a 2D data file and save it as a new 2D data file. Then, I read in the new 2D data file to make a plot as an image by using gnuplot. The plotting results are the 2D array from original data file can be plotted as an image while the 2D array from new data file cannot be. Could you please indicate the problem in my codes? My codes are below:

#define GNUPLOT_NAME “C:/gnuplot/bin/gnuplot -persist”

int main()
{
#ifdef WIN32
FILE *pipe = _popen(GNUPLOT_NAME, “w”);
#else
FILE *pipe = popen(GNUPLOT_NAME, “w”);
#endif

    int xsize = 20, ysize = 41; 
    float **din; 

    // alloc memory 
    din = (float **)malloc(sizeof(float *) * ysize); 


    for (int i = 0; i < ysize; i++) { 
            din[i] = (float *)malloc(sizeof(float) * xsize);                 
    } 


    // read in 2D data 
    FILE *fp, *fp1;         
    fp = fopen("wedgegas.txt", "r"); 
            
    // read in data 
    if (fp == NULL) { 
            printf("fp is NULL!\n"); 
            exit(1); 
    } 

    for (int i = 0; i < ysize; i++) { 
            printf("i=%d \n ", i); 
            for (int j = 0; j < xsize; j++) { 
                    fscanf(fp, "%f ", &din[i][j]); 
                    printf("%f ", din[i][j]); 
            } 
            fscanf(fp, "\n"); 
            printf("\n "); 
    } 

     
    /* write out as a new data file */ 
    fp1 = fopen("inversedwedgegas.txt", "w"); 
    if (fp1 == NULL) { 
            printf("fp1 is NULL!\n"); 
            exit(1); 
    } 
    for (int i = 0; i < ysize; i++) {                
            for (int j = 0; j < xsize; j++) {  

                    fprintf(fp1, "%f ", din[i][j]); 
                    
            } 
            fprintf(fp1, "\n");               
    } 


  /* ---------------------------- plot results ------------------------------------------------------- */ 

    if (pipe != NULL) 
    { 
            // plot original data file        
            fprintf(pipe, "set multiplot layout 1,2 \n");         // set multiple plots                                                                                                                                         


            fprintf(pipe, "plot 'wedgegas.txt' matrix with image \n"); // plot type           

                                                                                                                             

            fprintf(pipe, "%s \n", "e");             // termination character, end a plot 
            fflush(pipe);                           // flush the pipe 

             

            // plot the new data file                                                                                 
            fprintf(pipe, "plot 'inversedwedgegas.txt' matrix with image \n"); // plot type 

            fprintf(pipe, "%s\n", "e");             // termination character 
            fflush(pipe);                           // flush the pipe 
                                                                                                                       
            // wait for key press 
            std::cin.clear(); 
            std::cin.ignore(std::cin.rdbuf()->in_avail()); 
            std::cin.get(); 

#ifdef WIN32
_pclose(pipe);
#else
pclose(pipe);
#endif
}
else
std::cout << “Could not open pipe” << std::endl;

    fclose(fp); 
    fclose(fp1); 



    return 0; 

}


 The error message show me "matrix does not represent a grid".

 Thank you very much in advance!

Dawn