
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.

>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

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

>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:
> 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)

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