EigenでFFTする

はじめに

Eigenという行列演算ライブラリーにFFTが入っているので試してみました。

FFT

EigenのFFTはkissfft、FFTW、Intel Math Kernel Libraryをバックエンドとしていて、デフォルトだとLGPLのkissfftがバックエンドになるみたいです。インターフェースがEigenになるので元のライブラリーより使いやすくなります。

実装

実装といってもすごく簡単でこれだけです。

#include <Eigen/Core>
#include <unsupported/Eigen/FFT>

void dofft(const std::vector<float>& time, 
               std::vector<std::complex<float>>& freq)
{
	Eigen::FFT<float> fft;

	fft.fwd(freq, time);
}

void doinv(const std::vector<std::complex<float>>& freq, 
                std::vector<float>& time)
{
	Eigen::FFT<float> fft;

	fft.inv(time, freq);
}

以前C#のMath.NETを使ってFFTをしましたので比べてみます。C++/CLIでdofwd()を呼ぶコードを書きます。

List<float>^ FFT::ProcessList(List<float>^ data)
{
	std::vector<float> time;
	for each (float x in data)
	{
		time.push_back(x);
	}

	std::vector<std::complex<float>> freq;
	dofwd(time, freq);

	List<float>^ processed_data = gcnew List<float>();
	for (auto x : freq)
	{
		processed_data->Add(sqrt(x.real()*x.real()+x.imag()*x.imag()));
	}

	return processed_data;
}

同じ値が得られました(*ﾟｰﾟ)

2次元の場合

せっかく行列演算できるEigenを使っているので、２次元のデータをFFTしてみます。
一発でできる関数はないみたいなので、
(c++ - How to use Eigen FFT with MatrixXf? - Stack Overflow)
1行ごとにFFTしていきます。
ソースはこんな感じです。

void dofwd2d(MatrixXf& image, MatrixXcf& freq)
{
	FFT<float> fft;

	MatrixXcf temp(image.rows(), image.cols());
	for (int i = 0; i < image.cols(); i++)
	{
		temp.col(i) = fft.fwd(image.col(i));
	}

	temp.transposeInPlace();
	for (int i = 0; i < temp.cols(); i++)
	{
		freq.col(i) = fft.fwd(temp.col(i));
	}
}

// 表示用に結果をずらす
void shift(MatrixXcf& freq)
{
	MatrixXcf temp(freq.rows(), freq.cols());

	int cx = freq.cols() / 2;
	int cy = freq.rows() / 2;// 追記:rowとcol逆
	temp.block(0, 0, cx, cy) = freq.block(cx, cy, cx, cy);
	temp.block(cx, cy, cx, cy) = freq.block(0, 0, cx, cy);
	temp.block(cx, 0, cx, cy) = freq.block(0, cy, cx, cy);
	temp.block(0, cy, cx, cy) = freq.block(cx, 0, cx, cy);

	freq = temp;
}

col方向に演算した後でrow方向に演算するんですが、うまくいかなかったのでtransposeしています。
(追記：freq.row(i) =fft.fwd(temp.row(i)).transpose()みたいにすれば良かった)
適当な画像をFFTしてみました。Eigenの形式だと触りやすいので、いろいろ便利そうです。

cv::Mat_<uchar> image = cv::imread("image_Lena512.png", 0);

MatrixXf mat_image;
cv::cv2eigen(image, mat_image);

MatrixXcf freq(mat_image.rows(), mat_image.cols());
dofwd2d(mat_image, freq);
shift(freq);

追記:画像を逆変換で元に戻す

フーリエ変換してから逆変換すると画像が元に戻る確認です。

#include <Eigen/Dense>
#include <unsupported/Eigen/FFT>
using namespace Eigen;

#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/core/eigen.hpp>

MatrixXcf fft2(MatrixXf& in)
{
  FFT<float> fft;

  MatrixXcf temp(in.rows(), in.cols());
  for (int i = 0; i < in.cols(); i++)
  {
    temp.col(i).noalias() = fft.fwd(in.col(i));
  }

  MatrixXcf freq(in.rows(), in.cols());
  for (int i = 0; i < temp.rows(); i++)
  {
    freq.row(i).noalias() = fft.fwd(temp.row(i)).transpose();
  }

  return freq;
}

MatrixXf ifft2(MatrixXcf& in)
{
  FFT<float> fft;

  MatrixXcf temp(in.rows(), in.cols());
  for (int i = 0; i < in.cols(); i++)
  {
    temp.col(i).noalias() = fft.inv(in.col(i));
  }

  MatrixXf time(in.rows(), in.cols());
  for (int i = 0; i < temp.rows(); i++)
  {
    time.row(i).noalias() = fft.inv(temp.row(i)).transpose().real();
  }

  return time;
}

int main()
{
  cv::Mat1b image = cv::imread("a.jpg", 0);

  MatrixXf t;
  cv::cv2eigen(image, t);

  cv::Mat resultimage;
  cv::eigen2cv(ifft2(fft2(t)), resultimage);

  cv::imshow("Original", image);
  cv::imshow("Result", resultimage/255);
  cv::waitKey(0);
}