ロコリンの雑記

アニメ大好きニートのロコリンのブログ。2015年卒(修士)の社会人。学生時代(2010年)から続けてるブログなのでエントリによっては学生ブログと社会人ブログになっています。時系列から察して。
 
 
このブログについて
ブログ内検索
カテゴリ
プロフィール

ロコリン

Author:ロコリン
2018年5月だけニート(6月から会社員)。2015年3月まで大学院生でした。
趣味:アニメ/Twitter/ゲーム/ニコ動
今(2015年2月更新):プリキュア/プリパラ/アイカツ/ごちうさ/艦これ

外部リンク
Twitter

スポンサーサイト 

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

フーリエ解析 2: C++ で高速フーリエ変換 (FFT) 

1 まえがき

ブログの更新間隔は忙しさに比例する法則 (←経験則)。5月からアルバイトを始め、さらに同じく5月から外国語の文献を講読してスライドにまとめる必修科目に追われ、それ以外のことを何もしていませんでしたorz 久々のブログ更新です。

そろそろ高速フーリエ変換 (FFT) の自前ライブラリを用意しないとなぁ…と思っていました。しかし、FFT の原理は正直言ってまだ完全には理解していません。そのため、「高速フーリエ変換 C++」でググってみたのですが、すぐに理想的な実装が見つかりました。

めでたしめでたし。なのですが、とりあえず私も実装してみたいなぁと思って実装してみました。といっても、ほとんど上記ブログの実装と変わりありません。ついにネタ切れで他者ブログのパクりに走ってしまいましたorz 丸パクりだと意味がないので、自分なりに機能を追加しています。イテレータとコピーコンストラクタに対応しただけですけど。

2 フーリエ変換

2.1 フーリエ変換

フーリエ変換とは、複素関数 (または実関数) から別の複素関数への変換です。関数 \(f: \mathbb{C} \to \mathbb{C}\) のフーリエ変換を \(F: \mathbb{C} \to \mathbb{C}\) とすると、フーリエ変換は次の積分で表されます。ここで \(i^{2}=-1\) です。

\[F\left(\omega\right) = \int_{-\infty}^{\infty}f\left(t\right)e^{-i\omega t}dt\]

フーリエ変換には逆変換 (フーリエ逆変換) が存在し、次の積分によって \(F\left(\omega\right) \) から \(f\left(t\right)\) を得ることができます。 (例外はありますが細かいことは気にしない。)

\[f\left(t\right) = \frac{1}{2\pi}\int_{-\infty}^{\infty}F\left(\omega\right)e^{i\omega t}d\omega\]

この演算にどんな意味があるかというと、例えば \(f\left(t\right)\) が時間 \(t\) の関数 (時系列信号) であれば、\(F\left(\omega\right)\) は角周波数 \(\omega\) の関数 (スペクトル) になります。

2.2 離散フーリエ変換

離散フーリエ変換 (DFT) とは、離散信号に対して行うフーリエ変換みたいなものです。

時系列信号 \(f\left(t\right)\) のある時間区間 (半開区間) \(\left[t_{0}, t_{0}+\Delta t\right) = \left\{t \mid t_{0} \le t < t_{0}+\Delta t\right\}\) 内の \(N\) 点をサンプリングしたものを \(x_{n}\) とします (ここで \(n \in \left\{n \in \mathbb{N} \mid 0 \le n < N\right\}\))。

\[x_{n} = f\left(t_{0}+\frac{n}{N}\Delta t\right)\]

\(x_{n}\) の離散フーリエ変換 \(X_{k}\) は次のようになります (ここで \(k \in \left\{k \in \mathbb{N} \mid 0 \le k < N\right\}\))。

\[X_{k} = \sum_{n=0}^{N-1}x_{n}\exp\left(-i\frac{2\pi}{N}kn\right)\]

離散フーリエ逆変換 (IDFT) は次のようになります (ここで \(n \in \left\{n \in \mathbb{N} \mid 0 \le n < N\right\}\))。

