Amateur DFT (Discrete Fourier Transform) implementation in C++

dft.cc

/**
 *   (c) Copyright 2020 Tripulse.
 * 
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */

/// This DFT implementation is intended for educational purposes,
/// if you are planning to use it in production, DO NOT.
///
/// DFT is, by nature, O(n^2) so as this program. FFT, which is
/// algorithmically superior to DFT because it is O(n * log n).

#include <complex>
#include <vector>

using namespace std::literals;

/**
 * \brief Computes 1-dimensional DFT (discrete fourier transform) on given input.
 * 
 * \param input  an array of complex numbers, possibly a signal in time-domain.
 * \param output a modifiable array to store the output in, if larger than input
 *               only elements past input length are changed and if smaller the
 *               result is truncated to output length.
 */
void dft(const std::vector<std::complex<double>>& input,
               std::vector<std::complex<double>>& output)
{
    auto M_I_2PI_DL = -(6.28318530718i / (double)input.size());

    for (size_t k = 0; k < output.size(); ++k) {
        output[k] = 0;
        for (size_t n = 0; n < input.size(); ++n)
            output[k] += input[n] * pow(2.718281828459045, M_I_2PI_DL * (double)k * (double)n);
    }
}

/**
 * \brief Computes 1-dimensional IDFT (inverse discrete fourier transform) on given input.
 * 
 * \param input an array of complex numbers, possibly a DFT'd signal in frequency domain.
 * \param output a modifiable array to store the output in, if larger than input
 *               only elements past input length are changed and if smaller the
 *               result is truncated to output length.
 */
void idft(const std::vector<std::complex<double>>& input,
                std::vector<std::complex<double>>& output)
{
    for (size_t k = 0; k < output.size(); ++k) {
        output[k] = 0;
        for (size_t n = 0; n < input.size(); ++n)
            output[k] += input[n] * pow(2.718281828459045, (6.28318530718i * (double)k * (double)n) / (double)input.size());

        output[k] /= input.size();
    }
}

I don't even know if this even works, didn't compare the results against a working implementation.
It works quite fine as expected, I build up a prototype in Python and tested against NumPy.

Input: [0, 1, 2]

numpy.fft.fft(): [3+0j, -1.5+0.8660254j, -1.5-0.8660254j]
          dft(): [3+0j, -1.4999999999996418+0.8660254037847832j, -1.5000000000007163-0.8660254037837494j]

Don't use this in production, just use FFTW which is a lot better than this.

'-'

''
-

When DFT'ed input was IDFT'ed it produces (possibly garbage) values in imaginary part whilst input being purely real. Though the real part is completely accurate.

Is this a floating-point error or problem in implementation? This remains a mystery for now.

'-'

'-'