Implementation of 2-D FFT algorithm--base 2 fast two-dimensional Fourier transform based on GPU

implementation of 2-D FFT algorithm--base 2 fast two-dimensional Fourier transform based on GPU

The first one-dimensional FFT of the GPU implementation (FFT algorithm implementation-based on the GPU base 2 fast Fourier transform), and then I need to do a second-dimensional FFT, probably the following ideas.

The first thing to look at is definitely the formula:

As described in the above formula, the 2-D FFT only needs to be split into a line of FFT, and the column FFT on the line, where my implementation is assumed to be the origin in F (0,0), because my code needs the origin in the center, so in the end I moved the origin to the center.

Here is the pseudo-code for the 2-D FFT of the Origin F (0,0):

    C2dfft
    //Is executed 2DFFT is a n*n matrix, stored in row order in source_2d/
    /horizontal FFT for
    (int i=0;i<n;i++)
    {
        fft1 ( &source_2d[i*n],&source_2d_1[i*n],n);
    }
    Transpose column rows for
    (int i=0;i<n*n;i++)
    {
        int x = i%n;
        int y = i/n;
        int index = x*n+y;
        Source_2d[index] = Source_2d_1[i];
    }
    Vertical FFT for
    (int i=0;i<n;i++)
    {
        fft1 (&source_2d[i*n],&source_2d_1[i*n],n);
    }
    Transpose back for
    (int i=0;i<n*n;i++)
    {
        int x = i%n;
        int y = i/n;
        int index = x*n+y;
        Source_2d[index] = Source_2d_1[i];
    }

The GPU implementation simply transforms these things onto the GPU.

I calculate the FFT based on OpenGL's fragment shader, and the data is stored in textures or FBO. Unlike the 1-D FFT, the NXN data contains only one-dimensional FFT for the current column or the current row, so the bit reversal table requires only a 1*n buffer. The corresponding butterfly chart data also only need 1*n. So we have the following allocations:

Static OFFBO _fbo_bitrev_table;
Static OFFBO _origin_butterfly_2d;

_fbo_bitrev_table.allocate (n,1,gl_rgba32f);
_origin_butterfly_2d.allocate (n,1,gl_rgba32f);

The first thing to do is to find the bit reversal table of length n, which requires only one time, so at the very beginning, the CPU is calculated:

    for (int i=0;i<n;i++)
    {
          _bitrev_index_2d.setcolor (I,0,offloatcolor (Bit_rev (i,n-1), 0,0,0));
    }

    _bitrev_index_2d.update ();

    The inverted index
    _fbo_bitrev_table.begin ();
    _bitrev_index_2d.draw (0,0,n,1);
    _fbo_bitrev_table.end ();

Then initialize the original butterfly chart, which is the same as the 1-D FFT, but with a different length:

for (int i=0;i<n;i++)
    {
        //Initialize two-dimensional butterfly chart
        if (i%2==0)
        {
            _data_2d.setcolor (I,0,offloatcolor (0.f,2. f,0,i+1));
        }
        else
        {
            _data_2d.setcolor (I,0,offloatcolor (1.f,2.f,0,i-1));
        }
        
    }

    _data_2d.update ();

    /////////////////
    //Initialize 2D Butterfly chart
    _weight_index_2d[0].begin ();
    _data_2d.draw (0,0,n,1);
    _weight_index_2d[0].end ();
    Back up the 2D initial butterfly chart for the next new calculation
    _origin_butterfly_2d.begin ();
    _data_2d.draw (0,0,n,1);
    _origin_butterfly_2d.end ();

Auxiliary functions:

    static unsigned int bit_rev (unsigned int v, unsigned int maxv)
    {
        unsigned int t = log (MAXV + 1)/log (2);
        unsigned int ret = 0;
        unsigned int s = 0x80000000>> (+);
        for (unsigned int i = 0; i < T; ++i)
        {
            unsigned int r = v& (s << i);
            RET |= (R << (t-i-1)) >> (i);    
        }
        return ret;
    }

    static void Bit_reverse_copy (RBVector2 src[], RBVector2 des[], int len)
    {for
        (int i = 0; i < len;i++)
        {
            Des[bit_rev (i, len-1)] = Src[i];
        }
    }

The following defines the functions that calculate the 2-D Ifft:

void gpufft::ifft_2d (offbo& in,offbo& out,int size);

Where in is the input, out is the output, size is n, by the time the initialization was passed in once, write here is to facilitate debugging temporarily change the size.

The IFFT itself has the same code as the above, and the content becomes a variety of shader calculations:

void gpufft::ifft_2d (offbo& in,offbo& out,int size) {//disables alpha blending, otherwise drawing to FBO will mix alpha, resulting in data loss OFDISABLEALPHABL

    Ending ();
    Horizontal FFT _weight_index_2d[_cur_2d].begin ();
    _origin_butterfly_2d.draw (0,0,n,1);

    
    _weight_index_2d[_cur_2d].end ();
    _fbo_in_bitreved_2d.begin ();
    _bit_reverse_shader_2d.begin ();
    _bit_reverse_shader_2d.setuniform3f ("Iresolution", n,n,0);
    _bit_reverse_shader_2d.setuniform1i ("n", N);
    _bit_reverse_shader_2d.setuniform1i ("dir", 1);
    _bit_reverse_shader_2d.setuniformtexture ("Tex_origin", In.gettexturereference (), 1);
    _bit_reverse_shader_2d.setuniformtexture ("Tex_bitreverse_table", _fbo_bitrev_table.gettexturereference (), 2);
    Ofrect (0,0,n,n);
    _bit_reverse_shader_2d.end ();

    _fbo_in_bitreved_2d.end ();
    Data _res_back_2d[_cur_2d].begin after rollover ();
    _fbo_in_bitreved_2d.draw (0,0,n,n);


    _res_back_2d[_cur_2d].end ();
        for (int i = 1;i<n;i*=2) {_res_back_2d[1-_cur_2d].begin (); Ofclear (0,0,0,0);
        _gpu_fft_shader_2d.begin (); _gpu_fft_shader_2d.setuniform1i (