python - KISSFFT中2D阵列之间的逐元素乘法结果与SciPy FFT不同

不鼓励使用KISSFFT之后，我正在 C++ 中尝试FFTPACK to process 2D arrays。

我编写了一个逐元素乘法函数，用于将两个2D数组用kiss_fftnd()转换后相乘。然后，通过逆FFT函数将乘法的结果转换回去。不幸的是，我从 C 中的 kissfft 得到的结果与我在 python 中从 SciPy 得到的结果不同，如下图所示:



为了测试乘法功能，在转换完二维输入数组后，为了简单起见，我将其与自身相乘。这是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;
}

最佳答案

问题是在width中使用height和shape的顺序。稍后将此变量作为参数传递给kiss_fftnd_alloc()，并且必须首先定义height:

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

在fft2d()和ifft2d()中进行此更改之后，该应用程序显示了正确的结果。