Element-wise multiplication results between 2D arrays in KISSFFT are different than SciPy FFT

I'm experimenting with KISSFFT in C++ after being discouraged to use FFTPACK to process 2D arrays.

I wrote an element-wise multiplication function to multiply two 2D arrays after they've been transformed with kiss_fftnd(). The result of the multiplication is then transformed back with the inverse FFT function. Unfortunately, the results I get from kissfft in C are different from what I get with SciPy in python as you can see in the image below:

To test the multiplication function, after the 2D input array is transformed I multiply it with itself for simplicity purposes. Here is a simplified version in Python to illustrate the algorithm:

import numpy as np
from scipy import fft as scipy_fft

in1 = np.array([[  98,  92], \
                [   9,  21], \
                [ 130,   4]], dtype=np.uint8)

fft_out = scipy_fft.rfft2(in1)
fft_mult = fft_out * fft_out
ifft_data = scipy_fft.irfft2(fft_mult, in1.shape)
print('\nSciPy IRFFT2: shape=', ifft_data.shape, 'dtype=', ifft_data.dtype, '\n', ifft_data)

I can't think of a reason why this simple operation couldn't be done with kissfft, which means that my approach to multiplication is probably wrong. As the output of kiss_fftnd() is a 1D array and not 2D, maybe the logic I'm using to iterate on this array and perform element-wise multiplication is incorrect.

Why are these results different and how do I make kissfft return the same values as SciPy?

If you know a function in kissfft that already does the multiplication correctly, that would work for me too. Please don't suggest other libraries to do this job. I'm looking for an answer that specifically deals with kissfft.

This is the full source code in Python:

import numpy as np
from scipy import fft as scipy_fft

# complex_mult: multiplies two complex numbers
def complex_mult(n1, n2):
     real_part = n1.real*n2.real - n1.imag*n2.imag
     imag_part = n1.real*n2.imag + n2.real*n1.imag
     return complex(real_part, imag_part)

# fft2d_mult: multiplies two 2D arrays of complex numbers
def fft2d_mult(array1, array2):
    array_mult = np.empty(array1.shape, dtype=array1.dtype)
    h, w = in1.shape
    for j in range(h):
        for i in range(w):
            array_mult[j,i] = complex_mult(array1[j,i], array2[j,i])
    return array_mult


print("\n######################## SCIPY RFFT/MULT/IRFFT #######################\n");

# initialize input data
in1 = np.array([[  98,  92], \
                [   9,  21], \
                [ 130,   4]], dtype=np.uint8)

print('Original data: shape=', in1.shape, 'dtype=', in1.dtype, '\n', in1)

# perform 2D RFFT
fft_out = scipy_fft.rfft2(in1)
print('\nSciPy RFFT2: shape=', fft_out.shape, 'dtype=', fft_out.dtype, '\n', fft_out)

# perform element-wise multiplication
fft_mult = fft2d_mult(fft_out, fft_out) # equivalent to: fft_mult = fft_out * fft_out
print('\nMultiplication result: shape=', fft_mult.shape, 'dtype=', fft_mult.dtype, '\n', fft_mult)

# perform inverse 2D RFFT
ifft_data = scipy_fft.irfft2(fft_mult, in1.shape)
print('\nSciPy IRFFT2: shape=', ifft_data.shape, 'dtype=', ifft_data.dtype, '\n', ifft_data)

This is the full source code in C++:

// compile with: g++ so_issue.cpp -o so_issue -I kissfft kissfft/kiss_fft.c kissfft/tools/kiss_fftnd.c
#include "kissfft/kiss_fft.h"
#include "kissfft/tools/kiss_fftnd.h"

// fft2d: receives a 2D array of floats, performs the forward transform with kiss_fftnd() and converts it into a kiss_fft_cpx array
kiss_fft_cpx* fft2d(float* input, int width, int height)
{
    const int numDim = 2;
    int shape[numDim] = { width, height };
    int nfft = width * height;

    // allocate 2D input array for FFT
    kiss_fft_cpx* cin = new kiss_fft_cpx[nfft];
    memset(cin, 0, nfft * sizeof(kiss_fft_cpx));

    // allocate 2D output array for FFT
    kiss_fft_cpx* cout = new kiss_fft_cpx[nfft];
    memset(cout, 0, nfft * sizeof(kiss_fft_cpx));

    // copy the input data to cin
    int k = 0;
    int idx = 0;
    for (int j = 0; j < height; ++j)
    {
        for (int i = 0; i < width; ++i)
        {
            idx = i + width * j; // access 1D array as 2D
            cin[k++].r = input[idx];
        }
    }

    // execute 2D FFT
    bool inverse_fft = false;
    kiss_fftnd_cfg cfg_f = kiss_fftnd_alloc(shape, numDim, inverse_fft, 0, 0);
    kiss_fftnd(cfg_f, cin , cout);

    // release resources
    kiss_fft_free(cfg_f);
    delete[] cin;

    return cout;
}

