[Help-gsl] Help with fft...

[Bill Oliver]

I am playing with the fft routines and am having a bizarre result.  The
forward and inverse ffts seem to work just fine, and adding and
subtracting ffts works OK.  However, when I attempt a convolution by
multiplying two ffts, I get a reasonable result *except* that the right
and left halves of the inverse are shifted.  In other words, for a
gaussian, instead of


                      11
                    1$##$1
                   1$#@@#$1
                  1$#@@@@$$1
                   1$#@@#$1
                    1$##$1
                      11

I get:
                      1
                    1$#
                  1$#@@1
                 1$#@@@#$1
                  1$#@@@#$1
                   1$#@@@#$1
                     address@hidden
                      1#$1
                       1


The shifting is only vertical.  It's not a matter of the actual complex
multiplication; printing out the values of each multiplication result
seems OK.  It's not a matter of mirroring the quadrants incorrectly,
since the inverse works from the straight fft or addition or
subtraction.

Any pointers to what I am doing wrong would be appreciated.  The boring
code is below:

**********************************************************


Here's the code I've wrapped it in -- part of my personal image
processing library:

For doing the fft:

( There is a forward_fft(image<T>& inimage) method that is
just a loop that calls the following method for each
channel (color) )


#define REAL(x,i) ( (x) [2*(i)] )
#define IMAG(x,i) ( (x) [2*(i) + 1] )


template <class T>
image<mycomplex>
fft_factory<T>::forward_fft(image<T>& inimage, int in_channel){

// get the dimensions and number of colors (channels)
// in the image

// dimsize is a vector holding the number of
// pixels in each dimension
vector <int> local_dimsize = inimage.show_dimsize();

// ndims is the number of dimensions in the image
int ndims = inimage.show_ndims();

// channels is the number of generalized colors
int nchannels = inimage.show_nchannels();

// get the number of points in the image
int npts = inimage.data.show_npoints();

// create the data array for the fft
double *data = new double[2*npts];

gsl_fft_complex_wavetable *wavetable = NULL;
gsl_fft_complex_workspace *workspace = NULL;

// Make the complex version of the real image.
// Once I understand how this is working, I'll just
// use the real fft routines.
for (int i=0;i<npts;i++){
        REAL(data,i) = (double)(inimage.data.show(in_channel,i));
        IMAG(data,i) = (double)0.0;
        }

// make a result image
// "mycomplex" is a wrapper for the stl complex<double>
image<mycomplex> resultimage(ndims, local_dimsize, nchannels);

wavetable = gsl_fft_complex_wavetable_alloc(npts);

workspace = gsl_fft_complex_workspace_alloc(npts);

gsl_fft_complex_forward(data,1,npts,wavetable,workspace);

gsl_fft_complex_wavetable_free(wavetable);
gsl_fft_complex_workspace_free(workspace);

mycomplex tmpdat;
for (int i=0;i<npts;i++){
        tmpdat = mycomplex(((double)REAL(data,i)), ((double)IMAG(data,i)));
        resultimage.data.set( in_channel,i,tmpdat );
        }

delete [] data;
return resultimage;
}

The inverse is:

( There is an inverse_fft(image<T>& inimage) method that is
just a loop that calls the following method for each
channel (color) )

template <class T>
image<T>
fft_factory<T>::inverse_fft(image<mycomplex>& inimage, int in_channel){

// get the dimensions, size and number of colors
// in the image
vector <int> local_dimsize = inimage.show_dimsize();
int ndims = inimage.show_ndims();
int nchannels = inimage.show_nchannels();

// make the result image
image<T> resultimage(ndims, local_dimsize, nchannels);

// get the number of points
int npts = resultimage.data.show_npoints();

double *data = new double[2*npts];


// transfer the dagta from the image to the data array
for (int i=0;i<npts;i++){
        REAL(data,i) = (double)(inimage.data.show(in_channel,i)).real();
        IMAG(data,i) = (double)(inimage.data.show(in_channel,i)).imag();
        }

gsl_fft_complex_wavetable *wavetable = NULL;
gsl_fft_complex_workspace *workspace = NULL;
wavetable = gsl_fft_complex_wavetable_alloc(npts);
workspace = gsl_fft_complex_workspace_alloc(npts);

gsl_fft_complex_inverse(data,1,npts,wavetable,workspace);

gsl_fft_complex_wavetable_free(wavetable);
gsl_fft_complex_workspace_free(workspace);

// transfer back to my image class
T tmp_dat;
for (int i=0;i<npts;i++){
        tmp_dat = (T)(REAL(data,i)/npts);
        resultimage.data.set( in_channel,i,tmp_dat);
        }

delete [] data;
return resultimage;
}

The convolution is:

template<class T>
image<T>
convolve_factory<T>::frequency_convolve(image<T>& inimage, image<T>& 
kernelimage){

        fft_factory<T> newfft;

        // make sure the image and kernel are compatible
        int inimage_ndims = inimage.show_ndims();
        vector<int> inimage_dimsize = inimage.show_dimsize();
        int inimage_nchannels = inimage.show_nchannels();

        int kernel_ndims = kernelimage.show_ndims();
        vector<int> kernel_dimsize = kernelimage.show_dimsize();
        int kernel_nchannels = kernelimage.show_nchannels();

        assert(inimage_ndims == kernel_ndims);
        for (int i=0;i<inimage_ndims;i++)
                assert(inimage_dimsize[i] == kernel_dimsize[i]);
        assert(inimage_nchannels == kernel_nchannels);

        // in the gnu scientific library, it looks like doing a
        // forward -> inverse fft will result
        // in a flipped image.  This is for reversing this with flip()
        vector<int> flipdims(2);
        flipdims[0] = 1;
        flipdims[1] = 1;

        // do the forward ffts
        image<mycomplex> in_forward = newfft.forward_fft(inimage);
        image<mycomplex> ke_forward = newfft.forward_fft(kernelimage);

        // multiply the images
        // as an aside, the image multiplication is pretty well
        // validated.
        image<mycomplex> tmp_result = in_forward*ke_forward;

        // do the inverse fft
        image<T> result = newfft.inverse_fft(tmp_result);

        // make it pretty
        result = result.mirror_quads();
        result = result.flip(flipdims);

        return result;
}


Thanks, and I'm sorry about all the code.

billo