DFT demo

pi.cpp

/* g++ pi.cpp -o pi -std=c++14 -lSDL2 */
/* Grab pi.csv from https://dsp.stackexchange.com/a/59181/70343 */

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <complex>
#include <vector>
#include <algorithm>

#define SDL_MAIN_HANDLED
#include <SDL2/SDL.h>

const double pi = std::acos(-1.0);

std::vector<std::complex<double>> DFT(const std::vector<std::complex<double>>& x)
{
    int N = x.size();

    std::vector<std::complex<double>> X;

    for (int k = 0; k < N; ++k)
    {
        std::complex<double> X_k(0, 0);

        for (int t = 0; t < N; ++t)
        {
            X_k += x[t] * std::exp(std::complex<double>(0, -2 * pi * k * t / N));
        }

        X.push_back(X_k);
    }

    return X;
}

std::vector<std::complex<double>> IDFT(const std::vector<std::complex<double>>& X)
{
    int N = X.size();

    std::vector<std::complex<double>> x;

    for (int t = 0; t < N; ++t)
    {
        std::complex<double> x_t(0, 0);

        for (int k = 0; k < N; ++k)
        {
            x_t += X[k] * std::exp(std::complex<double>(0, 2 * pi * k * t / N));
        }

        x_t /= std::complex<double>(N, 0);

        x.push_back(x_t);
    }

    return x;
}

/* https://stackoverflow.com/a/48291620 */
void DrawCircle(SDL_Renderer* renderer, int32_t centreX, int32_t centreY, int32_t radius)
{
    const int32_t diameter = (radius * 2);

    int32_t x = (radius - 1);
    int32_t y = 0;
    int32_t tx = 1;
    int32_t ty = 1;
    int32_t error = (tx - diameter);

    while (x >= y)
    {
        //  Each of the following renders an octant of the circle
        SDL_RenderDrawPoint(renderer, centreX + x, centreY - y);
        SDL_RenderDrawPoint(renderer, centreX + x, centreY + y);
        SDL_RenderDrawPoint(renderer, centreX - x, centreY - y);
        SDL_RenderDrawPoint(renderer, centreX - x, centreY + y);
        SDL_RenderDrawPoint(renderer, centreX + y, centreY - x);
        SDL_RenderDrawPoint(renderer, centreX + y, centreY + x);
        SDL_RenderDrawPoint(renderer, centreX - y, centreY - x);
        SDL_RenderDrawPoint(renderer, centreX - y, centreY + x);

        if (error <= 0)
        {
            ++y;
            error += ty;
            ty += 2;
        }

        if (error > 0)
        {
            --x;
            tx += 2;
            error += (tx - diameter);
        }
    }
}

int main()
{
    std::vector<std::complex<double>> data;

    std::ifstream data_file("pi.csv");

    if (!data_file.is_open())
    {
        std::cerr << "[IO error] Unable to open pi.csv\n";
        std::exit(1);
    }

    std::string line;

    while (std::getline(data_file, line))
    {
        int x = std::stoi(line.substr(0, line.find(",")));
        int y = std::stoi(line.substr(line.find(",") + 1));

        data.push_back(std::complex<double>(x + 300, -y + 300));
    }

    data_file.close();

    int N = data.size();

    std::vector<std::complex<double>> data_dft = DFT(data);

    /* uncomment the code below for the basic demonstration of filtering */
/*
    std::complex<double> peak_component = *std::max_element(data_dft.begin(), data_dft.end(), [](const auto& lhs, const auto& rhs) {return std::abs(lhs) < std::abs(rhs);});

    for (int k = 0; k < N; ++k)
    {
        if (std::abs(data_dft[k]) < std::abs(peak_component) / 30)
        {
            data_dft[k] = 0;
        }
    }

    data = IDFT(data_dft);
*/
    if (SDL_Init(SDL_INIT_EVERYTHING) < 0)
    {
        SDL_LogError(SDL_LOG_CATEGORY_APPLICATION, "Couldn't initialize SDL: %s", SDL_GetError());
        std::exit(1);
    }

    SDL_Window *window = SDL_CreateWindow(
        "PI",
        SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED,
        600, 600,
        SDL_WINDOW_SHOWN
    );

    SDL_Renderer *renderer = SDL_CreateRenderer(
        window,
        -1,
        SDL_RENDERER_ACCELERATED | SDL_RENDERER_TARGETTEXTURE
    );

    SDL_Texture *texture = SDL_CreateTexture(
        renderer,
        SDL_PIXELFORMAT_RGBA8888,
        SDL_TEXTUREACCESS_TARGET,
        600, 600
    );

    bool quit = false;
    SDL_Event event;

    int t = 1;

    SDL_SetRenderTarget(renderer, texture);
    SDL_SetRenderDrawColor(renderer, 255, 255, 255, 255);
    SDL_RenderClear(renderer);
    SDL_SetRenderTarget(renderer, NULL);

    while (!quit)
    {
        SDL_Delay(30);
        SDL_PollEvent(&event);

        switch (event.type)
        {
            case SDL_QUIT:
                quit = true;
                break;
        }

        SDL_SetRenderDrawColor(renderer, 255, 255, 255, 255);
        SDL_RenderClear(renderer);

        SDL_SetRenderTarget(renderer, texture);
        SDL_SetRenderDrawColor(renderer, 0, 0, 255, 255);
        SDL_RenderDrawLine(renderer, data[(t - 1) % N].real(), data[(t - 1) % N].imag(), data[t % N].real(), data[t % N].imag());
        SDL_SetRenderTarget(renderer, NULL);

        SDL_RenderCopy(renderer, texture, NULL, NULL);

        SDL_SetRenderDrawColor(renderer, 255, 0, 0, 255);

        std::complex<double> prev(0, 0);
        std::complex<double> x_t(0, 0);

        for (int k = 0; k < N; ++k)
        {
            std::complex<double> k_contrib = data_dft[k] * std::exp(std::complex<double>(0, 2 * pi * k * t / N)) / std::complex<double>(N, 0);

            x_t += k_contrib;

            SDL_RenderDrawLine(renderer, prev.real(), prev.imag(), x_t.real(), x_t.imag());
            DrawCircle(renderer, prev.real(), prev.imag(), std::abs(k_contrib));

            prev = x_t;
        }

        SDL_RenderPresent(renderer);

        t++;
    }

    SDL_DestroyTexture(texture);
    SDL_DestroyRenderer(renderer);
    SDL_DestroyWindow(window);
    SDL_Quit();

    return 0;
}