0
$\begingroup$

I am implementing radix-2 DIT FFT (real only) in C++. For input sizes > 8 I can't seem to match the output of FFTW. Here's my implementation of the FFT:

std::vector<std::complex<double>> MyFFT::transform (const std::vector<double>& samples) {
    compute_twiddles (samples.size ());
    std::vector<double> samples_copy {samples};
    bit_reversal (samples_copy);

    // Copy real samples to internal complex vector.
    m_samples.reserve (samples.size ());
    for (double sample : samples_copy)
        m_samples.emplace_back (sample, 0);

    fft (m_samples, 0, samples.size () - 1);
    return m_samples;
}


// Places all the twiddle factors in a flat 1D array
// Twiddle factors for 2 point FFT reside in indices 0 and 1
// Twiddle factors for 4 point FFT reside in indices 2, 3, 4 and 5 and so on.
void MyFFT::compute_twiddles (size_t fft_size) {
    // Geometric series summation formula
    size_t size = -(2 * ((1 - pow (2, log2 (fft_size)))));
    m_twiddles.reserve (size);
    m_twiddles.emplace_back (1, 0);
    m_twiddles.emplace_back (-1, 0);
    for (size_t m = 4; m <= fft_size; m *= 2) {
        for (size_t i = 0; i < m; i++) {
            std::complex<double> exponent = {0, (-2 * PI * i) / m};
            auto tmp = exp (exponent);
            m_twiddles.emplace_back (tmp);
        }
    }
}

void MyFFT::bit_reversal (std::vector<double>& samples) {
    // TODO: check if sample vector size is a power of 2.
    size_t size = samples.size ();
    int bit_space = (int) log2 (size);
    for (int i = 0; i < size / 2; i++) {
        int complement_i = get_complement (i, bit_space);
        auto tmp = samples[i];
        samples[i] = samples[complement_i];
        samples[complement_i] = tmp;
    }
}

// Naive beyond measures
size_t MyFFT::get_complement (size_t number, size_t width) {
    // Unpack bits
    std::vector<u8> bits;
    bits.reserve (width);
    for (size_t i = 0; i < width; i++)
        bits.push_back ((number & (1u << i)) > 0);

    // Reverse bits
    for (size_t i = 0, j = width - 1; i < j; i++, j--) {
        u8 tmp = bits.at (i);
        bits[i] = bits[j];
        bits[j] = tmp;
    }

    // Compute "complement"
    size_t accum = 0;
    for (int i = 0; i < width; i++)
        accum += bits[i] * (pow (2, i));
    return accum;
}

void MyFFT::fft (std::vector<std::complex<double>>& samples, size_t start, size_t end) {
    size_t mid = (start + end) / 2;
    size_t fft_size = (end - start) + 1;
    if (fft_size == 2) {
        synth (samples, start, end);
        return;
    }
    fft (samples, start, mid);
    fft (samples, mid + 1, end);
    synth (samples, start, end);
}

void MyFFT::synth (std::vector<std::complex<double>>& samples, size_t start, size_t end) {
    size_t fft_size = (end - start) + 1;
    for (size_t k = 0, m = start; k < fft_size / 2; k++, m++) {
        auto out1 = samples[m] + (twiddle (fft_size, k) * samples[m + (fft_size / 2)]);
        auto out2 = samples[m] + (twiddle (fft_size, k + (fft_size / 2)) * samples[m + (fft_size / 2)]);
        samples[m] = out1;
        samples[m + (fft_size / 2)] = out2;
    }
}

std::complex<double> MyFFT::twiddle (size_t fft_size, size_t k) {
    double numerator = 1 - pow (2, log2 (fft_size / 2));
    size_t segment = -(2 * numerator);
    return m_twiddles[segment + k];
}

Input: [1, 1.26007, -0.221232, -1.3968, -0.64204, 1, 1.26007, -0.221232, -1.3968, -0.64204, 1, 1.26007, -0.221232, -1.3968, -0.64204, 1]

