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);
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.)
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;
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With