// fft2d: receives an array of kiss_fft_cpx elements, performs the inverse transform with kiss_fftnd() and returns the result in a new kiss_fft_cpx array
kiss_fft_cpx* ifft2d(kiss_fft_cpx* input, int width, int height)
{
    const int numDim = 2;
    int shape[numDim] = { width, height };
    int nfft = width * height;

    // allocate 2D output array for FFT
    kiss_fft_cpx* cout = new kiss_fft_cpx[nfft];
    memset(cout, 0, nfft * sizeof(kiss_fft_cpx));

    // execute inverse 2D FFT
    bool inverse_fft = true;
    kiss_fftnd_cfg cfg_i = kiss_fftnd_alloc(shape, numDim, inverse_fft, 0, 0);
    kiss_fftnd(cfg_i, input , cout);

    // release resources
    kiss_fft_free(cfg_i);

    return cout;
}

// complex_mult: performs element-wise multiplication between two complex numbers
kiss_fft_cpx complex_mult(const kiss_fft_cpx& a, const kiss_fft_cpx& b)
{
    kiss_fft_cpx c;

    // real_part = a.real*b.real - a.imag*b.imag
    c.r = a.r*b.r - a.i*b.i;

    // imag_part = a.real*b.imag + b.real*a.imag
    c.i = a.r*b.i + b.r*a.i;

    return c;
}

// complex_mult: performs element-wise multiplication between two kiss_fft_cpx arrays
kiss_fft_cpx* fft2d_mult(kiss_fft_cpx* input1, kiss_fft_cpx* input2, int width, int height)
{
    int nfft = width * height;
    kiss_fft_cpx* output = new kiss_fft_cpx[nfft];
    memset(output, 0, nfft * sizeof(kiss_fft_cpx));

    int idx = 0;
    for (int j = 0; j < height; ++j)
    {
        for (int i = 0; i < width; ++i)
        {
            idx = i + width * j; // access 1D array as 2D
            output[idx] = complex_mult(input1[idx], input2[idx]);
        }
    }

    return output;
}

void run_test(float* in1, const int& w, const int& h)
{
printf("\n#######################  KISSFFT FFT/MULT/IFFT  #######################\n\n");

    printf("Original data:\n");
    int idx = 0;
    for (int j = 0; j < h; ++j)
    {
        for (int i = 0; i < w; ++i)
        {
            idx = i + w * j;
            printf("%.4f  \t", in1[idx]);
        }
        printf("\n");
    }

    /* perform FFT */

    kiss_fft_cpx* cout = fft2d((float*)in1, w, h);

    printf("\nkissfft FFT2D:\n");
    for (int j = 0; j < h; ++j)
    {
        for (int i = 0; i < w; ++i)
        {
            idx = i + w * j;
            printf("%.4f %.4fj  \t", cout[idx].r,  cout[idx].i);
        }
        printf("\n");
    }

    /* perform element-wise multiplication */

    kiss_fft_cpx* cout_mult = fft2d_mult(cout, cout, w, h);

    printf("\nMultiplication result:\n");
    for (int j = 0; j < h; ++j)
    {
        for (int i = 0; i < w; ++i)
        {
            idx = i + w * j;
            printf("%.4f %.4fj  \t", cout_mult[idx].r,  cout_mult[idx].i);
        }
        printf("\n");
    }

    /* perform inverse FFT */

    kiss_fft_cpx* cinput = ifft2d(cout_mult, w, h);

    printf("\nkissfft IFFT2D:\n");

    int nfft = w * h;
    for (int j = 0; j < h; ++j)
    {
        for (int i = 0; i < w; ++i)
        {
            idx = i + w * j;
            printf("%.4f  \t", cinput[idx].r / nfft); // div by N to scale data back to the original range
        }
        printf("\n");
    }

    // release resources
    delete[] cout_mult;
    delete[] cinput;
    delete[] cout;
}

int main()
{
    int h = 3,  w = 2;
    float in1[h][w] =
    {
        {  98,  92 },
        {   9,  21 },
        { 130,  4  }
    };

    run_test((float*)in1, w, h);

    return 0;
}

The problem was the order in which width and height were used in shape. This variable is later passed to kiss_fftnd_alloc() as an argument and height must be defined first:

const int numDim = 2;
int shape[numDim] = { height, width };

After making this change inside fft2d() and ifft2d() the application displayed the correct results.