用FFTW函數庫計算FFT和IFFT C++

Question

我想計算FFT然後IFFT只是爲了試驗如果我可以得到相同的信號但我不確定如何實現它。這是我如何做FFT：

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

Answer 1

你有沒有至少嘗試閱讀更多體面的文件？

他們有一個完整的教程想爲你去FFTW知道：

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

更新：我假設你知道如何與C-陣列工作，因爲那正是用作輸入並輸出。

This page有您需要的FFT與IFFT相關的信息（請參閱參數 - >符號）。文檔還說輸入 - > FFT-> IFFT-> n *輸入。所以你必須小心地正確擴展數據。

Answer 2

這裏是一個例子。它有兩件事。首先，準備一個輸入數組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 3

用於圖像的magnify2x線非常有用的

如在下面的例子中通過因子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; 
} 

Answer 4

這裏是我做了。我的設計專門爲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); 
}