\[x_{n} = \frac{1}{N}\sum_{k=0}^{N-1}X_{k}\exp\left(i\frac{2\pi}{N}kn\right)\]

2.3 高速フーリエ変換

高速フーリエ変換 (FFT) とは、離散フーリエ変換 (DFT) を高速に計算するアルゴリズムです。DFT の計算量は \(O\left(N^{2}\right)\) となり、大きなデータに対しては非常に時間がかかります。FFT は計算手順を工夫することによって計算量を \(O\left(N\log N\right)\) に下げました。ただし、FFT においては、サンプル点数 \(N\) は 2 のべき乗に制限されます (\(N \in \left\{2^{n} \mid n \in \mathbb{N}\right\}\))。


C++ における FFT の実装例を示します。参考文献 [1] のブログを参考にしました。

fft.hpp

/**
 * FFT
 * coder: キュアガウス
 * date: 2012/06/19
 * message: 高速フーリエ変換 (FFT) ライブラリです。
 **

 ====
 クラス fft::FftArray
 複素数固定長配列です。FFT ができます。

 -- 使用例
 fft::FftArray z(16);	// サイズ16の複素数配列
 for (int i = 0; i < z.size(); ++i) {
 	z[i] = sin(2. * M_PI * i / z.size());
 }
 z.fft();	// フーリエ変換
 for (fft::FftArray::iterator it = z.begin(); it != z.end(); ++it) {
 	std::cout << *it << std::endl;
 }
 z.ifft();	// フーリエ逆変換
 for (int i = 0; i < z.size(); ++i) {
 	std::cout << z[i] << std::endl;
 }

 */

#ifndef	FFT_HPP
#define	FFT_HPP

#include <complex>
#include <vector>
#include <algorithm>
#ifndef	M_PI
#define	M_PI	3.14159265358979323846
#endif	// M_PI

namespace fft {
	class FftArray {
		typedef std::complex<double> Complex;
		std::vector<Complex> m_x;

		static size_t next2n(size_t x)
		{	// x 以上の最小の 2^n (n は自然数) を返す
			size_t y = 1u;
			for (; y < x; y <<= 1u);
			return y;
		}

		void bitReverse()
		{	// 要素の位置を、添え字をビット反転した位置に移動
			for (size_t i = 0u; i < m_x.size(); ++i) {
				size_t j = 0u;
				size_t bL = 1u, bR = m_x.size() >> 1u;
				for (; bL < m_x.size(); bL <<= 1u, bR >>= 1u) {	// bit reverse
					if ((i & bL) != 0u) {
						j |= bR;
					}
				}
				if (i < j) {
					std::swap(m_x[i], m_x[j]);
				}
			}
		}

		void doFft(bool isInverse = false) {
			bitReverse();
			for (size_t m = 2u; m <= m_x.size(); m *= 2u) {
				const double arg = (isInverse ? 2. * M_PI : -2. * M_PI) / static_cast<double>(m);
				const Complex w(std::cos(arg), std::sin(arg));
				for (size_t i = 0u; i < m_x.size(); i += m) {
					Complex wj(1., 0.);
					for (size_t j = 0u; j < m / 2u; ++j) {
						const size_t a = i + j;
						const size_t b = i + j + m / 2u;
						const Complex t(wj * m_x[b]);
						m_x[b] = m_x[a] - t;
						m_x[a] = m_x[a] + t;
						wj *= w;
					}
				}
			}
		}

	public:
		typedef std::vector<Complex>::iterator iterator;

		FftArray(size_t n) :
			m_x(next2n(n))
		{
		}

		template <class Iterator>
		FftArray(Iterator first, Iterator last)
			: m_x(next2n(std::distance(first, last)))
		{
			for (iterator it = m_x.begin(); first != last; *it++ = *first++);
		}

		Complex& operator[](int index) {
			if (index < 0) {
				const int add = (-index / static_cast<int>(m_x.size()) + 1) * static_cast<int>(m_x.size());
				index = add + index;
			}
			return m_x[static_cast<size_t>(index) % m_x.size()];
		}

