Результаты поэлементного умножения между двумерными массивами в KISSFFT отличаются от SciPy FFT.

Я экспериментирую с KISSFFT в C++ после того, как мне не рекомендовали использовать FFTPACK для обработки двумерных массивов.

Я написал функцию поэлементного умножения для умножения двух двумерных массивов после их преобразования с помощью kiss_fftnd(). Затем результат умножения преобразуется обратно с помощью функции обратного БПФ. К сожалению, результаты, которые я получаю с kissfft на C, отличаются от тех, которые я получаю с SciPy на python, как вы можете увидеть на изображении ниже:

Чтобы проверить функцию умножения, после преобразования двумерного входного массива я умножаю его сам на себя для простоты. Вот упрощенная версия на Python для иллюстрации алгоритма:

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)

Я не могу придумать причину, по которой эту простую операцию нельзя было выполнить с помощью kissfft, а это значит, что мой подход к умножению, вероятно, неверен. Поскольку вывод kiss_fftnd() представляет собой массив 1D, а не 2D, возможно, логика, которую я использую для итерации этого массива и выполнения поэлементного умножения, неверна.

Почему эти результаты отличаются и как сделать, чтобы kissfft возвращал те же значения, что и SciPy?

Если вы знаете функцию в kissfft, которая уже правильно выполняет умножение, она мне тоже подойдет. Пожалуйста, не предлагайте другие библиотеки для выполнения этой работы. Я ищу ответ, который конкретно касается kissfft.

Это полный исходный код на 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)

Это полный исходный код на 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;
}