Hello,

I just download the fftw3 library and read the tutorial.

In fact, my goal is to make the convolution of two huge double arrays, and I want to use the Fourier transform for its processing speed.

I know that is necessary to follow the algorithm below to do convolution with FFT:

Step 1

TF (A) = FFT (A) with one of our arrays - size(A)= M

TF (B) = FFT (B) B with one of our arrays - size(A)= N

Step 2

Then doing TF (A) * TF (B)

Step 3

And eventually be reversed by getting Conv (A, B) = IFFT (TF (A) * TF (B)) which has a size equal to M + N-1.

My problem is in step 2, I really do not know how to implement this step.

I know there must be a zero padding of the vectors A and B having a size equal to M + N-1 before doing their Fourier transform

Help me, it's the only step that remains for me to finish a project.

I put my code below and line is the function that enable the multiplication of the Fourier transform (help me on it)

[code]

void confftW (fftw_complex * A * B fftw_complex, int M, int N) {

fftw_complex Apadding *, * Bpadding, ApaddingFFT, BpaddingFFT * R * IR;

fftw_plan plan_a, plan_b, plan_R;

int i;

/ * Allocate input & output array * /

Apadding = (fftw_complex *) fftw_malloc (sizeof (fftw_complex) * M + N-1);

Bpadding = (fftw_complex *) fftw_malloc (sizeof (fftw_complex) * M + N-1);

ApaddingFFT = (fftw_complex *) fftw_malloc (sizeof (fftw_complex) * M + N-1);

BpaddingFFT = (fftw_complex *) fftw_malloc (sizeof (fftw_complex) * M + N-1);

IR = (fftw_complex *) fftw_malloc (sizeof (fftw_complex) * M + N-1);

/ / 0-padding for A and B

for (i = 0; i <M + N-1, i + +)

{

/ / For simplicity it is assumed that A and B are the same size

/ / M = N

if (i <M)

{

Apadding [i] [0] = A [i] [0];

Bpadding [i] [0] = B [i] [0];

}

else

{

Apadding [i] [0] = 0.0;

Bpadding [i] [0] = 0.0;

}

Apadding [i] [1] = 0.0;

Bpadding [i] [1] = 0.0;

}

/ * Create plans * /

plan_a fftw_plan_dft_1d = (M + N-1, Apadding, ApaddingFFT, FFTW_FORWARD, FFTW_ESTIMATE);

plan_b fftw_plan_dft_1d = (M + N-1, Bpadding, BpaddingFFT, FFTW_FORWARD, FFTW_ESTIMATE);

plan_R fftw_plan_dft_1d = (M + N-1, R, IR, FFTW_BACKWARD, FFTW_ESTIMATE);

/ * Implementation of TF A and B * /

fftw_execute (plan_a);

fftw_execute (plan_b);

R = Multiply (BpaddingFFT, ApaddingFFT) / / FUNCTION THAT I WANT TO IMPLEMENT

/ / Run the inverse Fourier transform of R, which will store

/ / The result in IR

fftw_execute (plan_R);

/ * Free memory * /

fftw_destroy_plan (plan_a);

fftw_destroy_plan (plan_b);

fftw_destroy_plan (plan_R);

fftw_free (Apadding);

fftw_free (Bpadding);

fftw_free (ApaddingFFT);

fftw_free (BpaddingFFT);

fftw_free (IR);

}

[/ code]

Thank you in advance