Convolution using FFTW c++

[pachalcs]

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

[egawtry]

You are using too many of your own defined routines for anyone to get a good handle on what you are doing.

For a good example of this take a look at "Digital Signal Processing" by Steven Smith, p. 233.

-Erik