FFTW transform of a real and symmetric image is not real and symmetric?

Question

I have an issue regarding the FFT implementation in FFTW. Or maybe the problem is my knowledge on FFT. The point is, as far as I know, if I transform forward a symmetric and real signal, I should get a symmetric and real signal also. However, this is not what I am finding in FFTW. I am working in 2d, with the fftwf_plan_dft_r2c_2d routine.

I leave here a simple piece of code that creates a matrix that stores in each pixel (x,y) the result of computing 1/(sqrt(x^2+y^2)), which is a symmetric signal. There is also a loop that prints out these values, just to check that I built the matrix properly.

#include <fftw3.h>

int main(int argc, char* argv[]){

width=5;height=5;
float *omega;
omega = (float *) malloc(sizeof(float)*width*height);
float *omegasym = (float *) malloc(sizeof(float)*2*width*2*height);
int y,  x;

//Build a symmetric real signal
for(y = 0, i = 0; y < height; y++) /* omega = 1/sqrt(x^2 + y^2) */
    for(x = 0; x < width; x++, i++)
        omega[i] = (x == 0 && y == 0) ? 0 : 1.0f/(float) std::sqrt((float) x*x + y*y);

//Just check if we did well 
for(int i=0; i<5; i++){
    for(int j=0; j<5; j++){
        std::cout<<" " << omega[j+5*i] << " ";
        }
        std::cout<<endl;
    }


fftwf_complex *omega_trans_complexx;
omega_trans_complexx = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*2*height*(2*width/2+1));


fftwf_plan plan;
plan=fftwf_plan_dft_r2c_2d(height, width, omega, omega_trans_complexx, FFTW_ESTIMATE);

//Should not this has imaginary part = 0??? :S
for(int i=0; i<25; i++){
    std::cout<<"Real part of omega_trans is :" << omega_trans_complexx[i][0] << endl;
    std::cout<<"Imaginary part of omega_trans is :" << omega_trans_complexx[i][1] << endl;
}

return 0;
}

After computing the forward transform with FFTW, i am getting a non-zero imaginary part. At first I thought it could be the fact that the signal is not symmetric, but for a symmetrized version of it, it does not work, neither. I wonder if I am doing something wrong, or FFTW requires some kind of zero-padding, or I am missing something else..

Answer 1

This question seems to be partly about fft's properties

x,y-> 1/(sqrt(x^2+y^2)) is neither even nor odd. It is normal that you get a non-zero imaginary part. You could expect a zero imaginary part if you try :

 omega[i]=sin(2.0*M_PI/height*y)*sin(2.0*M_PI/width*x);

or :

 omega[i]=cos(2.0*M_PI/height*y)*cos(2.0*M_PI/width*x);

And something purely imaginary if you try :

 omega[i]=cos(2.0*M_PI/height*y)*sin(2.0*M_PI/width*x);

The effect of using real as input is that you don't have to compute the fft for negative frequencies since it is the conjugate of the corresponding positive frequency : \hat{f}(-k)=(\hat{f}(k))*

On the code itself :

• You won't get anything useful if you don't call fftwf_execute(plan); somewhere !
• sizeof(fftwf_complex)*2*height*(2*width/2+1) is too large. You can save some memory by using sizeof(fftwf_complex)*2*height*(width/2+1) instead.
• you are right about padding. In fact, fftw stores an "extra" frequency in the x direction. That why you should allocate width/2+1 complex numbers in the x direction. As soon as you correct the printing operations, you will get symetric values between x and y

Here is the result :

#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <stdlib.h>

int main(int argc, char* argv[]){

    int width=5;int height=5;
    float *omega;
    omega = (float *) fftwf_malloc(sizeof(float)*width*height);
    float *omegasym = (float *) malloc(sizeof(float)*2*width*2*height);
    int y,  x;

    //Build a symmetric real signal
    int i;
    for(y = 0, i = 0; y < height; y++) /* omega = 1/sqrt(x^2 + y^2) */
        for(x = 0; x < width; x++, i++)
            //omega[i]=sin(2.0*M_PI/height*y)*cos(2.0*M_PI/width*x);
            omega[i] = (x == 0 && y == 0) ? 0 : 1.0f/(float) std::sqrt((float) x*x + y*y);

    //Just check if we did well 
    for(i=0; i<5; i++){
        for(int j=0; j<5; j++){
            std::cout<<" " << omega[j+5*i] << " ";
        }
        std::cout<<std::endl;
    }


    fftwf_complex *omega_trans_complexx;
    omega_trans_complexx = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*2*height*(width/2+1));


    fftwf_plan plan;
    plan=fftwf_plan_dft_r2c_2d(height, width, omega, omega_trans_complexx, FFTW_ESTIMATE);

    fftwf_execute(plan);

    //Should not this has imaginary part = 0??? :S : imaginary part is ok :)
    for(y = 0, i = 0; y < height; y++) /* omega = 1/sqrt(x^2 + y^2) */
        for(x = 0; x < width/2+1; x++, i++)
            std::cout<<"freqx "<<x<<" freqy "<<y<<" real "<<omega_trans_complexx[i][0]<< " imag "<<omega_trans_complexx[i][1]<<std::endl;

    return 0;
}

Bye,

Francis