<!--親の顔より見た光景-->

日々の発見を残していきます。

【C++】離散フーリエ級数展開

理解するためにすべて手書きで書いた。
乱数で作ったグラフを、
三角関数の足し合わせのみでキレイに再現できている。
f:id:akagi13213:20190203015531p:plain

(あと一息で離散フーリエ変換ができそう)

#include "pch.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>

#define M_PI    3.14159265358979323846

int main(void)
{
    std::ifstream ifs( "./source.csv" );
    if ( ifs.fail() ) {
        std::cout << "file read error..." << std::endl;
        return -1;
    }

    std::string str;
    std::vector<double> f;
    while ( getline( ifs, str ) ) {
        f.push_back(std::stod( str ));
    }

    double a0 = 0;
    std::vector<double> an;
    std::vector<double> bn;

    double sum = 0;
    int N = f.size();
    for ( int x = 0; x < N; x++ ) {
        sum += f[x];
    }
    a0 = sum / N;

    for ( int n = 1; n < N; n++ ) {
        double a_sum = 0;
        double b_sum = 0;
        
        for ( int x = 0; x < N; x++ ) {
            a_sum += f[x] * cos( 2 * M_PI * n * x  / N );
            b_sum += f[x] * sin( 2 * M_PI * n * x / N );
        }
        
        if ( n == ( N / 2 - 1 ) ) {
            an.push_back( a_sum / N );
            bn.push_back( b_sum / N );
        }
        else {
            an.push_back( 2 * a_sum / N );
            bn.push_back( 2 * b_sum / N );
        }
    }
    
    std::vector<double> fx_fourier;
    for (int x = 0; x < N; x++) {
        double a_sum = 0;
        double b_sum = 0;
        for ( int n = 1; n < N; n++ ) {
            a_sum += an[n - 1] * cos( 2 * M_PI * n * x / N );
            b_sum += bn[n - 1] * sin( 2 * M_PI * n * x / N );
        }
        fx_fourier.push_back( a0 + a_sum + b_sum );
    }

    std::ofstream dist;
    dist.open("dist.csv", std::ios::out);
    for ( int i = 0; i < fx_fourier.size(); i++ ) {
        dist << fx_fourier[i] << std::endl;
    }
    
    return 0;
}