1
Я реализовал рекурсивный фрейм-код DDR FFT на Java и регулярный DFT, чтобы проверить мои результаты с помощью FFT, но результаты этих двух отличаются, и я не могу понять это , Оба передают весь массив с помощью методов apply(), индекс start и stop равен 0 и data.length соответственно. Версия DFT выглядит корректно с хорошим пиком в бункере 50, а FFT - мусором. Что я делаю не так?Устранение неисправностей DIT FFT Radix-2 Algorithm
Это реализация FFT (адаптировано из http://www.engineeringproductivitytools.com/stuff/T0001/PT04.HTM, я проверил, сравнивая его с псевдокодом на https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#Pseudocode «рутинный Рекурсивном DIT FFT.»):
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));
}
}
И это реализация ДПФА:
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));
}
Реализация комплекса для справки:
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;
}
}
Спасибо, что решить все! :) – phnom