コンパイル時離散フーリエ変換(DFT)

2012-04-30

Sprout.Numeric.DFT は、constexpr DFT アルゴリズムを提供します。

https://github.com/bolero-MURAKAMI/Sprout/tree/master/sprout/numeric/dft

DFT は周波数解析などに用いられますが、ここでは詳しくは述べません。

とりあえずサンプルコードを見てください。

単純な正弦波の DFT

#include <iostream>
#include <fstream>
#include <sprout/array.hpp>
#include <sprout/pit.hpp>
#include <sprout/complex.hpp>
#include <sprout/numeric/dft.hpp>
#include <sprout/range/adaptor.hpp>
#include <sprout/functional.hpp>

template<typename Elem, typename Traits, typename InputRange>
std::basic_ostream<Elem, Traits>& output_plot(std::basic_ostream<Elem, Traits>& os, InputRange const& range) {
    os << std::fixed;
    int x = 0;
    for (auto const& e : range) {
        os << x++ << ',' << e << '\n';
    }
    return os;
}
template<typename Elem, typename Traits, typename InputRange>
std::basic_ostream<Elem, Traits>& output_plot_real(std::basic_ostream<Elem, Traits>& os, InputRange const& range) {
    os << std::fixed;
    int x = 0;
    for (auto const& e : range) {
        os << x++ << ',' << real(e) << '\n';
    }
    return os;
}

int main() {
    using namespace sprout;

    constexpr std::size_t size = 256;
    typedef array<complex<double>, size> container_t;
    typedef array<double, size> spectrum_t;

    // 周波数 10 の正弦波を生成
    constexpr container_t src_data = adaptors::sinusoidal(10. / size, 10000.) | adaptors::copied;
    {
        std::ofstream os("src.txt");
        output_plot_real(os, src_data);
    }

    // DFT(離散フーリエ変換)
    constexpr auto dft_data = sprout::dft(begin(src_data), end(src_data), pit<container_t>());

    // IDFT(逆離散フーリエ変換)
    constexpr auto idft_data = sprout::idft(begin(dft_data), end(dft_data), pit<container_t>());
    {
        std::ofstream os("dft_idft.txt");
        output_plot_real(os, idft_data);
    }

    // 周波数スペクトルに変換
    constexpr auto spec_data = sprout::spectrum(begin(dft_data), end(dft_data), pit<spectrum_t>());
    {
        std::ofstream os("spec.txt");
        output_plot(os, spec_data);
    }
}

出力を gnuplot でグラフ化
- src.txt (入力ソース)

dft_idft.txt ([DFT→IDFT]で変換→逆変換したもの)

spec.txt (周波数スペクトル)

周波数スペクトルのピークが X 軸上の 10 のところにあるので、周波数 10 の正弦波であることがわかります。

変換→逆変換の結果が入力ソースと等しくなっているので、正しく相互変換されていることも確かめられました。

主な機能は、

　　　sprout::dft →DFT

　　　sprout::idft →IDFT(逆変換)

　　　sprout::spectrum →周波数スペクトルに変換

また、正弦波の生成のところで使っているのは、最近 Sprout に追加された Range アダプタです。

　　　adaptors::sinusoidal →正弦波を返す Range アダプタ。そのものを Range として扱うこともできる。

　　　adaptors::copied →他のコンテナへ暗黙の変換をする Range アダプタ。Pstade.Oven の copied と同じ。

もうすこし複雑な例も。

重ね合わせた正弦波の DFT

#include <iostream>
#include <fstream>
#include <sprout/array.hpp>
#include <sprout/pit.hpp>
#include <sprout/complex.hpp>
#include <sprout/numeric/dft.hpp>
#include <sprout/range/adaptor.hpp>

template<typename Elem, typename Traits, typename InputRange>
std::basic_ostream<Elem, Traits>& output_plot(std::basic_ostream<Elem, Traits>& os, InputRange const& range) {
    os << std::fixed;
    int x = 0;
    for (auto const& e : range) {
        os << x++ << ',' << e << '\n';
    }
    return os;
}
template<typename Elem, typename Traits, typename InputRange>
std::basic_ostream<Elem, Traits>& output_plot_real(std::basic_ostream<Elem, Traits>& os, InputRange const& range) {
    os << std::fixed;
    int x = 0;
    for (auto const& e : range) {
        os << x++ << ',' << real(e) << '\n';
    }
    return os;
}

int main() {
    using namespace sprout;

    constexpr std::size_t size = 256;
    typedef array<complex<double>, size> container_t;
    typedef array<double, size> spectrum_t;

    // 周波数 10, 25, 35 の正弦波の重ね合わせを生成
    constexpr container_t src_data
        = adaptors::sinusoidal(10. / size, 10000.)
        | adaptors::transformed(adaptors::sinusoidal(25. / size, 5000.), plus<double>())
        | adaptors::transformed(adaptors::sinusoidal(35. / size, 7500.), plus<double>())
        | adaptors::copied
        ;
    {
        std::ofstream os("src.txt");
        output_plot_real(os, src_data);
    }

    // DFT(離散フーリエ変換)
    constexpr auto dft_data = sprout::dft(begin(src_data), end(src_data), pit<container_t>());

    // IDFT(逆離散フーリエ変換)
    constexpr auto idft_data = sprout::idft(begin(dft_data), end(dft_data), pit<container_t>());
    {
        std::ofstream os("dft_idft.txt");
        output_plot_real(os, idft_data);
    }

    // 周波数スペクトルに変換
    constexpr auto spec_data = sprout::spectrum(begin(dft_data), end(dft_data), pit<spectrum_t>());
    {
        std::ofstream os("spec.txt");
        output_plot(os, spec_data);
    }
}