Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to do inverse real to real FFT in FFTW library

Tags:

c++

c

fft

inverse

fftw

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);

like image 641
Kirill Dubovikov Avatar asked Dec 12 '10 15:12

Kirill Dubovikov


2 Answers

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.)

like image 200
wnoise Avatar answered Nov 08 '22 23:11

wnoise


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;
}
like image 31
hbp Avatar answered Nov 09 '22 01:11

hbp