fft test

gistfile1.cpp

#define _USE_MATH_DEFINES
#include <memory>
#include <complex>
#include <iostream>
#include <cmath>
#include "bmp.hpp"

using complex_t = std::complex<float>;

int bit_rev(int a){
    if(a == 0){
        return a;
    }
    int n = 0, t = a;
    while(t > 0){
        ++n;
        t >>= 1;
    }
    int r = 0;
    if(n % 2 == 1){
        r |= ((a >> (n / 2)) & 1) << (n / 2);
    }
    for(int i = 0; i < n / 2; ++i){
        r |= ((a >> i) & 1) << (n - 1 - i);
        r |= ((a >> (n - 1 - i)) & 1) << i;
    }
    return r;
}

void bit_rev_copy(const float *a, complex_t *A, int n){
    for(int i = 0; i < n; ++i){
        A[bit_rev(i)] = a[i];
    }
}

// a[1 << lg_n] : input
// A[1 << lg_n] : output
void fft(const float *a, complex_t *A, int lg_n){
    int n = 1 << lg_n;
    bit_rev_copy(a, A, n);
    for(int s = 1; s <= lg_n; ++s){
        int m = 1 << s;
        complex_t omega_m(std::cos(2.0 * M_PI / m), std::sin(2.0 * M_PI / m)), omega = 1.0;
        for(int j = 0; j < m / 2; ++j){
            for(int k = j; k < n; k += m){
                complex_t t = omega * A[k + m / 2], u = A[k];
                A[k] = u + t;
                A[k + m / 2] = u - t;
            }
            omega *= omega_m;
        }
    }
}

int main(){
    std::unique_ptr<float[]> input(new float[128]{ 0 });
    std::unique_ptr<complex_t[]> output(new complex_t[128]{ complex_t(0, 0) });
    for(int i = 0; i < 128; ++i){
        // 32Hz
        input.get()[i] += std::sin(32 * 2.0 * M_PI * static_cast<float>(i) / 128);
        input.get()[i] += std::sin(16 * 2.0 * M_PI * static_cast<float>(i) / 128);
    }
    fft(input.get(), output.get(), 7);
    float max = 0.0;
    for(int i = 0; i < 128; ++i){
        if(std::abs(output.get()[i].real()) > max){
            max = std::abs(output.get()[i].real());
        }
    }
    tt_legacy::bmp bmp(24, 64, 100);
    for(int i = 0; i < 100; ++i){
        for(int j = 0; j < 64; ++j){
            bmp.clr(j, i, bmp.rgb(0xFF, 0xFF, 0xFF));
        }
    }
    for(int i = 0; i < 64 - 1; ++i){
        bmp.line(i, 99 - static_cast<int>(std::abs(output.get()[i].real()) / max * 100.0), i + 1, 99 - static_cast<int>(std::abs(output.get()[i + 1].real()) / max * 100.0), bmp.rgb(0xFF, 0x00, 0x00));
    }
    bmp.write("put.bmp");
    return 0;
}