Outputs (Format is [i: index, myfft_output: {real, imag}, fftw_output: {real, imag}]:

[i: 0, out1: {1, 0},                    out2: {1, 0}]
[i: 1, out1: {2.87583, -0.204185},      out2: {1.1399, -0.0809334}]
[i: 2, out1: {1.41061, 3.94722},        out2: {1.83142, -0.230978}]
[i: 3, out1: {5.85787, -0.912478},      out2: {10.4696, -1.63084}]
[i: 4, out1: {-1.81526, -4.89294},      out2: {-2.65688, 0.420808}]
[i: 5, out1: {3.94792, -0.544678},      out2: {-1.2589, 0.173685}]
[i: 6, out1: {-1.31929, 3.42707},       out2: {-0.898483, 0.0904908}]
[i: 7, out1: {-3.09441, 0.163614},      out2: {-0.763371, 0.0403626}]
[i: 8, out1: {-0.726543, -3.09897e-09}, out2: {-0.726543, 0}]
[i: 9, out1: {-3.09441, -0.163614},     out2: {-0.763371, -0.0403626}]
[i: 10, out1: {-1.31929, -3.42707},     out2: {-0.898483, -0.0904908}]
[i: 11, out1: {3.94792, 0.544678},      out2: {-1.2589, -0.173685}]
[i: 12, out1: {-1.81526, 4.89294},      out2: {-2.65688, -0.420808}]
[i: 13, out1: {5.85787, 0.912478},      out2: {10.4696, 1.63084}]
[i: 14, out1: {1.41061, -3.94722},      out2: {1.83142, 0.230978}]
[i: 15, out1: {2.87583, 0.204185},      out2: {1.1399, 0.0809334}]

2, 4 and 8 point FFTs don't cause trouble but going beyond that, as in this case a 16-point FFT outputs results that are completely different than those of FFTW.

$\endgroup$
6
  • $\begingroup$ I am implementing FFT purely for learning purposes. I don't intend to use my implementation in production. $\endgroup$ – Devashish Jaiswal Mar 6 at 18:45
  • $\begingroup$ Very good but, in its presented form, this sounds more like a programming problem than a signal processing question... It seems you made an error at the implementation stage... $\endgroup$ – Fat32 Mar 6 at 19:03
  • $\begingroup$ Standard debug procedure: 1) verify each individual step by itself:(twiddle factors, ordering, indexing. 2) Choose a simple signal with a known answer: dirac delta, delayed dirac delta, single-bin sine wave. Single step through the ones that don't check out. $\endgroup$ – Hilmar Mar 6 at 19:21
  • $\begingroup$ Thanks. Since I'll be performing a lot of tests trying out different signals, I was wondering if there is any library that eases the process of generating "mixed" signal samples. I know it's trivial to generate such samples but I would rather not waste time writing code for that. $\endgroup$ – Devashish Jaiswal Mar 6 at 19:28
  • $\begingroup$ @Fat32 Alright! I think I fixed the issue. Turns out the bit reversal algorithm was not only naive, but also wrong. Anyway, how much deviation from fftw's output should be tolerated? Computing 2 to 2^17 point FFTs and comparing real and imaginary parts to that of FFTW's I see deviations less than 0.005. But computing 2^18 point FFT, there are 2 bins where the (real component) errors are slightly greater than 0.005. Should I just ignore that? I want to believe that those differences arise because of floating point arithmetic and there's nothing wrong with my algorithm implementation. $\endgroup$ – Devashish Jaiswal Mar 7 at 4:16
2
$\begingroup$

Reading the comments it seems that your question is not simplified to "how to measure accuracy of the FFT"

I would say that instead of comparing with the accuracy of a different implementation you could compare with itself. After all it is a good exercise to be able to test your own implementations, you will need this skill if one day you attempt to implement any original numeric algorithm.

With the FFT is implementing the following sum

$$ FFT(x)_l = \sum_{k=0}^{N-1} x_k exp\left(-\frac{2\pi (-1)^{1/2} l k}{N}\right) $$

The inverse FFT would be

$$ IFFT(x)_k = \frac{1}{N} \sum_{k=0}^{N-1} X_l exp\left(\frac{2\pi (-1)^{1/2} l k}{N}\right) $$

Comparing two expressions we have what for me is the first handy property of FFT

$$ \begin{eqnarray}IFFT(x)_k &=&\frac{1}{N} conj\left(\sum_{k=0}^{N-1} conj\left(X_l exp\left(\frac{2\pi (-1)^{1/2} l k}{N}\right) \right) \right) \\ &=&\frac{1}{N} conj\left(\sum_{k=0}^{N-1} conj\left(X_l\right) exp\left(-\frac{2\pi (-1)^{1/2} l k}{N}\right) \right) \end{eqnarray} $$

Thus ifft(X) = conj(fft(conj(X)))/N,

vector<complex<double>> ifft(vector<complex<double>> &samples, size_t start, size_t end){
  const size_t N = (end - start);
  for(int i = 0; i < N; ++i){
    samples[i] = conj(samples[i]);
  }
  fft(samples, start, end);
  for(int i = 0; i < N; ++i){
    samples[i] = conj(samples[i])/N;
  }
}

Using this you can generate random data and compute x - ifft(fft(x)) and test your implementation independent of any other implementation.

Theoretically the error must be less than $2^{1-s}N\max(|x|)$, where s is the number of significative bits (53 for double, 24 for float, 11 for half). If the 0.005 you mentioned is greater than one billionth of your maximum input I would say that you need to review the algorithm.

Typically the error will be close to the errors in a multiplication or addition, and will grow with the (square root of) logarithm of the length of the input, in a FFT implementation. With double precision the relative error should be around 1e-15.

The following plot shows the maximum error for a white normal noise input. max error for a random normal input

$\endgroup$
10
  • 1
    $\begingroup$ Hey thanks! I had plans to implement ifft next but I wanted to test the forward fft first. I assume the graph is showing relative error between x(n) and ifft(fft(x(n)))? If errors shouldn't exceed 1e-15 even across implementations (for example, between my fft and fftw) then 0.005 feels too much. $\endgroup$ – Devashish Jaiswal Mar 7 at 13:37
  • $\begingroup$ Slightly off-topic, but I would love to get some constructive criticism on my implementation. Here's the code if you don't mind going through it once: github.com/devsh0/fft-playground $\endgroup$ – Devashish Jaiswal Mar 7 at 13:40
  • 1
    $\begingroup$ @DevashishJaiswal: if you want feedback on code quality, I suggest asking a separate question like "How to write a good FFT" and specify what criteria you associate with "good". Sorry, at first blush your code is very inefficient and probably not great numerically either. $\endgroup$ – Hilmar Mar 7 at 17:03
  • 1
    $\begingroup$ Yes, I think 0.005 is too much, but I cannot say because it depends on the input you are using. In terms of making the code fast there is a lot of modern CPU instructions, and challenges to make better use of cache. But I think that if you are doing this as an exercise, your challenge should be "how to write the most readable possible implementation" $\endgroup$ – Bob Mar 8 at 9:21
  • $\begingroup$ @user12750353. Readability is in the eyes of the beholder, I guess. Personally I would think that something that follows the standard block diagram of an FFT is both either to read and way more efficient than this recursive approach. $\endgroup$ – Hilmar Mar 8 at 13:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.