Как сделать обратное реальное к реальному FFT в библиотеке FFTW - PullRequest
Новая
       20

Как сделать обратное реальное к реальному FFT в библиотеке FFTW

5 голосов
/ 12 декабря 2010

Я пытаюсь сделать фильтрацию с помощью БПФ. Я использую план r2r_1d, и я не знаю, как сделать обратное преобразование ... 

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

}

Я все делаю правильно? Я в замешательстве, потому что в сложных БПФ FFTW вы можете установить направление преобразования с помощью флага, например так:
p = fftw_plan_dft_1d (N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
или
p = fftw_plan_dft_1d (N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);

Ответы [ 2 ]

4 голосов
/ 13 декабря 2010

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

«Вид» указывает направление.

(Обратите внимание, что вы, вероятно, захотите перенормировать ваш сигнал путем деления на n. Соглашение о нормализации FFTW умножается на n после преобразования и его инверсии.)

1 голос
/ 26 февраля 2014

Вы сделали это правильно.FFTW_REDFT00 означает косинусное преобразование, которое является его собственным обратным.Поэтому нет необходимости различать «вперед» и «назад».Однако будьте осторожны с размером массива.Если вы хотите определить частоту 10, и ваши данные содержат 100 значащих точек, то массив data должен содержать 101 точек данных и установить n = 101 вместо 100. Нормализация должна быть 2*(n-1).См. Пример ниже, скомпилируйте с 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;
}
...