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();
    }
}

'-'

'-'