2011-08-17 71 views
1

我想將Matlab的快速傅里葉變換函數fft()轉換爲本地Java代碼。將Matlab的FFT轉換爲本地Java

由於我使用其中FFT作爲實施的JMathLib代碼起點如下:

// given double[] x as the input signal 

    n = x.length; // assume n is a power of 2 
    nu = (int)(Math.log(n)/Math.log(2)); 
    int n2 = n/2; 
    int nu1 = nu - 1; 
    double[] xre = new double[n]; 
    double[] xim = new double[n]; 
    double[] mag = new double[n2]; 
    double tr, ti, p, arg, c, s; 
    for (int i = 0; i < n; i++) { 
     xre[i] = x[i]; 
     xim[i] = 0.0; 
    } 
    int k = 0; 

    for (int l = 1; l <= nu; l++) { 
     while (k < n) { 
      for (int i = 1; i <= n2; i++) { 
       p = bitrev (k >> nu1); 
       arg = 2 * (double) Math.PI * p/n; 
       c = (double) Math.cos (arg); 
       s = (double) Math.sin (arg); 
       tr = xre[k+n2]*c + xim[k+n2]*s; 
       ti = xim[k+n2]*c - xre[k+n2]*s; 
       xre[k+n2] = xre[k] - tr; 
       xim[k+n2] = xim[k] - ti; 
       xre[k] += tr; 
       xim[k] += ti; 
       k++; 
      } 
      k += n2; 
     } 
     k = 0; 
     nu1--; 
     n2 = n2/2; 
    } 
    k = 0; 
    int r; 
    while (k < n) { 
     r = bitrev (k); 
     if (r > k) { 
      tr = xre[k]; 
      ti = xim[k]; 
      xre[k] = xre[r]; 
      xim[k] = xim[r]; 
      xre[r] = tr; 
      xim[r] = ti; 
     } 
     k++; 
    } 
    // The result 
    // -> real part stored in xre 
    // -> imaginary part stored in xim 

不幸的是,不給我正確的結果,當我單元測試,例如與array

double [] x = {1.0d,5.0d,9.0d,13.0d};

結果在Matlab:

28.0
-8.0 - 8.0i的
-8.0
-8.0 + 8.0i的

結果在我的實現:

28.0
-8.0 + 8.0i的
-8.0
-8.0 - 8.0i的

注意的跡象是如何錯在複雜的部分。

當我使用更長,更復雜的信號時,實現之間的差異也會影響數字。所以實現差異不僅與某些符號 - 「錯誤」有關。

我的問題:我如何調整自己的實現以使其與Matlab的「相等」? 或者:是否已經有一個圖書館完全做到這一點?

+2

如果我在Matlab中運行這個,我得到第二組數字。 –

+2

不是'本地Java'是一個矛盾嗎? :-) – Praetorian

+2

你能舉出一個8號尺寸差異的例子嗎?大小16?謝謝。 – Nayuki

回答

1

爲了在矩陣上使用Jtransforms進行FFT,您需要按col執行fft col,然後將它們連接成矩陣。這裏是我的代碼,我與Matlab比較fft

double [][] newRes = new double[samplesPerWindow*2][Matrixres.numberOfSegments]; 
double [] colForFFT = new double [samplesPerWindow*2]; 
DoubleFFT_1D fft = new DoubleFFT_1D(samplesPerWindow); 

for(int y = 0; y < Matrixres.numberOfSegments; y++) 
{ 
    //copy the original col into a col and and a col of zeros before FFT 
    for(int x = 0; x < samplesPerWindow; x++) 
    { 
     colForFFT[x] = Matrixres.res[x][y];   
    } 

    //fft on each col of the matrix 
    fft.realForwardFull(colForFFT);              //Y=fft(y,nfft); 

    //copy the output of col*2 size into a new matrix 
    for(int x = 0; x < samplesPerWindow*2; x++) 
    { 
     newRes[x][y] = colForFFT[x];  
    } 

} 

希望這是你在找什麼。注意Jtransforms表示複數爲

array[2*k] = Re[k], array[2*k+1] = Im[k]