		inline iterator begin() {
			return m_x.begin();
		}

		inline iterator end() {
			return m_x.end();
		}

		inline size_t size() const {
			return m_x.size();
		}

		void resize(size_t s, Complex v = Complex()) {
			m_x.resize(next2n(s), v);
		}

		inline void fft() {
			doFft();
		}

		void ifft() {
			doFft(true);
			for (iterator it = begin(); it != end(); ++it) {
				*it /= static_cast<double>(m_x.size());
			}
		}
	};
}

#endif	// FFT_HPP

main.cpp

#include "fft.hpp"
#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <cstdlib>	// std::exit

static const char *USAGE_TEXT = "\
Usage\n\
\n\
Options\n\
 -i	Inverse FFT\n\
 -c	Complex input\n\
 -a	Absolute output\n\
 -r	Output only real number\n\
 -h	This page";

class Interpret {
	bool flagInverse;
	bool flagComplex;
	bool flagAbsolute;
	bool flagReal;

	fft::FftArray input() const
	{
		typedef std::complex<double> Complex;

		std::vector<Complex> v;
		if (flagComplex) {
			double re, im;
			while (std::cin >> re >> im) {
				v.push_back(Complex(re, im));
			}
		} else {
			double re;
			while (std::cin >> re) {
				v.push_back(Complex(re, 0.));
			}
		}

		return fft::FftArray(v.begin(), v.end());
	}

	template <class Iterator>
	void output(Iterator first, Iterator last) const
	{
		std::cout << std::setprecision(15);	// 表示桁数 (精度) を15桁に変更
		for (Iterator it = first; it != last; ++it) {
			if (flagAbsolute) {
				std::cout << abs(*it) << std::endl;
			} else if (flagReal) {
				std::cout << it->real() << std::endl;
			} else {
				std::cout << it->real() << '\t' << it->imag() << std::endl;
			}
		}
	}

public:
	Interpret(const int argc, const char *argv[])
		: flagInverse(false), flagComplex(false), flagAbsolute(false), flagReal(false)
	{
		const std::vector<std::string> args(argv, argv + argc);
		std::vector<std::string>::const_iterator it = args.begin();
		for (++it; it != args.end(); ++it) {
			const std::string str = *it;
			if (str.find("i") != std::string::npos) {
				flagInverse = true;
			}
			if (str.find("c") != std::string::npos) {
				flagComplex = true;
			}
			if (str.find("a") != std::string::npos) {
				flagAbsolute = true;
			}
			if (str.find("r") != std::string::npos) {
				flagReal = true;
			}
			if (str.find("h") != std::string::npos) {
				std::cout << USAGE_TEXT << std::endl;
				std::exit(0);
			}
		}
	}

	void operator()() const
	{
		fft::FftArray z(input());

		// FFT or IFFT
		if (flagInverse) {
			z.ifft();
		} else {
			z.fft();
		}

		output(z.begin(), z.end());
	}
};

int main(const int argc, const char *argv[])
{
	Interpret ipr(argc, argv);
	ipr();
	return 0;
}

脚注

  • \(\mathbb{N} = \left\{n \mid n \text{は自然数}\right\}\)

    ただし本ブログでは \(0 \in \mathbb{N}\) とします (0 は自然数)。0 を除いた自然数を \(\mathbb{N}^{*} = \mathbb{N}\setminus\left\{0\right\}\) とします。

  • \(\mathbb{C} = \left\{z \mid z \text{は複素数}\right\}\)

参考文献

  1. c++で高速フーリエ変換, fft - yattの日記
  2. 高速フーリエ変換 - Wikipedia
  3. 小野測器-FFTアナライザについて

関連記事

  1. フーリエ解析 1: フーリエ級数 (2011.02.25)
スポンサーサイト
コメント















 管理者にだけ表示を許可する

トラックバック
 
http://rexpit.blog29.fc2.com/tb.php/78-8a5f232c
最新記事
最新コメント
FC2カウンタ
欲しい
最近買ったもの
Amazon 検索
 
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。