java - DIT FFT Radix-2 算法故障排除

Question

我已经在 J​​ava 中实现了一个递归基数 2 DIT FFT，以及一个常规 DFT 来验证我的 FFT 结果，但两者的结果不同，我似乎无法弄清楚。两者都使用 apply() 方法提供整个数组，开始和停止索引分别为 0 和 data.length。 DFT 版本看起来正确，在 bin 50 处有一个漂亮的峰值，而 FFT 版本充满了垃圾。我做错了什么？

这是 FFT 实现(改编自 http://www.engineeringproductivitytools.com/stuff/T0001/PT04.HTM “A Recursive DIT FFT Routine.”，我通过与 https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#Pseudocode 的伪代码进行比较验证):

public class DITFFT2 extends Transform {
public float[] apply(float[] data, int startIndex, int stopIndex) throws IllegalArgumentException {
    int N;
    float[] filteredData;
    Complex[] complexData;
    Complex[] filteredComplexData;

    if (stopIndex < startIndex) {
        throw new IllegalArgumentException("stopIndex cannot be lower than startIndex!");
    }

    if (stopIndex < 0 || startIndex < 0) {
        throw new IllegalArgumentException("Index cannot be negative!");
    }

    N = stopIndex - startIndex;
    filteredData = new float[N];
    complexData = new Complex[N];

    for (int i = startIndex; i < stopIndex; i++) {
        complexData[i-startIndex] = new Complex(data[i], 0.0f);
    }

    filteredComplexData = transform(complexData, N);

    for (int i = 0; i < N; i++) {
        filteredData[i] = filteredComplexData[i].abs();
    }

    return filteredData;
}

public Complex[] transform(Complex[] data, int N) {
    Complex x;
    Complex[] result = new Complex[N];

    if (N == 1) {
        result[0] = data[0];
    } else {
        Complex[] fe = new Complex[N/2];
        Complex[] fo = new Complex[N/2];

        for (int i = 0; i < N/2; i++) {
            fe[i] = data[2*i];
            fo[i] = data[2*i+1];
        }

        Complex[] Fe = transform(fe, N / 2);
        Complex[] Fo = transform(fo, N / 2);

        for (int k = 0; k < N/2; k++) {
            x = Fo[k].copy();
            x.mul(getTwiddleFactor(k, N));

            result[k] = Fe[k].copy();
            result[k].add(x);

            result[k+N/2] = Fe[k].copy();
            result[k+N/2].sub(x);
        }
    }

    return result;
}

private Complex getTwiddleFactor(int k, int N) {
    return new Complex(1.0f, (float)(-2.0f * Math.PI * k / (float)N));
}
}

这是 DFT 实现:

public class DFT extends Transform {
public float[] apply(float[] data, int startIndex, int stopIndex) throws IllegalArgumentException {
    int N;
    float[] filteredData;
    Complex[] complexData;
    Complex[] filteredComplexData;

    if (stopIndex < startIndex) {
        throw new IllegalArgumentException("stopIndex cannot be lower than startIndex!");
    }

    if (stopIndex < 0 || startIndex < 0) {
        throw new IllegalArgumentException("Index cannot be negative!");
    }

    N = stopIndex - startIndex;
    filteredData = new float[N];
    complexData = new Complex[N];
    filteredComplexData = new Complex[N];

    for (int i = startIndex; i < stopIndex; i++) {
        complexData[i-startIndex] = new Complex(data[i], 0.0f);
        filteredComplexData[i-startIndex] = new Complex(0.0f, 0.0f);
    }

    for (int k = 0; k < N; k++) {
        for (int n = 0; n < N; n++) {
            Complex c = complexData[n].copy();
            filteredComplexData[k].add(c.mul(new Complex(1.0f, (float)(-2*Math.PI*n*k/(float)N))));
        }
    }

    for (int i = 0; i < N; i++) {
        filteredData[i] = filteredComplexData[i].abs();
    }

    return filteredData;
}
}

现在，两者似乎都给出了 [8.0, 4.0, 8.0, 0.0] 的正确答案，即 [20.0, 4.0j, 12.0, -4.0j]。但是，如果我给他们提供由以下各项产生的正弦波:

mBuffer = new float[1024];
float sampleRate = 1000.0f;
float frequency = 50.0f;

for (int i = 0; i < mBuffer.length; i++) {
    mBuffer[i] = (float)(0.5*Math.sin(2*Math.PI*i*frequency/sampleRate));
}

Complex的实现供引用:

public final class Complex {
public float mR, mTheta;

public Complex() {
    mR = 0.0f;
    mTheta = 0.0f;
}

public Complex(float r, float theta) {
    mR = r;
    mTheta = theta;
}

public Complex copy() {
    return new Complex(mR, mTheta);
}

public Complex add(Complex c) {
    float real, imag;
    real = (float)(mR * Math.cos(mTheta) + c.mR * Math.cos(c.mTheta));
    imag = (float)(mR * Math.sin(mTheta) + c.mR * Math.sin(c.mTheta));

    mR = (float)Math.sqrt(Math.pow(real, 2) + Math.pow(imag, 2));

    if (real != 0.0f) {
        mTheta = (float)Math.atan(imag / real);
    } else {
        mTheta = (float)(imag > 0.0f ? Math.PI/2.0f : Math.PI*3.0f/2.0f);
    }
    return this;
}

public Complex sub(Complex c) {
    float real, imag;
    real = (float)(mR * Math.cos(mTheta) - c.mR * Math.cos(c.mTheta));
    imag = (float)(mR * Math.sin(mTheta) - c.mR * Math.sin(c.mTheta));

    mR = (float)Math.sqrt(Math.pow(real, 2) + Math.pow(imag, 2));

    if (real != 0.0f) {
        mTheta = (float)Math.atan(imag / real);
    } else {
        mTheta = (float)(imag > 0.0f ? Math.PI/2.0f : Math.PI*3.0f/2.0f);
    }
    return this;
}

public Complex mul(Complex c) {
    mR = mR * c.mR;
    mTheta = mTheta + c.mTheta;
    return this;
}

public Complex div(Complex c) {
    mR = mR / c.mR;
    mTheta = mTheta - c.mTheta;
    return this;
}

public Complex pow(float exp) {
    mTheta = mTheta * exp;
    mR = (float)Math.pow(mR, exp);
    return this;
}

public float abs() {
    return mR;
}

public float getRealPart() {
    return (float)(mR * Math.cos(mTheta));
}

public float getImagPart() {
    return (float)(mR * Math.sin(mTheta));
}

public String toStringRectangular() {
    float real, imag;
    StringBuilder sb = new StringBuilder();

    real = (float)(mR * Math.cos(mTheta));
    imag = (float)(mR * Math.sin(mTheta));

    sb.append(real);
    if (imag >= 0) {
        sb.append(" + ");
    } else {
        sb.append(" - ");
    }
    sb.append(Math.abs(imag));
    sb.append("i");

    return sb.toString();

}

public String toStringExponential() {
    StringBuilder sb = new StringBuilder();

    sb.append(mR);
    sb.append(" * e ^ ");
    sb.append(mTheta);
    sb.append("i");

    return sb.toString();

}

public String toString() {
    return toStringExponential() + " [ " + toStringRectangular() + " ] ";
}

public static Complex[] getInitializedArray(int size) {
    Complex[] arr = new Complex[size];

    for (int i = 0; i < arr.length; i++) {
        arr[i] = new Complex(0.0f, 0.0f);
    }

    return arr;
}
}

Answer 1

您的 FFT 实现似乎是合理的。但是，使用 Math.atan 存在问题(返回 [-pi/2,pi/2] 内的值，而不是整个 [-pi,pi] 范围)在 Complex 的 add 和 sub。

要解决此问题，您应该使用:

mTheta = (float)Math.atan2(imag, real);