【C++】【Eigen】フーリエ変換を利用したスムージング
はじめに
機器から出力された信号データに対してフーリエ変換を実施すると、周波数領域のデータに変換できます。この変換したデータに対して、低周波数成分のみ残して逆変換を行うと、ノイズを低減したデータを得ることができます。今回は、C++とEigenを使用して、ノイズを含む信号データをスムージングしてみます。
Savitzky-Golay法やWhittaker Smootherによるスムージングについてはこちらをご覧ください。
https://izadori.net/math-savitzky-golayEigenによる高速フーリエ変換のやり方についてはこちらをご覧ください。
Eigenによる実装
データ
今回スムージングを行うデータをこちらに用意しました。疑似データで分光スペクトル等を想定しています。
このデータファイルがプログラムを同じフォルダに置いてあるという前提で、データを読み込み処理します。
フーリエ変換
読み込んだデータをフーリエ変換します。変換にはEigen::FFT<>::fwd()を呼び出します。
Eigen::VectorXd y_vector;
Eigen::VectorXcd y_fft;
Eigen::FFT<double> fft;
// y_vectorには読み込んだデータが入っているものとします
fft.fwd(y_fft, y_vector);フーリエ変換した結果は以下のようになります。今回、データの間隔が0.01ごととなっているので、表示範囲を1 / 0.01 = 100の半分の50までとします。
実数部分
虚数部分
高周波成分の除去
フーリエ変換したデータから高周波成分を除去します。今回は、図から（適当ですが）空間周波数15をしきい値として、それより大きい空間周波数の値を0に置き換えることとします。実数部分、虚数部分ともに0に置き換えます。
// dxはデータの間隔
double dx = 0.01;
// data_xはｘ軸のデータが入っているものとします
for(unsigned int i = 0; i < data_x.size(); i++){
// 空間周波数の値を計算
double freq = i / (dx * data_x.size());
if(freq > 15){
y_fft(i) = std::complex<double>(0, 0);
}
}逆変換
高周波成分を除去したデータに対して逆変換を実施します。これはEigen::FFT<>::inv()を呼び出すだけです。
Eigen::VectorXd y_ifft;
fft.inv(y_ifft, y_fft);全コード
全コードを以下に示します。データを読み込んで、元のデータとスムージング後のデータ、周波数領域へ変換したデータを標準出力へ出力します。
上の説明とコードを多少変えているところがあります。
#include <complex>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
#include <Eigen/Eigen>
#include <unsupported/Eigen/FFT>
int main()
{
std::vector<double> data_x, data_y;
std::ifstream fin;
std::string s;
double x, y;
double dx = 0.01;
fin.open("./test_data.txt");
if(!fin)
return 1;
fin >> s >> s;
while(1){
fin >> x >> y;
if(fin.eof()){
break;
}
data_x.push_back(x);
data_y.push_back(y);
}
fin.close();
Eigen::VectorXd x_freq(data_x.size()), y_vector, y_ifft;
Eigen::VectorXcd y_fft, y_fft_tmp(data_y.size());
Eigen::FFT<double> fft;
double threshold = 15;
y_vector = Eigen::VectorXd::Map(data_y.data(), data_y.size());
fft.fwd(y_fft, y_vector);
for(unsigned int i = 0; i < data_x.size(); i++){
x_freq(i) = i / (dx * data_x.size());
if(x_freq(i) > threshold){
y_fft_tmp(i) = std::complex<double>(0, 0);
}
else{
y_fft_tmp(i) = y_fft(i);
}
}
fft.inv(y_ifft, y_fft_tmp);
for(unsigned int i = 0; i < data_y.size(); i++){
std::cout << data_x[i] << '\t' << data_y[i] << '\t' << y_ifft(i) << '\t';
std::cout << x_freq(i) << '\t' << y_fft(i).real() << '\t' << y_fft(i).imag() << std::endl;
}
return 0;
}動作を確認した環境
- OS
- Windows 10
- Microsoft C/C++ Optimizing Compiler 19.28.29337 for x64
- MinGW-W64 GCC version 8.1.0
- Clang version 11.0.1
- Linux
- GCC version 9.3.0
- Windows 10
- Eigen version 3.3.9
結果
上記コードによる結果を以下に図示してみました。グラフはJavaScriptのChart.js+chartjs-plugin-zoomを使っています。$x$軸、$y$軸上でマウスのホイールを回転させることで拡大・縮小を、ドラッグでスライドさせることが可能です。
平滑化後の元のデータはこちらにあります。
まとめ
フーリエ変換を使用したスムージングを、C++とEigenを使用して実装してみました。Eigenを使うと、比較的簡単にフーリエ変換によるノイズ低減を実現できます。