I am trying to implement a Cooley-Tukey FFT algorithm in its iterative version. Specifically, I am implementing this Wikipedia pseudocode (please look at the section Data reordering, bit reversal, and in-place algorithms which I couldn't link properly).
The test I am carrying out:
- Simulate a 10 Hz sinusoidal wave, with 100Hz sampling rate.
- Putting this data in a buffer
- Applying the algorithm to the buffer
- Comparing to a 3rd party algorithm to verify that the data was created successfully.
Since the FFT is working with an imaginary part, I should expect a mirroring of the frequencies, and windeed When I ran a 3rd party algorithm (Apache Math) on the synthetic data, I got correct, accurate results.
When I ran my own implementation on the buffer I got correct frequency peaks at 10 HZ and at 90Hz, but these peaks weren't symmetric. Each had a different height.
My question is if I should expect such results, or there is a mistake in my implementation or in my understanding of the algorithm.
Specifically I wonder if the terms ω ωm, ω A[k + j + m/2], etc in the reference pseudocode indeed represent a simple two complex numbers multiplication?
The results from the 3rd party algorithm:
The results from my implementation:
The pseudo code from Wikipedia:
algorithm iterative-fft is
input: Array a of n complex values where n is a power of 2.
output: Array A the DFT of a.
bit-reverse-copy(a, A)
n ← a.length
for s = 1 to log(n) do
m ← 2s
ωm ← exp(−2πi/m)
for k = 0 to n-1 by m do
ω ← 1
for j = 0 to m/2 – 1 do
t ← ω A[k + j + m/2]
u ← A[k + j]
A[k + j] ← u + t
A[k + j + m/2] ← u – t
ω ← ω ωm
return A
My own implementation:
public void calcFft(double[] timeDomainData, double[][] frequencyDomainData) {
int bitsLength = calcLog2(timeDomainData.length);
bitReverseCopy(timeDomainData, frequencyDomainData, bitsLength);
int sampleSize = timeDomainData.length;
for (int iteration = 1; iteration <= bitsLength; iteration++) {
int m = 1 << iteration;
double twiddleReal = Math.cos(2 * Math.PI / ((double) m));
double twiddleIm = -Math.sin(2 * Math.PI / ((double) m));
for(int k = 0; k < sampleSize; k += m){
double omegaReal = 1.0;
double omegaIm = 0.0;
for (int j = 0; j < m / 2; j++){
double currentFrequencyForTReal = frequencyDomainData[REAL][k + j + m / 2];
double currentFrequencyForTIm = frequencyDomainData[IM][k + j + m / 2];
double currentFrequencyForUReal = frequencyDomainData[REAL][k + j];
double currentFrequencyForUIm = frequencyDomainData[IM][k + j];
double tReal = omegaReal * currentFrequencyForTReal - omegaIm * currentFrequencyForTIm;
double tIm = omegaIm * currentFrequencyForTReal + omegaReal * currentFrequencyForTIm;
double uReal = currentFrequencyForUReal;
double uIm = currentFrequencyForUIm;
frequencyDomainData[REAL][k + j] = uReal + tReal;
frequencyDomainData[IM][k + j] = uIm + tIm;
frequencyDomainData[REAL][k + j + m / 2] = uReal - tReal;
frequencyDomainData[IM][k + j + m / 2] = uIm - tIm;
omegaReal = omegaReal * twiddleReal - omegaIm * twiddleIm;
omegaIm = omegaReal * twiddleIm + omegaIm * twiddleReal;
}
}
}
}
private void bitReverseCopy(double[] timeDomainData, double[][] frequencyDomainData, int bitsLength){
for (int timeDomainIndex = 0; timeDomainIndex < timeDomainData.length; timeDomainIndex++) {
int frequencyDomainIndex = reverse(timeDomainIndex, bitsLength);
frequencyDomainData[REAL][frequencyDomainIndex] = timeDomainData[timeDomainIndex];
frequencyDomainData[IM][frequencyDomainIndex] = 0;
}
}
public int reverse(int value, int bitsLength) {
int reversed = 0;
for (int step = 0; step < bitsLength; step++) {
int bit = value & 1;
reversed <<= 1;
reversed |= bit;
value >>>= 1;
}
return reversed;
}
public int calcLog2(int length) {
switch (length){
case 0:
return 0;
case 1 << 1:
return 1;
case 1 << 2:
return 2;
case 1 << 3:
return 3;
case 1 << 4:
return 4;
case 1 << 5:
return 5;
case 1 << 6:
return 6;
case 1 << 7:
return 7;
case 1 << 8:
return 8;
case 1 << 9:
return 9;
case 1 << 10:
return 10;
case 1 << 11:
return 11;
case 1 << 12:
return 12;
case 1 << 13:
return 13;
case 1 << 14:
return 14;
case 1 << 15:
return 15;
case 1 << 16:
return 16;
case 1 << 17:
return 17;
case 1 << 18:
return 18;
case 1 << 19:
return 19;
case 1 << 20:
return 20;
case 1 << 21:
return 21;
case 1 << 22:
return 22;
case 1 << 23:
return 23;
case 1 << 24:
return 24;
case 1 << 25:
return 25;
case 1 << 26:
return 26;
case 1 << 27:
return 27;
case 1 << 28:
return 28;
case 1 << 29:
return 29;
case 1 << 30:
return 30;
case 1 << 31:
return 31;
default:
throw new IllegalArgumentException("Array size must be a power of 2");
}
}
My Sinusoidal signal generation:
private double[] createSamplesWithExtrapolation(double frequencyHz, double sampleRateHz, double sampleTimeSeconds) {
int numberOfSamples = ((int) (sampleTimeSeconds * sampleRateHz)) + 1;
numberOfSamples = shiftToPowerOf2(numberOfSamples);
double[] samples = new double[numberOfSamples];
for (int timeStep = 0; timeStep < numberOfSamples; timeStep++) {
double t = 1.0 / sampleRateHz * timeStep;
double sample = Math.sin(2 * Math.PI * frequencyHz * t);
samples[timeStep] = sample;
}
return samples;
}

