Fast Fourier Transformation

In mathematics, the Fourier transform (often abbreviated FT) is an operation that transforms one complex-valued function of a real variable into another. In such applications as signal processing, the domain of the original function is typically time and is accordingly called the time domain. The domain of the new function is typically called the frequency domain, and the new function itself is called the frequency domain representation  of the original function. It describes which frequencies are present in the original function. This is analogous to describing a musical chord in terms of the notes being played. In effect, the Fourier transform decomposes a function into oscillatory  functions. The term Fourier transform refers both to the frequency domain representation of a function, and to the process or formula that "transforms" one function into the other.

The Fourier transform and its generalizations are the subject of Fourier analysis. In this specific case, both the time and frequency domains are unbounded linear continua. It is possible to define the Fourier transform of a function of several variables, which is important for instance in the physical study of wave motion and optics. It is also possible to generalize the Fourier transform on discrete structures such as finite groups, efficient computation of which through a fast Fourier transform is essential for high-speed computing.
[Source:Wikipedia]


A short OpenCL based demonstration version of the FFT


// ----------------- OpenCL Kernel -----------------
const char *KernelSource = "                                        
#pragma OPENCL EXTENSION cl_khr_fp64 : enable                       
__kernel void square(                                               
    __global double* input,                                          
    __global double* RealOutput,
    __global double* ImgOutput,
    const double k,
    const int N)  
{
    const double PI = 3.14159265;
    int n = get_global_id(0); 

    if(n < N) {
        double value = -2.0*PI*n*k/(N/2.0);

        RealOutput[n] = input[n] * cos(value);                                               
        ImgOutput[n]  = input[n] * sin(value);
    }
}                                                                   
\n";

// ----------------- Host Calculation -----------------
const double PI = 3.14159265;

double RealSumOdd  = 0.0;
double RealSumEven = 0.0;

double ImgSumOdd   = 0.0;
double ImgSumEven  = 0.0;

#pragma omp parallel sections num_threads(2)
{
    #pragma omp for reduction(+: RealSumOdd,ImgSumOdd)
    for (int n=0;n<N;n+=2) { // odd
        RealSumOdd  += RealOutput[i];
        ImgSumOdd   += ImgOutput[i];
    }

    #pragma omp section

    #pragma omp for reduction(+: RealSumEven,ImgSumEven)
    for (int n=1;n<N;n+=2) { // even
        RealSumEven += RealOutput[i];
        ImgSumEven  += ImgOutput[i];
    }
}

double RealSum = RealSumEven + cos(-2.0*PI*k/N) * RealSumOdd;
double ImgSum  = ImgSumEven  + sin(-2.0*PI*k/N) * ImgSumOdd;