Forums

FFTW?

Started by tcharles December 14, 2006
Hello,
can anyone explain me how to perform the convolution of two dimentional
array with fftw?
Thanks

tcharles wrote:
> Hello, > can anyone explain me how to perform the convolution of two dimentional > array with fftw? > Thanks
This method of performing convolution is a relatively straightforward application of the FFT; see these links: http://www.dspguide.com/ch18/2.htm http://cnx.org/content/m10963/latest/ This would be no different with FFTW than any other environment; just the mechanics of the source code used to perform the FFT would be different. Jason
tcharles wrote:
> can anyone explain me how to perform the convolution of two dimentional > array with fftw?
Sometime back I wrote a FFTW based plugin for Image Apprentice. You may consult that. Here's the link: http://www.adislindia.com/pdn.html or alternately (direct link) http://divyarathore.googlepages.com/fftwpluginforimageapprentice.htm I didn't do convolution there but it is self explainatory. In case you still need help, just write back. regards, Divya Rathore www.adislindia.com (remove underscores for email ID)
>Sometime back I wrote a FFTW based plugin for Image Apprentice. You may >consult that. Here's the link: >http://www.adislindia.com/pdn.html >or alternately (direct link) >http://divyarathore.googlepages.com/fftwpluginforimageapprentice.htm > >I didn't do convolution there but it is self explainatory. >In case you still need help, just write back. > >regards, >Divya Rathore >www.adislindia.com >(remove underscores for email ID) >
Thank you for your answer. I still have a problem with the fftw, to perform the convolution I wrote the following functions: fftw_complex * fft::forward(double * in) { forward_result = (fftw_complex *) fftw_malloc ( sizeof ( fftw_complex ) * nx *(ny/2+1) ); plan_forward = fftw_plan_dft_r2c_2d( nx, ny, in, forward_result, FFTW_ESTIMATE); fftw_execute ( plan_forward ); return(forward_result); } double * fft::backward(fftw_complex * in) { out = new double[nx*(ny/2+1)]; /* Backward transformation */ plan_backward = fftw_plan_dft_c2r_2d(nx, ny, in, out, FFTW_ESTIMATE); fftw_execute ( plan_backward ); return(out); } //double * fft::multiply(double * A, double * B) { fftw_complex * fft::multiply(fftw_complex * A, fftw_complex * B) { double scale = 1.0 / (nx * ny); R = (fftw_complex *) fftw_malloc ( sizeof ( fftw_complex ) * nx *(ny/2+1) ); for(int i= 0; i < nx; i++) for(int j= 0; j < ny/2+1; j++){ int ij = i*(ny/2+1) + j; R[ij][0] = (A[ij][0] * B[ij][0] - A[ij][1] * B[ij][1]) * scale; R[ij][1] = (A[ij][0] * B[ij][0] + A[ij][1] * B[ij][0]) * scale; } return(R); } and I have a segmentation fault in the backward function, the debugger output is at 0x82175DA: apply_hc2r (in /home/tchoutat/Work/projekt_tp1/nina/fdp/fdp) ==27739== by 0x821A687: apply (in /home/tchoutat/Work/projekt_tp1/nina/fdp/fdp) ==27739== Address 0x459FD68 is 0 bytes after a block of size 10,608 alloc'd ==27739== at 0x401AD0D: operator new[](unsigned) (vg_replace_malloc.c:195) ==27739== by 0x8096B14: fft::backward(double (*) [2]) (fft_l.C:66) And I don't understand this size problem. Thanks for your help. Best regards. Charles Tchoutat
tcharles wrote:
> double * fft::backward(fftw_complex * in) { > out = new double[nx*(ny/2+1)]; > /* Backward transformation */ > plan_backward = fftw_plan_dft_c2r_2d(nx, ny, in, out, FFTW_ESTIMATE); > fftw_execute ( plan_backward ); > return(out); > }
> and I have a segmentation fault in the backward function, the debugger
The real output array of an out-of-place c2r transform should be of size nx*ny, not nx*(ny/2+1). You are confusing the size of the complex array with the size of the real array. See the FFTW manual. You also have a memory leak because you don't call fftw_destroy_plan. Note also that using "new" instead of "fftw_malloc" may prevent FFTW from exploiting SIMD instructions (e.g. SSE/SSE2), since "new" is probably not 16-byte aligned. Steven G. Johnson
>tcharles wrote: >> double * fft::backward(fftw_complex * in) { >> out = new double[nx*(ny/2+1)]; >> /* Backward transformation */ >> plan_backward = fftw_plan_dft_c2r_2d(nx, ny, in, out,
FFTW_ESTIMATE);
>> fftw_execute ( plan_backward ); >> return(out); >> } > >> and I have a segmentation fault in the backward function, the debugger > >The real output array of an out-of-place c2r transform should be of >size nx*ny, not nx*(ny/2+1). You are confusing the size of the complex >array with the size of the real array. See the FFTW manual. > >You also have a memory leak because you don't call fftw_destroy_plan. > >Note also that using "new" instead of "fftw_malloc" may prevent FFTW >from exploiting SIMD instructions (e.g. SSE/SSE2), since "new" is >probably not 16-byte aligned. > >Steven G. Johnson
Thanks for your answer, I wrote a destructor to call fftw_destroy_plan und fftw_free, und made the following changes: fft::~fft() { fftw_destroy_plan (plan_backward ); fftw_destroy_plan (plan_forward ); fftw_free ( out ); fftw_free ( R ); } double * fft::backward(fftw_complex * in) { out = (double *)fftw_malloc ( sizeof ( double ) *nx*ny ); /* Backward transformation */ plan_backward = fftw_plan_dft_c2r_2d(nx, ny, in, out, FFTW_ESTIMATE); fftw_execute ( plan_backward ); return(out); } but I still have size fault. The debugger output is: Address 0x464C028 is 0 bytes inside a block of size 20,808 alloc'd ==2095== at 0x401A4CE: malloc (vg_replace_malloc.c:149) ==2095== by 0x8096AF4: fft::backward(double (*) [2]) (fft_l.C:68) the line 68 is out = (double *)fftw_malloc ( sizeof ( double ) *nx*ny ); Thanks for your help. Charles Tchoutat

On Dec 15 2006, 5:48 pm, "divya_ratho...@gmail.com" 
<divyarath...@gmail.com> wrote:
Sometime back I wrote a FFTW based plugin for Image Apprentice. You 
may
> consult that. Here's the link:http://www.adislindia.com/pdn.html > or alternately (direct link)http://divyarathore.googlepages.com/fftwpluginforimageapprentice.htm > > I didn't do convolution there but it is self explainatory. > In case you still need help, just write back. > > regards, > Divya Rathorewww.adislindia.com > (remove underscores for email ID)
A plugin alongwith source code has just been published. It does Gaussian High Pass Filtering using convolution (utilizes FFTW) http://www.adislindia.com/login/pluginlinks.php One may use it as a reference.