2012-03-02 64 views
1

我想讓FFTW在C中工作,它曾經爲另一個項目(這是在JNI中)工作,我或多或少地從該項目複製代碼,可悲的是沒有結果。FFTW結果爲零

首先,我產生一個正弦信號,就像這樣:

double* generateSignal() { 
    int fs=44100; 
    double fsd = 44100.0; // fs in double format 
    double f1=1000.0; 
    int i; 
    double PI = 3.141592653589793238462643; 

    double t[fs]; 
    double value = 0.0; 
    for (i = 0; i < fs; i++) { 
     t[i] = value; 
     value += 1.0/fsd; 
    } 

    double* signal = (double*) malloc(sizeof(double) * fs); 
    for (i = 0; i < fs; i++) { 
     signal[i] = sqrt(2) * sin(2 * PI * f1 * t[i]); 
    } 

    return signal; 
}  

這工作,我只張貼的完整性。

接下來,我想變換使用FFTW的信號,這是我在下面的方法(基於FFTW documentation)做:

void processSignal(double* signal) { 
    int size = 44100; 
    int i; 

    fftw_complex* in = fftw_malloc(sizeof(fftw_complex) * size); 
    fftw_complex* out = fftw_malloc(sizeof(fftw_complex) * size); 

    for (i = 0; i < size; i++) { 
     double* ptr = in[i]; 
     *ptr = signal[i];  // set first double, real part 
     *(ptr + 1) = 0.0;  // set second double, imaginary part 
    } 

    fftw_plan p = fftw_plan_dft_1d(size, in, out, FFTW_FORWARD, FFTW_ESTIMATE); 
    fftw_execute(p); 

    for (i = 0; i < size; i++) { 
     double* ptr = out[i]; 
     signal[i] = *ptr;  // get real part 
    } 

    fftw_destroy_plan(p); 
    fftw_free(in); 
    fftw_free(out); 
} 

請從FFTW文檔注意這一點:typedef double fftw_complex[2];

現在,這導致signal-陣列的所有值爲-0.000000。我真的不能看到這個代碼的問題,所以請指出我做錯了什麼?

謝謝。

PS:爲了清晰起見,從我的代碼中刪除了打印語句。

+0

您是否收到任何警告/錯誤?這行不應該編譯:'double * ptr = in [i];' – Mysticial 2012-03-02 15:47:52

+0

沒有任何,甚至沒有設置-Wall標誌(我意識到,這很奇怪)。我看到你的觀點,爲什麼它不會編譯,但如果是這種情況,我會通過將它轉換爲雙指針來解決它。無論如何,我想這不是導致問題的原因。 – gleerman 2012-03-02 15:54:50

+0

您沒有投射fftw_malloc(sizeof(fftw_complex)* size)的輸出; – L7ColWinters 2012-03-02 15:56:28

回答

0

天哪,這是愚蠢的。

這些值具有-2.0E-9的大小。結果是四捨五入的,因此打印了-0.00000。

非常感謝Fred。他的問題是「你的預期產量是多少?」使我對前面提到的JNI項目的輸出有了很好的瞭解,並意識到了我的錯誤。

也感謝您的其他答案!

0

我認爲這個問題可能是你用前一個迭代的虛部覆蓋了當前實部的實部。

for (i = 0; i < size; i++) { 
     double* ptr = in[i]; // <-- here's a problem 
     *ptr = signal[i];  // set first double, real part 
     *(ptr + 1) = 0.0;  // set second double, imaginary part 
    } 

i被遞增逐一在每次迭代,所以在第一次迭代中,ptr點的輸入複數[0]和ptr + 1實部是指向的複數[0]的虛部,但在第二次迭代時,ptr指向複數[0]的虛數部分,而ptr + 1指向複數[1]的實數部分。

一些建議,以解決此問題可能是:

for (i = 0, j = 0; i < size; i++, j+= 2) { 
    double* ptr = in[j]; // j increments by 2 making ptr alays point to real part 
    *ptr = signal[i];  // set first double, real part 
    *(ptr + 1) = 0.0;  // set second double, imaginary part 
} 

double* ptr = in[0] 
for (i = 0; i < size; i++) { 
     *ptr++ = signal[i]; // set first double, real part 
     *ptr++ = 0.0;  // set second double, imaginary part 
    }