2011-08-17 70 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]