[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[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
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Help-gsl] Help with fft...,
Bill Oliver <=