1
$\begingroup$

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:

  1. Simulate a 10 Hz sinusoidal wave, with 100Hz sampling rate.
  2. Putting this data in a buffer
  3. Applying the algorithm to the buffer
  4. 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 the 3rd party algorithm:

The results from my implementation:

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;
}
$\endgroup$
2
  • 2
    $\begingroup$ I suspect you are off by a sample count somewhere (did you used the exact same waveform for each, and they both returned the same number of samples?) I recommend comparing your algorithm and the 3rd party fft with short sequences such as [1 2 3 4 5 8 7 6] and [j -2 2j -1 5 2 3 5]. You should then be able to step through each of the stages to identify the error. $\endgroup$ Jul 7, 2020 at 10:46
  • 2
    $\begingroup$ Adding to Dan's point: it looks like you are close, so debug with something as simple as possilble: start with all zeros and verify that the output is all zero. Then try [1 0 0 0 ] and make sure the output is all ones.. Then do try [0 1 0 0] and make sure you get a complex exponential. . Work your way up to larger lengths and more complex signal until you find the issue $\endgroup$
    – Hilmar
    Jul 7, 2020 at 10:56

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.