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 会这样吗？！有点奇怪……

我的代码看起来像这样：

/***** 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 char*。查看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 似乎没问题。