Вычисление FFT и IFFT с библиотекой FFTW C ++

Question

Я пытаюсь вычислить БПФ и затем БПФ, просто чтобы проверить, смогу ли я получить тот же сигнал обратно, но я не совсем уверен, как это сделать.Вот как я делаю БПФ:

    plan = fftw_plan_r2r_1d(blockSize, datas, out, FFTW_R2HC, FFTW_ESTIMATE);
    fftw_execute(plan);

Answer 1

Вот пример. Это делает две вещи. Сначала он подготавливает входной массив in[N] в виде косинусной волны, частота которого равна 3, а величина равна 1,0, и Фурье преобразует ее. Итак, на выходе вы должны увидеть пик на out[3], а другой на out[N-3]. Поскольку величина волны косинуса равна 1,0, вы получите N / 2 при out[3] и out[N-3].

Во-вторых, обратное преобразование Фурье превращает массив out[N] обратно в in2[N]. И после правильной нормализации вы можете увидеть, что in2[N] идентичен in[N].

#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#define N 16
int main(void) {
  fftw_complex in[N], out[N], in2[N]; /* double [2] */
  fftw_plan p, q;
  int i;

  /* prepare a cosine wave */
  for (i = 0; i < N; i++) {
    in[i][0] = cos(3 * 2*M_PI*i/N);
    in[i][1] = 0;
  }

  /* forward Fourier transform, save the result in 'out' */
  p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute(p);
  for (i = 0; i < N; i++)
    printf("freq: %3d %+9.5f %+9.5f I\n", i, out[i][0], out[i][1]);
  fftw_destroy_plan(p);

  /* backward Fourier transform, save the result in 'in2' */
  printf("\nInverse transform:\n");
  q = fftw_plan_dft_1d(N, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(q);
  /* normalize */
  for (i = 0; i < N; i++) {
    in2[i][0] *= 1./N;
    in2[i][1] *= 1./N;
  }
  for (i = 0; i < N; i++)
    printf("recover: %3d %+9.5f %+9.5f I vs. %+9.5f %+9.5f I\n",
        i, in[i][0], in[i][1], in2[i][0], in2[i][1]);
  fftw_destroy_plan(q);

  fftw_cleanup();
  return 0;
}

Это вывод:

freq:   0  -0.00000  +0.00000 I
freq:   1  +0.00000  +0.00000 I
freq:   2  -0.00000  +0.00000 I
freq:   3  +8.00000  -0.00000 I
freq:   4  +0.00000  +0.00000 I
freq:   5  -0.00000  +0.00000 I
freq:   6  +0.00000  -0.00000 I
freq:   7  -0.00000  +0.00000 I
freq:   8  +0.00000  +0.00000 I
freq:   9  -0.00000  -0.00000 I
freq:  10  +0.00000  +0.00000 I
freq:  11  -0.00000  -0.00000 I
freq:  12  +0.00000  -0.00000 I
freq:  13  +8.00000  +0.00000 I
freq:  14  -0.00000  -0.00000 I
freq:  15  +0.00000  -0.00000 I

Inverse transform:
recover:   0  +1.00000  +0.00000 I vs.  +1.00000  +0.00000 I
recover:   1  +0.38268  +0.00000 I vs.  +0.38268  +0.00000 I
recover:   2  -0.70711  +0.00000 I vs.  -0.70711  +0.00000 I
recover:   3  -0.92388  +0.00000 I vs.  -0.92388  +0.00000 I
recover:   4  -0.00000  +0.00000 I vs.  -0.00000  +0.00000 I
recover:   5  +0.92388  +0.00000 I vs.  +0.92388  +0.00000 I
recover:   6  +0.70711  +0.00000 I vs.  +0.70711  +0.00000 I
recover:   7  -0.38268  +0.00000 I vs.  -0.38268  +0.00000 I
recover:   8  -1.00000  +0.00000 I vs.  -1.00000  +0.00000 I
recover:   9  -0.38268  +0.00000 I vs.  -0.38268  +0.00000 I
recover:  10  +0.70711  +0.00000 I vs.  +0.70711  +0.00000 I
recover:  11  +0.92388  +0.00000 I vs.  +0.92388  +0.00000 I
recover:  12  +0.00000  +0.00000 I vs.  +0.00000  +0.00000 I
recover:  13  -0.92388  +0.00000 I vs.  -0.92388  +0.00000 I
recover:  14  -0.70711  +0.00000 I vs.  -0.70711  +0.00000 I
recover:  15  +0.38268  +0.00000 I vs.  +0.38268  +0.00000 I

Answer 2

Вы хотя бы пытались прочитать более чем приличную документацию?

У них есть целое учебное пособие для вас, чтобы вы могли познакомиться с FFTW:

http://fftw.org/fftw3_doc/Tutorial.html#Tutorial

ОБНОВЛЕНИЕ: Полагаю, вы знаете, как работать с массивами Си, поскольку именно это используется для ввода и вывода.

Эта страница содержит необходимую информациюдля БПФ против IFFT (см. Аргументы-> знак).В документах также говорится, что input-> FFT-> IFFT-> n * input.Поэтому вам нужно быть осторожным, чтобы правильно масштабировать данные.

Answer 3

Вот что я сделал.Шахта разработана специально, чтобы дать тот же результат, что и Matlab.Это работает только для матрицы столбцов (вы можете заменить AMatrix на std::vector).

AMatrix<complex<double>> FastFourierTransform::Forward1d(const AMatrix<double> &inMat)
{
    AMatrix<complex<double>> outMat{ inMat.NumRows(), inMat.NumCols() };
    size_t N = inMat.NumElements();
    bool isOdd = N % 2 == 1;
    size_t outSize = (isOdd) ? ceil(N / 2 + 1) : N / 2;
    fftw_plan plan = fftw_plan_dft_r2c_1d(
        inMat.NumRows(),
        reinterpret_cast<double*>(&inMat.mutData()[0]), // mutData here is as same v.data() for std::vector
        reinterpret_cast<fftw_complex*>(&outMat.mutData()[0]),
        FFTW_ESTIMATE);
    fftw_execute(plan);

    // mirror
    auto halfWayPoint = outMat.begin() + (outMat.NumRows() / 2) + 1;
    auto startCopyDest = (isOdd) ? halfWayPoint : halfWayPoint - 1;
    std::reverse_copy(outMat.begin() + 1, halfWayPoint, startCopyDest); // skip DC (+1)
    std::for_each(halfWayPoint, outMat.end(), [](auto &c) { return c = std::conj(c);});

    return std::move(outMat);
}

AMatrix<complex<double>> FastFourierTransform::Forward1d(const AMatrix<double> &inMat, size_t points)
{
    // append zeros to matrix until it's the same size as points
    AMatrix<double> sig = inMat;
    sig.Resize(points, sig.NumCols());
    for (size_t i = inMat.NumRows(); i < points; i++)
    {
        sig(i, 0) = 0;
    }
    return Forward1d(sig);
}

AMatrix<complex<double>> FastFourierTransform::Inverse1d(const AMatrix<complex<double>> &inMat)
{
    AMatrix<complex<double>> outMat{ inMat.NumRows(), inMat.NumCols() };
    size_t N = inMat.NumElements();

    fftw_plan plan = fftw_plan_dft_1d(
        inMat.NumRows(), 
        reinterpret_cast<fftw_complex*>(&inMat.mutData()[0]),
        reinterpret_cast<fftw_complex*>(&outMat.mutData()[0]),
        FFTW_BACKWARD, 
        FFTW_ESTIMATE);
    fftw_execute(plan);

    // Matlab normalizes
    auto normalize = [=](auto &c) { return c *= 1. / N; };
    std::for_each(outMat.begin(), outMat.end(), normalize);

    fftw_destroy_plan(plan);
    return std::move(outMat);
}

// From Matlab documentation: "ifft tests X to see whether vectors in X along the active dimension 
// are conjugate symmetric. If so, the computation is faster and the output is real. 
// An N-element vector x is conjugate symmetric if x(i) = conj(x(mod(N-i+1,N)+1)) for each element of x."
// http://uk.mathworks.com/help/matlab/ref/ifft.html
AMatrix<double> FastFourierTransform::Inverse1dConjSym(const AMatrix<complex<double>> &inMat)
{
    AMatrix<double> outMat{ inMat.NumRows(), inMat.NumCols() };
    size_t N = inMat.NumElements();

    fftw_plan plan = fftw_plan_dft_c2r_1d(
        inMat.NumRows(), 
        reinterpret_cast<fftw_complex*>(&inMat.mutData()[0]),
        reinterpret_cast<double*>(&outMat.mutData()[0]),
        FFTW_BACKWARD | FFTW_ESTIMATE);
    fftw_execute(plan);

    auto normalize = [=](auto &c) { return c *= 1. / N; };
    std::for_each(outMat.begin(), outMat.end(), normalize);

    fftw_destroy_plan(plan);
    return std::move(outMat);
}

Answer 4

Очень полезно для увеличения линии изображения

Как и в следующем примере, увеличить в 2 раза прямой сигнал

int main (void)
{
  fftw_complex in[N], out[N], in2[N2], out2[N2];        /* double [2] */
  fftw_plan p, q;
  int i;
  int half;

  half=(N/2+1);
  /* prepare a cosine wave */
  for (i = 0; i < N; i++)
    {
      //in[i][0] = cos ( 3 * 2 * M_PI * i / N);
      in[i][0] = (i > 3 && i< 12 )?1:0;
      in[i][1] = (i > 3 && i< 12 )?1:0;
      //in[i][1] = 0;
    }

  /* forward Fourier transform, save the result in 'out' */
  p = fftw_plan_dft_1d (N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute (p);
  for (i = 0; i < N; i++)
    printf ("input: %3d %+9.5f %+9.5f I   %+9.5f %+9.5f I\n", i, in[i][0], in[i][1],out[i][0],out[i][1]);
  fftw_destroy_plan (p);

  for (i = 0; i<N2; i++) {out2[i][0]=0.;out2[i][1]=0.;}

  for (i = 0; i<half; i++) {
    out2[i][0]=2.*out[i][0];
    out2[i][1]=2.*out[i][1];
  }
  for (i = half;i<N;i++) {
    out2[N+i][0]=2.*out[i][0];
    out2[N+i][1]=2.*out[i][1];
  }



  /* backward Fourier transform, save the result in 'in2' */
  printf ("\nInverse transform:\n");
  q = fftw_plan_dft_1d (N2, out2, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute (q);
  /* normalize */
  for (i = 0; i < N2; i++)
    {
      in2[i][0] /= N2;
      in2[i][1] /= N2;
    }
  for (i = 0; i < N2; i++)
    printf ("recover: %3d %+9.1f %+9.1f I\n",
            i, in2[i][0], in2[i][1]);
  fftw_destroy_plan (q);

  fftw_cleanup ();
  return 0;
}