Convolution using FFTW c++
CodeGuru Home VC++ / MFC / C++ .NET / C# Visual Basic VB Forums Developer.com
Results 1 to 2 of 2

Thread: Convolution using FFTW c++

  1. #1
    Join Date
    Oct 2011
    Posts
    0

    Convolution using FFTW c++

    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

  2. #2
    Join Date
    Oct 2005
    Location
    Minnesota, U.S.A.
    Posts
    680

    Re: Convolution using FFTW c++

    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

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  


Azure Activities Information Page

Windows Mobile Development Center


Click Here to Expand Forum to Full Width

This is a CodeGuru survey question.


Featured


HTML5 Development Center