How do I use inverse FFT from FFTPack?

Question

I'm trying to use FFT from the ANSI C (C89) compatible FFT library FFTPack. So far, I can transform time domin data into frequecy domain data

/* Fast Fourier Transform 
 * x[n] - Will contains the FFT values. The first value will have no imaginary part
 */
void fft(float x[], size_t n) {

    /* Init */
    float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
    int* ifac = (int*)malloc((n + 15) * sizeof(float));
    __ogg_fdrffti(n, wsave, ifac);

    /* Forward transform */
    __ogg_fdrfftf(n, x, wsave, ifac);

    /* Free */
    free(wsave);
    free(ifac);
}

The output float x[]

0.5010999       0.4075914       0.4787459       -0.5674312      -3.2158968      -2.5260105      2.5675464       -2.4416344      -1.6614010      1.3101853       3.0937805

Compare with the output from Octave/MATLAB

>> x
x =

  -0.6485   0.2851   1.3475  -0.6743  -0.1499  -1.5549   1.4951   0.4504  -0.4982   0.3488   0.1000

>> fft(x)
ans =

 Columns 1 through 5:

   0.5011 +      0i   0.4076 + 0.4787i  -0.5674 - 3.2159i  -2.5260 + 2.5675i  -2.4416 - 1.6614i

 Columns 6 through 10:

   1.3102 + 3.0938i   1.3102 - 3.0938i  -2.4416 + 1.6614i  -2.5260 - 2.5675i  -0.5674 + 3.2159i

 Column 11:

   0.4076 - 0.4787i

>>

Now I'm trying backward FFT, so called Inverse FFT.

/* Fast Fourier Transform
 * x[n] - Will contains the inverse FFT values. The first value will have no imaginary part
 */
void ifft(float x[], size_t n) {

    /* Init */
    float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
    int* ifac = (int*)malloc((n + 15) * sizeof(float));
    __ogg_fdrffti(n, wsave, ifac);

    /* Backward transform */
    __ogg_fdrfftb(n, x, wsave, ifac);

    /* Free */
    free(wsave);
    free(ifac);
}

The output float x[]

-2.9383001      -4.9931164      1.6534307       -0.0710702      -7.5080967      -4.3042197      4.2668128       -0.5331497      -0.7820582      6.1745853       1.9016848

Compare with the output from Octave/MATLAB

>> x
x =

  -0.6485   0.2851   1.3475  -0.6743  -0.1499  -1.5549   1.4951   0.4504  -0.4982   0.3488   0.1000

>> ifft(x)
ans =

 Columns 1 through 5:

   0.0456 +      0i   0.0371 - 0.0435i  -0.0516 + 0.2924i  -0.2296 - 0.2334i  -0.2220 + 0.1510i

 Columns 6 through 10:

   0.1191 - 0.2813i   0.1191 + 0.2813i  -0.2220 - 0.1510i  -0.2296 + 0.2334i  -0.0516 - 0.2924i

 Column 11:

   0.0371 + 0.0435i
>>

So the regular FFT function __ogg_fdrfftf is working. But the inverse function __ogg_fdrfftb is not working as it should be.

Question:

Does anyone know how to use __ogg_fdrfftb inside FFTPack library?

How do you call these functions? This is critical. `n` for the inverse transform must have the same value as for the forward transform. It’s the size of the output. — Cris Luengo, Aug 27 '23 at 23:23
Oh, your `x` in the frequency domain is a real-valued half-transform, you cannot do `ifft(x)` on that in MATLAB and get the right result. — Cris Luengo, Aug 27 '23 at 23:24

euraad · Answer 1 · 2023-08-28T14:16:03.313

1

I found the answer.

/* Fast Fourier Transform
 * x[n] - Will contains the inverse FFT values. The first value will have no imaginary part
 */
void ifft(float x[], size_t n) {

    /* Init */
    float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
    int* ifac = (int*)malloc((n + 15) * sizeof(float));
    __ogg_fdrffti(n, wsave, ifac);

    /* Backward transform */
    __ogg_fdrfftb(n, x, wsave, ifac);
    
    /* Normalize */
    size_t i;
    for (i = 0; i < n; i++) {
        x[i] /= n;
    }

    /* Free */
    free(wsave);
    free(ifac);
}

And the input x[] will be the same as output x[]

#define n 10

int main() {
    clock_t start, end;
    float cpu_time_used;
    start = clock();

    /* Create x */
    float x[n] = {1, 2, 5, 10, 5, 3.5, 1, 2, 3, 10};

    /* Do Inverse FFT */
    fft(x, n);
    ifft(x, n);

    /* Print */
    print(x, 1, n);

    end = clock();
    cpu_time_used = ((float)(end - start)) / CLOCKS_PER_SEC;
    printf("\nTotal speed  was %f\n", cpu_time_used);

    return EXIT_SUCCESS;
}

edited Aug 28 '23 at 14:16

answered Aug 28 '23 at 13:39

euraad

2,467
5
30
51

Instead of hard-coding `n 10`, you could do `float x[] = {...}` and `fft(x, sizeof x / sizeof *x)`; let the compiler do the counting. – Neil Aug 28 '23 at 14:24
@Neil Thank you. But someone need to tell the length of `x` at the beginning? – euraad Aug 28 '23 at 14:39
An _initialized_ array with empty brackets `x[] = ...` infers the size in bytes; this is all the way back to C89. Using `sizeof`, if you add elements, one doesn't have to change anything else. – Neil Aug 28 '23 at 14:58
@Neil Please, show me an example. – euraad Aug 28 '23 at 18:37

How do I use inverse FFT from FFTPack?

1 Answers1