1

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;
}
Guy Yafe
  • 111
  • 1
  • 2
    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. – Dan Boschen Jul 07 '20 at 10:46
  • 2
    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 – Hilmar Jul 07 '20 at 10:56

0 Answers0