How to do inverse real to real FFT in FFTW library
Asked Answered
M

2

6

I'm trying to do some filtering with FFT. I'm using r2r_1d plan and I have no idea how to do the inverse transform...

    void PerformFiltering(double* data, int n)
    {
                    /* FFT */
        double* spectrum = new double[n];

        fftw_plan plan;

        plan = fftw_plan_r2r_1d(n, data, spectrum, FFTW_REDFT00, FFTW_ESTIMATE);

        fftw_execute(plan); // signal to spectrum
        fftw_destroy_plan(plan); 


                    /* some filtering here */


                    /* Inverse FFT */
        plan = fftw_plan_r2r_1d(n, spectrum, data, FFTW_REDFT00, FFTW_ESTIMATE);
        fftw_execute(plan); // spectrum to signal (inverse FFT)
        fftw_destroy_plan(plan);

}

Am I doing all the things right? I'm confused because in FFTW complex DFT's you can set a transform direction by flag like this:
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
or
p = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);

Mating answered 12/12, 2010 at 15:2 Comment(0)
S
4

http://www.fftw.org/fftw3_doc/Real_002dto_002dReal-Transform-Kinds.html

http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html

The "kind" specifies the direction.

(Note also that you'll probably want to renormalize your signal by dividing by n. The normalization convention of FFTW multiplies by n after a transform and its inverse.)

Siddra answered 12/12, 2010 at 22:46 Comment(0)
B
1

You did it correctly. FFTW_REDFT00 means the cosine transform, which is its own inverse. So there is no need to distinguish "forward" and "backward." However, be careful about the array size. If you want to detect a frequency of 10, and your data contains 100 meaningful points, then the array data should hold 101 data points, and set n = 101 instead of 100. The normalization should be 2*(n-1). See the example below, compile with gcc a.c -lfftw3.

#include <stdio.h>
#include <math.h>
#include <fftw3.h>
#define n 101 /* note the 1 */
int main(void) {
  double in[n], in2[n], out[n];
  fftw_plan p, q;
  int i;
  p = fftw_plan_r2r_1d(n, in, out, FFTW_REDFT00, FFTW_ESTIMATE);
  for (i = 0; i < n; i++) in[i] = cos(2*M_PI*10*i/(n - 1)); /* n - 1 instead of n */
  fftw_execute(p);
  q = fftw_plan_r2r_1d(n, out, in2, FFTW_REDFT00, FFTW_ESTIMATE);
  fftw_execute(q);
  for (i = 0; i < n; i++)
    printf("%3d %9.5f %9.5f\n", i, in[i], in2[i]/(2*(n - 1))); /* n - 1 instead of n */
  fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup();
  return 0;
}
Bartz answered 12/12, 2010 at 15:3 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.