KISSFFT中2D阵列之间的逐元素乘法结果与SciPy FFT不同 [英] Element-wise multiplication results between 2D arrays in KISSFFT are different than SciPy FFT
问题描述
不鼓励使用 C ++ 中的 KISSFFT 使用 FFTPACK处理2D阵列.
我编写了一个逐元素乘法函数，用于将两个二维数组用kiss_fftnd()转换后相乘.然后，通过逆FFT函数将乘法结果转换回去.不幸的是，我在 C 中从 kissfft 获得的结果与您在 python 中从 SciPy 获得的结果不同可以在下图中看到:
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:
为了测试乘法功能，在转换2D输入数组后，为了简单起见，我将其与自身相乘.这是Python的简化版本，用于说明算法:
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)
我想不出使用 kissfft 无法完成此简单操作的原因，这意味着我的乘法方法可能是错误的.由于kiss_fftnd()的输出是一个1D数组而不是2D数组，所以也许我在该数组上进行迭代并执行逐元素乘法的逻辑是错误的.
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.
为什么这些结果不同，如何使Kissfft返回与SciPy相同的值?
如果您知道 kissfft 中的一个函数，该函数已经正确执行了乘法运算，那么对我也同样适用.请不要建议其他库来完成此工作.我正在寻找专门针对 kissfft 的答案.
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.
这是Python中的完整源代码:
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)
这是C ++中的完整源代码:
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;
}
推荐答案
问题是在shape中使用width和height的顺序.此变量随后作为参数传递给kiss_fftnd_alloc()，并且必须首先定义height:
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 };
在fft2d()和ifft2d()中进行此更改之后，该应用程序显示了正确的结果.
After making this change inside fft2d() and ifft2d() the application displayed the correct results.
这篇关于KISSFFT中2D阵列之间的逐元素乘法结果与SciPy FFT不同的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
