Informática FFT e IFFT con biblioteca FFTW C++

Question

Estoy tratando de calcular la FFT y luego la IFFT solo para probar si puedo obtener la misma señal pero no estoy seguro de cómo lograrla. Así es como hago la FFT:

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

Answer 1

¿Al menos ha intentado leer la documentación más que decente?

Ellos tienen todo un tutorial pensado para que usted pueda llegar a conocer FFTW:

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

ACTUALIZACIÓN: Asumo que sabe cómo trabajar con C-arrays, ya que eso es lo que se utiliza como entrada y salida.

This page tiene la información que necesita para FFT vs IFFT (ver Argumentos-> signo). Los documentos también dicen que entrada-> FFT-> IFFT-> n * entrada. Por lo tanto, deberá tener cuidado de escalar los datos correctamente.

Answer 2

Aquí hay un ejemplo. Hace dos cosas. Primero, prepara una matriz de entrada in[N] como una onda de coseno, cuya frecuencia es 3 y su magnitud es 1.0, y Fourier la transforma. Por lo tanto, en la salida, debería ver un pico en out[3] y otra en out[N-3]. Como la magnitud de la onda del coseno es 1.0, obtienes N/2 en out[3] y out[N-3].

En segundo lugar, el inverso de Fourier transforma la matriz out[N] de nuevo en in2[N]. Y después de una normalización adecuada, puede ver que in2[N] es idéntico a 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; 
} 

Esta es la salida:

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 3

Muy útil para la línea magnify2x de imagen

como en el siguiente ejemplo engrandezcas por el factor de señal 2 rect

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

Answer 4

Aquí es lo Hice. El mío está diseñado específicamente para dar el mismo resultado que Matlab. Esto solo funciona en una matriz de columna (puede reemplazar AMatrix con un 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); 
}