C++ FFTW 前向后向 DFT 值被包裹

Question

你好 StackOverflow 社区， 我对 fftw 库的 dft 算法有疑问。 我想做的就是将某个模式向前和向后转换以再次接收输入模式，当然稍后在转换之间会有某种过滤。

所以，我的程序在 atm 上做的是:

1. 创建测试信号
2. 用 1.0 或 0.5 的值过滤或“窗口化”测试信号
3. 将测试信号复制到 fftw_complex 数据类型
4. 执行正向和反向 dft
5. 计算幅度，这里称为相位
6. 复制和调整数据以供显示，最后通过OpenCV显示图像

我的问题是，当不使用过滤时，我的反向变换图像以某种方式被包裹，我无法计算正确的幅度，这应该与我的输入图像/测试信号相同。 当我将 fitler/“window”的值设置为 0.5 时，向后转换工作正常，但我的输入图像只有应有的亮度的一半。

下图说明了我的问题:(从左上到右下) 1. 输入信号， 2. 反向变换的实部， 3. 从反向变换数据计算幅度， 4. 输入信号乘以 0.5， 5. 反向变换实部， 6. 从反向变换数据计算幅度。 http://imageshack.com/a/img538/5426/nbL9YZ.png

有没有人知道为什么 dft 以这种方式执行？!有点奇怪...

我的代码看起来像这样的 atm:

/***** parameters **************************************************************************/
int     imSize                                          = 256;
int     imN                                             = imSize * imSize;

char*   interferogram                                   = new char[imN];
double* spectrumReal                                    = new double[imN];
double* spectrumImaginary                               = new double[imN];
double* outputReal                                      = new double[imN];
double* outputImaginary                                 = new double[imN];
double* phase                                           = new double[imN];

char*   spectrumRealChar                                = new char[imN];
char*   spectrumImaginaryChar                           = new char[imN];
char*   outputRealChar                                  = new char[imN];
char*   outputImaginaryChar                             = new char[imN];
char*   phaseChar                                       = new char[imN];

Mat     interferogramMat                                = Mat(imSize, imSize, CV_8U, interferogram);
Mat     spectrumRealCharMat                             = Mat(imSize, imSize, CV_8U, spectrumRealChar);
Mat     spectrumImaginaryCharMat                        = Mat(imSize, imSize, CV_8U, spectrumImaginaryChar);
Mat     outputRealCharMat                               = Mat(imSize, imSize, CV_8U, outputRealChar);
Mat     outputImaginaryCharMat                          = Mat(imSize, imSize, CV_8U, outputImaginaryChar);
Mat     phaseCharMat                                    = Mat(imSize, imSize, CV_8U, phaseChar);


/***** compute interferogram ****************************************************************/
fill_n(interferogram, imN, 0);
double value = 0;
double window = 0;

for (int y = 0; y < imSize; y++)
{
    for (int x = 0; x < imSize; x++)
    {
        value = 127.5 + 127.5 * cos((2*PI) / 10000 * (pow(double(x - imSize/2), 2) + pow(double(y - imSize/2), 2)));

        window = 1;
        value *= window;

        interferogram[y * imSize + x] = (unsigned char)value;
    }
}


/***** create fftw arays and plans **********************************************************/
fftw_complex*       input;
fftw_complex*       spectrum;
fftw_complex*       output;
fftw_plan           p_fw;
fftw_plan           p_bw;

input               = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
spectrum            = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
output              = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
p_fw                = fftw_plan_dft_2d(imSize, imSize, input, spectrum, FFTW_FORWARD, FFTW_ESTIMATE);
p_bw                = fftw_plan_dft_2d(imSize, imSize, spectrum, output, FFTW_BACKWARD, FFTW_ESTIMATE);


/***** copy data ****************************************************************************/
for (int i = 0; i < imN; i++)
{
    input[i][0] = double(interferogram[i]) / 255.;
    input[i][1] = 0.;
    spectrum[i][0] = 0.;
    spectrum[i][1] = 0.;
    output[i][0] = 0.;
    output[i][1] = 0.;
}


/***** FPS algorithm ************************************************************************/
fftw_execute(p_fw);

fftw_execute(p_bw);

for (int i = 0; i < imN; i++)
{
    phase[i] = sqrt(pow(output[i][0], 2) + pow(output[i][1], 2));
}


/***** copy data ****************************************************************************/
for (int i = 0; i < imN; i++)
{
    spectrumReal[i] = spectrum[i][0];
    spectrumImaginary[i] = spectrum[i][1];

    outputReal[i] = output[i][0] / imN;
    outputImaginary[i] = output[i][1];
}

SaveCharImage(interferogram, imN, "01_interferogram_512px_8bit.raw");
SaveDoubleImage(spectrumReal, imN, "02_spectrum_real_512px_64bit.raw");
SaveDoubleImage(spectrumImaginary, imN, "03_spectrum_imaginary_512px_64bit.raw");
SaveDoubleImage(outputReal, imN, "03_output_real_512px_64bit.raw");

DoubleToCharArray(spectrumReal, spectrumRealChar, imSize);
DoubleToCharArray(spectrumImaginary, spectrumImaginaryChar, imSize);

DoubleToCharArray(outputReal, outputRealChar, imSize);
DoubleToCharArray(outputImaginary, outputImaginaryChar, imSize);

DoubleToCharArray(phase, phaseChar, imSize);


/***** show images **************************************************************************/

imshow("interferogram", interferogramMat);
imshow("spectrum real", spectrumRealCharMat);
imshow("spectrum imaginary", spectrumImaginaryCharMat);
imshow("out real", outputRealCharMat);
imshow("out imaginary", outputImaginaryCharMat);
imshow("phase", phaseCharMat);

int key = waitKey(0);

Answer 1

以下是您的几行代码:

char*   interferogram  = new char[imN];
...
double value = 0;
double window = 0;

for (int y = 0; y < imSize; y++)
{
   for (int x = 0; x < imSize; x++)
   {
      value = 127.5 + 127.5 * cos((2*PI) / 10000 * (pow(double(x - imSize/2), 2) + pow(double(y - imSize/2), 2)));

      window = 1;
      value *= window;

      interferogram[y * imSize + x] = (unsigned char)value;
   }
}

问题是 char 介于 -128 和 127 之间，而 unsigned char 介于 0 到 255 之间。在 interferogram[y * imSize + x ] = (unsigned char)value;，隐式转换为 char。

如果 window=0.5 不会影响输出，但如果 window=1 因为 value 变得高于 127，它会触发更改. 这正是您在问题中注意到的问题!

它不会影响第一个显示的图像，因为 CV_8U 对应于 unsigned char :interferogram 因此被转换回 unsigned字符*。看看Can I turn unsigned char into char and vice versa?了解有关 char 到 unsigned char 转换的更多信息。

问题出现在 input[i][0] = double(interferogram[i])/255.; :如果 window=1，interferogram [i] 可能为负数，input[i][0] 变为负数。

将所有的 char 更改为 unsigned char 应该可以解决问题。

你也可以改变

outputReal[i] = output[i][0] / imN;
outputImaginary[i] = output[i][1];

对于

outputReal[i] = output[i][0];
outputImaginary[i] = output[i][1];

调用 fftw 似乎没问题。