# FFTW?

Started by 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

I didn't do convolution there but it is self explainatory.
In case you still need help, just write back.

regards,
Divya Rathore
(remove underscores for email ID)

```
```>Sometime back I wrote a FFTW based plugin for Image Apprentice. You may
>
>I didn't do convolution there but it is self explainatory.
>In case you still need help, just write back.
>
>regards,
>Divya Rathore
>(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] = (A[ij] * B[ij]
- A[ij] * B[ij]) * scale;
R[ij] = (A[ij] * B[ij]
+ A[ij] * B[ij]) * 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
(vg_replace_malloc.c:195)
==27739==    by 0x8096B14: fft::backward(double (*) ) (fft_l.C:66)
And I don't understand this size problem.
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 (*) ) (fft_l.C:68)
the line 68 is out = (double *)fftw_malloc ( sizeof ( double ) *nx*ny );
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
>
> I didn't do convolution there but it is self explainatory.
> In case you still need help, just write back.
>
> regards,