Problem when performing IFFT for Tessendorf's ocean waves

Question

I am implementing Jerry Tessendorf's ocean waves as described in his paper in c++ and OpenGL.

I implemented two compute shaders, one for the h_tilde_0 and its conjugate, and one for the frequency function for the IFFT.

After performing the real to complex IFFT using GLFFT I obtain the following height map after normalizing the obtained values between 0 and 255:

The first 20 values of the height map before normalizing are:

• [0] = {float} -24667.7344
• [1] = {float} 23987.6855
• [2] = {float} -22951.1094
• [3] = {float} 22144.1914
• [4] = {float} -21596.2676
• [5] = {float} 22228.2578
• [6] = {float} -23422.1855
• [7] = {float} 23971.5117
• [8] = {float} -24944.8633
• [9] = {float} 26373.8594
• [10] = {float} -28082.5977
• [11] = {float} 29593.1387
• [12] = {float} -30536.7344
• [13] = {float} 31605.7578
• [14] = {float} -33134.7813
• [15] = {float} 34946.293
• [16] = {float} -36430.6445
• [17] = {float} 37443.6523
• [18] = {float} -38190.6758
• [19] = {float} 39565.5781
• [20] = {float} -41436.6563

So the values I get after the inverse fourier transform are with inverted signs in an interleaved fashion. And if I were to skip every odd value I get the following height map:

And skipping even values (same IFFT):

It looks like the previous image contains the height map and the opposite sign version of itself. Since it looks like the height field is somehow encoded I don't think that the error is on the generation of the spectrum. I think of two possibilities:

1. I have wrong coordinates in the frequency domain and somehow produce this result.
2. I am not understanding how GLFFT expects the input and how the output is given. I have been reading the glsl files in the GLFFT repository and I can't find any mistakes.

I understand that the question is very broad and I am maybe not giving enough information (please tell me if I should upload any other code of result). I have been trying to debug the problem but I am not making any progress.

Thanks in advance for any help!

Appendice:

How I am performing the IFFT:

std::shared_ptr<abstractions::SSBO> update_fft_texture(ssbo_pointer& h_k_t, int N){
    GLFFT::FFTOptions options;
    options.type.fp16 = false;
    options.type.output_fp16 = false;
    options.type.input_fp16 = false;
    GLFFT::GLContext context;

    GLFFT::FFT fft(&context, N, N, GLFFT::ComplexToReal, GLFFT::Inverse, GLFFT::SSBO, GLFFT::SSBO, std::make_shared<GLFFT::ProgramCache>(), options);

    GLuint output_texture, input_texture;

    std::shared_ptr<abstractions::SSBO> out(new abstractions::SSBO(nullptr, 4*N*N, GL_DYNAMIC_COPY));

    input_texture = (*h_k_t).GetRendererId();
    output_texture = (*out).GetRendererId();

    // Adapt raw GL types to types which GLContext uses internally.
    GLFFT::GLBuffer adaptor_output(output_texture);
    GLFFT::GLBuffer adaptor_input(input_texture);

    GLCall(
            {
                // Do the FFT
                GLFFT::CommandBuffer *cmd = context.request_command_buffer();
                fft.process(cmd, &adaptor_output, &adaptor_input);
                context.submit_command_buffer(cmd);
            }
            )


    GLCall(glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT));

    return out;
}

Frequency domain function where red channel represents the real part and green the imaginary (it is not normalized and the values exceed the 0-255 range):

I omitted the h_tidle_0 and frequency shaders' code since it looks to me that the problem is not there.

Two relevant links:

This implementation that uses GLFFT: https://arm-software.github.io/opengl-es-sdk-for-android/ocean_f_f_t.html

This paper that obtains spectrum and frequency textures similar to mines: https://tore.tuhh.de/bitstream/11420/1439/1/GPGPU_FFT_Ocean_Simulation.pdf

Answer 1

I found the problem. Maybe this answer helps someone since there isn't much discussion on the topic.

I was doing the following change of coordinates in the frequency domain to fit the values on the 256x256 matrix:

(n, m) = (x-N/2, y-N/2)

where N, n and m appear in the ocean waves' paper and x, y are the coordinates of the matrix.

Then I modified the change of coordinates as is done in this implementation:

(n, m) = (alias(x), alias(y))

where

int alias(int a){
    if(a > N/2) a -= N;
    return N;
}

Then immediately got the correct result (different seed):