validationFFT算法的正确性

今天我写了一个算法来计算表示离散函数的给定点arrays的快速傅里叶变换。 现在我正在尝试测试它是否有效。 我已经尝试了十几个不同的输入集，它们似乎与我在网上找到的例子相匹配。 然而，对于我的最终测试，我给它输入cos（i / 2），i从0到31，并且根据我使用的求解器得到了3个不同的结果。 我的解决方案似乎是最不准确的：

这是否表明我的算法有问题，还是仅仅是相对较小的数据集的结果？

我的代码如下，如果它有帮助：

/** * Slices the original array, starting with start, grabbing every stride elements. * For example, slice(A, 3, 4, 5) would return elements 3, 8, 13, and 18 from array A. * @param array The array to be sliced * @param start The starting index * @param newLength The length of the final array * @param stride The spacing between elements to be selected * @return A sliced copy of the input array */ public double[] slice(double[] array, int start, int newLength, int stride) { double[] newArray = new double[newLength]; int count = 0; for (int i = start; count < newLength && i < array.length; i += stride) { newArray[count++] = array[i]; } return newArray; } /** * Calculates the fast fourier transform of the given function. The parameters are updated with the calculated values * To ignore all imaginary output, leave imaginary null * @param real An array representing the real part of a discrete-time function * @param imaginary An array representing the imaginary part of a discrete-time function * Pre: If imaginary is not null, the two arrays must be the same length, which must be a power of 2 */ public void fft(double[] real, double[] imaginary) throws IllegalArgumentException { if (real == null) { throw new NullPointerException("Real array cannot be null"); } int N = real.length; // Make sure the length is a power of 2 if ((Math.log(N) / Math.log(2)) % 1 != 0) { throw new IllegalArgumentException("The array length must be a power of 2"); } if (imaginary != null && imaginary.length != N) { throw new IllegalArgumentException("The two arrays must be the same length"); } if (N == 1) { return; } double[] even_re = slice(real, 0, N/2, 2); double[] odd_re = slice(real, 1, N/2, 2); double[] even_im = null; double[] odd_im = null; if (imaginary != null) { even_im = slice(imaginary, 0, N/2, 2); odd_im = slice(imaginary, 1, N/2, 2); } fft(even_re, even_im); fft(odd_re, odd_im); // F[k] = real[k] + imaginary[k] // even odd // F[k] = E[k] + O[k] * e^(-i*2*pi*k/N) // F[k + N/2] = E[k] - O[k] * e^(-i*2*pi*k/N) // Split complex arrays into component arrays: // E[k] = er[k] + i*ei[k] // O[k] = or[k] + i*oi[k] // e^ix = cos(x) + i*sin(x) // Let x = -2*pi*k/N // F[k] = er[k] + i*ei[k] + (or[k] + i*oi[k])(cos(x) + i*sin(x)) // = er[k] + i*ei[k] + or[k]cos(x) + i*or[k]sin(x) + i*oi[k]cos(x) - oi[k]sin(x) // = (er[k] + or[k]cos(x) - oi[k]sin(x)) + i*(ei[k] + or[k]sin(x) + oi[k]cos(x)) // { real } { imaginary } // F[k + N/2] = (er[k] - or[k]cos(x) + oi[k]sin(x)) + i*(ei[k] - or[k]sin(x) - oi[k]cos(x)) // { real } { imaginary } // Ignoring all imaginary parts (oi = 0): // F[k] = er[k] + or[k]cos(x) // F[k + N/2] = er[k] - or[k]cos(x) for (int k = 0; k < N/2; ++k) { double t = odd_re[k] * Math.cos(-2 * Math.PI * k/N); real[k] = even_re[k] + t; real[k + N/2] = even_re[k] - t; if (imaginary != null) { t = odd_im[k] * Math.sin(-2 * Math.PI * k/N); real[k] -= t; real[k + N/2] += t; double t1 = odd_re[k] * Math.sin(-2 * Math.PI * k/N); double t2 = odd_im[k] * Math.cos(-2 * Math.PI * k/N); imaginary[k] = even_im[k] + t1 + t2; imaginary[k + N/2] = even_im[k] - t1 - t2; } } } 

Answer 1

1. validation

   看看这里： 慢DFT，iDFT最后是我的DFT和iDFT的慢速实施它们经过测试和纠正。 我过去也使用它们进行快速实现validation。

2. 你的代码

   停止递归是错误的（你忘了设置返回元素）我看起来像这样：

    if (n<=1) { if (n==1) { dst[0]=src[0]*2.0; dst[1]=src[1]*2.0; } return; } 

   所以当你的N==1 ，输出元素设置为Re=2.0*real[0], Im=2.0*imaginary[0] return之前。 我在复杂的数学(t,t1,t2)有点迷失，也懒得分析。

我要确保这是我的快速实施。 它需要来自类层次结构的太多东西，所以它对你来说不会是另一种用途，然后与你的代码进行视觉比较。

我的快速实现（cc意味着复杂的输出和输入）：

 //--------------------------------------------------------------------------- void transform::DFFTcc(double *dst,double *src,int n) { if (n>N) init(n); if (n<=1) { if (n==1) { dst[0]=src[0]*2.0; dst[1]=src[1]*2.0; } return; } int i,j,n2=n>>1,q,dq=+N/n,mq=N-1; // reorder even,odd (buterfly) for (j=0,i=0;i 

dst[]和src[]不重叠!!! 所以你不能将数组转换为自身。
_cos和_sin是cos和sin值的预计算表（由init（）函数计算如下：

  double a,da; int i; da=2.0*M_PI/double(N); for (a=0.0,i=0;i 

N是2幂（数据集的零填充大小）（来自init(n)调用的最后n init(n) ）

只是要完成这里是我的复杂到复杂的慢版本：

 //--------------------------------------------------------------------------- void transform::DFTcc(double *dst,double *src,int n) { int i,j; double a,b,a0,a1,_n,b0,b1,q,qq,dq; dq=+2.0*M_PI/double(n); _n=2.0/double(n); for (q=0.0,j=0;j 

Answer 2

我会使用像Wolfram Alpha这样权威的东西进行validation。

如果我为0 <= i < 32评估cos(i/2) ，我得到这个数组：

 [1,0.878,0.540,0.071,-0.416,-0.801,-0.990,-0.936,-0.654,-0.211,0.284,0.709,0.960,0.977,0.754,0.347,-0.146,-0.602,-0.911,-0.997,-0.839,-0.476,0.004,0.483,0.844,0.998,0.907,0.595,0.137,-0.355,-0.760,-0.978] 

如果我把它作为Wolfram Alpha的FFT函数的输入，我会得到这个结果 。

我得到的情节看起来是对称的，这是有道理的。 情节看起来与您提供的情节完全不同。