help-gsl
[Top][All Lists]

## [Help-gsl] Help with fft...

 From: Bill Oliver Subject: [Help-gsl] Help with fft... Date: Wed, 8 Sep 2004 22:41:20 -0400 (EDT)

```
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

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

I get:
1
1\$#
1\$#@@1
1\$#@@@#\$1
1\$#@@@#\$1
1\$#@@@#\$1
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 = 1;
flipdims = 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.flip(flipdims);

return result;
}

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

billo

```