FFT

components/fft.hpp · GitHub で見る

共通ページ

#include "lib/components/fft.hpp"

std::vector<std::complex<double>> a{1, 2, 3, 4};
fft(a);
fft(a, true);

操作	内容・返り値	計算量
fft(a,inverse = false)	複素数配列をその場で FFT。true は正規化済みの逆変換	\(O(n\log n)\)

a は std::vector<std::complex<Real>>、Real は浮動小数点型、長さは 2 の冪です。空配列では何もしません。入力を変換結果で上書きし、返り値は void です。順変換は正の偏角、逆変換は負の偏角を使い、逆変換には長さでの除算を含みます。浮動小数点の丸め誤差が生じます。

長さを \(n\) として時間は \(O(n\log n)\)、追加空間は \(O(1)\) です。

検証コード

• test/convolution.test.cpp
• test/convolution_fft.test.cpp
• test/convolution_light.test.cpp
• test/cyclic_convolution.test.cpp
• test/cyclic_convolution_light.test.cpp
• test/multiplicative_convolution.test.cpp
• test/multiplicative_convolution_light.test.cpp
• test/multivariate_convolution.test.cpp
• test/multivariate_convolution_light.test.cpp
• test/ordered_convolution.test.cpp
• test/ordered_convolution_light.test.cpp

コード

#ifndef FFT_HPP
#define FFT_HPP

#include <algorithm>
#include <bit>
#include <cassert>
#include <cmath>
#include <complex>
#include <cstddef>
#include <numbers>
#include <vector>

template <class Real>
void fft(std::vector<std::complex<Real>> &a, bool inverse = false) {
	std::size_t n = a.size();
	assert(!n || std::has_single_bit(n));
	for (std::size_t i = 1, j = 0; i < n; ++i) {
		std::size_t bit = n >> 1;
		for (; j & bit; bit >>= 1) j ^= bit;
		j ^= bit;
		if (i < j) std::swap(a[i], a[j]);
	}
	for (std::size_t width = 2; width <= n; width *= 2) {
		Real angle = (inverse ? -2 : 2) * std::numbers::pi_v<Real> / Real(width);
		std::complex<Real> step(std::cos(angle), std::sin(angle));
		for (std::size_t start = 0; start < n; start += width) {
			std::complex<Real> w(1);
			for (std::size_t j = 0; j < width / 2; ++j, w *= step) {
				auto x = a[start + j], y = a[start + j + width / 2] * w;
				a[start + j] = x + y;
				a[start + j + width / 2] = x - y;
			}
		}
	}
	if (inverse)
		for (auto &x : a) x /= Real(n);
}

#endif

Back to top page