サビツキーゴーレイの係数を求める@ C++/Eigen

昔Rubyで書いたサビツキーゴーレイ係数を求めるプログラムをC++に書き直してみた。

Rubyでは処理時間気にしてゴリゴリ各要素計算した部分があったけど、

読みにくいだけなので、書き直しました。

元々はhttp://www.empitsu.com/wp/?p=112の記事を参考に実装しましたが今回は昔Ruby で書いたコードを信じていの移植です。

使い方

#include "savi.h"
.......
Savi savi;
savi.init(何次関数を使うか,フィルタ長);
savi.coeff(0);//微分なしのFIRフィルタ係数
savi.coeff(1);//1次微分したときのFIRフィルタ係数

フィルタ長は指定した長さがnだとすると2n+1の長さのFIRフィルタ係数が得られます

以下ライブラリ本体(savi.h)

#ifndef _SAVI_H
#define _SAVI_H
//OriginalCode is 
//Copyright 2010 akira_you niconico-tech.
//Rewrite to c++ by akira noda 2014

#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wconversion"
#pragma clang diagnostic ignored "-Wsign-conversion"
#include <Eigen/Dense>
#pragma clang diagnostic pop
#include<cmath>
namespace SAVI{


template<typename T>
class Savi {
private:
int hlen;
int order;
public:
Savi(){
hlen=-1;
}
typedef Eigen::Matrix<T,-1,1> Vec;
typedef Eigen::Matrix<T,-1,-1> Mat;
Mat svm;//savizky golay matrix

void init(int in_order,int in_hlen){ //actual filter len is hlen*2+1
if(hlen==in_hlen && order==in_order)return;
hlen=in_hlen;
order=in_order;
assert(0<order);
assert(1<in_hlen);
assert(order <1+2*in_hlen -1 );
Mat xm(order+1,order+1);
for(int i=0;i<=order;i++)
for(int j=i;j<=order;j++){
T sum=0;
for(int a=-hlen;a<=hlen;a++){
sum+=(T)pow(a,i+j);
}
xm(j,i)=xm(i,j)=sum;
}
Mat xm_inv=xm.inverse();
Mat xx(order+1,hlen*2+1);
xx.row(0)=Vec::Ones(hlen*2+1);
for(int i=0;i<hlen*2+1;i++)xx(1,i)=i-hlen;
for(int di=2;di<=order;di++){
xx.row(di)=xx.row(di-1).array()*xx.row(1).array();
}
svm= xm_inv*xx;
}
Vec coeff(int diffOrder=0){
double a=1;
for(int i=2;i<=diffOrder;i++)a*=i;
return a*svm.row(diffOrder);
}
};
};
#endif

ちなみに、フィルタの中心部分を求めるように処理をしていますが、以下の２つを書き換えると中心位置をずらせます。（0が入る場所が中心）

中心位置は0.5だけずらすという事もできて、そうすると内挿ができます。（参考http://d.hatena.ne.jp/akira_you/20110428）