3
今天我寫了一個算法來計算表示離散函數的給定數組點的快速傅里葉變換。現在我試圖測試它，看它是否工作。我已經嘗試過十幾個不同的輸入集，而且他們似乎與我在網上找到的示例相匹配。然而，對於我的最終測試，我給出了cos（i/2）的輸入，其中i從0到31，並且根據我使用的求解器得到3個不同的結果。我的解決方案似乎是最不準確的：驗證FFT算法的正確性
這是否表明我的算法有問題，或者是它只是相對較小的數據集的結果？
我的代碼如下，在情況下，它可以幫助：
/**
* 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;
}
}
}
第二張圖顯然是錯的;實際輸入應該導致對稱輸出。 –
如果你想驗證你的實現（儘管浮點數值問題），只需編寫一個天真的DFT（約五行代碼），並進行比較。 –