購読中です 読者をやめる 読者になる 読者になる

ボレロ村上 - ENiyGmaA Code

中3女子です。

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

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



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]で変換→逆変換したもの)




周波数スペクトルのピークが 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);
    }
}
    • 出力を gnuplot でグラフ化
      • src.txt (入力ソース)


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




これもやはりピークが 10, 25, 35 のところにあり、正しく変換されています。


ここでは、正弦波の合成に transformed Range アダプタを使っています。
   adaptors::transformed →Range の各要素をファンクタで変換する。1引数版と2引数版がある。


もし Range アダプタを使わなければ、一々テンポラリのコンテナを生成しなければならず非効率なので、
Range アダプタの有用性を示すケースであるとも言えます。


sprout::dft の実装は、離散フーリエ変換の数学的定義を素直にコードに落としたものです。
最初は FFT(高速フーリエ変換) のアルゴリズムを実装しようと思ったのですが、
FFT アルゴリズムは空間の再利用とループを必要とするため、constexpr で同様の実装をすることは、
不可能であるかまたは却って非効率な処理にならざるをえず断念しました。


以上のサンプルではコードによって生成した波形を解析していますが、
例のごとく #include でテキストを取り込めば外部のソースを解析することもできるでしょう。


さて、コンパイル時に DFT を走らせて何の役に立つかと言えば、
あなたのコーヒータイムが伸びることで人生にほんの少しゆとりが生まれるかもしれません。
もっと楽しそうな使い方があればぜひ教えてください。