【GAS】アセンブラで名前付きローカル変数を利用できるようにするm4マクロ

アセンブラ開発でしばしば使われるm4マクロを使って
高級言語のように名前付き変数を実現する方法を紹介します。以下のコードはLinux GCC専用です。

define(`tmpcount',0)
define(`setVars',`ifelse($#,1,,eval($#%2),0,`dnl
        define(`count$1',eval(tmpcount+$2)) dnl
        define(`tmpcount',eval(tmpcount+$2)) dnl
        define($1,`-count$1 (%ebp)') dnl
        setVars(shift($@))',`setVars(shift($@))') dnl
')

上記のマクロを定義して、以下のように、setvarsマクロによって変数の宣言を行います。

hello: .string "hello!"
.global main
main:
        pushl %ebp
        movl  %esp,%ebp
        subl  $8,%esp
        setVars(var1,4,var2,4)
        movl    $1,var1
        movl    $2,var2
        pushl $hello
        call puts
        leave
        ret

setVarsマクロはsetVars(変数1の名前,変数1のバイト数,変数2の名前,変数2のバイト数,...)
といった具合に宣言を行います。

【C言語】カーネル回帰分析の実装と検証

非線形な関係を持つデータ間をカーネル法をつかって回帰分析してみます。 また、後半で本手法の問題点も示します。

このようなデータを当てはめてみることにします。

このようになります。高い表現力があることがわかります。

コードはこちら。行列計算にはmeschというライブラリを使いました。

// kernelRegression
// x:data
// dimention: x's dimention
// y: label of x 
// num: number of data
// result: coeffience
double *kernelRegression(double **x,int dimention,double *y,int num,double gamma,double regilarization){
 VEC *label;
 MAT *gram,*ident,*inv;
 int i,j;
 double *result;

 ident = m_get(num,num);
 label = v_get(num);
 memcpy(label->ve,y,sizeof(double)*num);
 gram = m_get(num,num);
 ident = m_get(num,num);
 for(i=0;i<num;i++){
  for(j=0;j<num;j++){
   gram->me[i][j]=kernel(x[i],x[j],dimention,gamma);
  }
 }
 inv=m_add(gram,sm_mlt(regilarization,m_ident(ident),NULL),NULL);
 inv=m_inverse(inv,NULL); 
 VEC *coefficient=v_get(num);
 mv_mlt(inv,label,coefficient);
 result=malloc(sizeof(double)*num);
 memcpy(result,coefficient->ve,sizeof(double)*num);

 m_free(ident);
 m_free(gram);
 m_free(inv);
 v_free(label);
 v_free(coefficient);
 return result;
}}

カーネルがこちら。

double kernel(double *x1,double *x2,int num,double gamma){
 VEC *xVector,*yVector,*diffVector;
 int i;

 xVector=v_get(num);
 memcpy(xVector->ve,x1,sizeof(double)*num);
 yVector=v_get(num);
 memcpy(yVector->ve,x2,sizeof(double)*num);

 diffVector=v_sub(xVector,yVector,NULL);

 double sum=0.0;
 for(i=0;i<num;i++){
  sum+=diffVector->ve[i]*diffVector->ve[i];
 }

 v_free(diffVector);
 v_free(xVector);
 v_free(yVector);
 
 return exp(-gamma*sum);
}

全ソースはgithubに上げています。
https://github.com/y-mitsui/kernelRegression
このように高い表現力を持つカーネル回帰ですが、多項式の次数がサンプル数と同じなのでサンプル数が増えると計算量が増えます。
また、線形の回帰分析と同様、誤差関数に二乗誤差を用いているので外れ値の影響を受けやすく、外れ値に対して積極的に当てはめようと するので、他のデータの当てはまりが悪ることがあります。
これらの問題を抑制させる回帰分析がサポートベクトル回帰で、そちらのほうがよく使われるように思います。

【C言語】ミーンシフトやKDEのバンド幅を自動調整

カーネル密度推定(KDE)はデータからノンパラメトリック(特定の確率分布を想定しない)に確率密度を推定する方法です。
また、ミーンシフトはKDEの理論を利用して関数の極値を求める方法で、カーネル数を事前に設定しなくてもよいクラスタリングとしてよく利用されていると思います。
しかし、どちらもバンド幅を設定する必要があります。
そのバンド幅をデータから自動調整する方法の一つにPlugin法があります。
Plugin法はデータの確率分布が正規分布に従うことを想定してます。

double variance(double *data,int num){
 int i;
 double mean=0,sum=0;
 
 for(i=0;i<num;i++){
  mean+=data[i];
 }
 mean/=num;

 for(i=0;i<num;i++){
  sum+=pow(mean-data[i],2);
 }

 return sum/num;
}
//プラグイン法によるバンド幅
//data:データ配列
//num:配列の要素数
//戻り値：バンド幅
double calcBandWidth(double *data,int num){
 return 1.06*sqrt(variance(data,num))*pow(num,-0.2);
}

【C言語】線形予測分析でスペクトル包括を求める

線形予測分析を使って音声ファイルのスペクトル包括を求めてみました。
手順は次のとおりになります。
1.音声を読み込む
2.指定フレームを切り出す
3.プリエンファシスフィルタをかける
4.自己相関関数を求める
5.自己相関関数をもとに次数が低いARモデルを構築する
6.ARモデルでFIRフィルタの周波数応答を求める
※プリエンファシスフィルタをかけたあと窓関数をかける場合もあります。

//autocorr 自己相関関数
// x:データ配列
// len:配列の要素数
// lagSize:最大のラグ(1～lagSizeまでの相関を求める)
// 戻り値：自己相関関数
double *autocorr(const double *x,const int len,const int lagSize){
 double *r;
 int lag,j;

 r=calloc(1,lagSize*sizeof(double));
 
 for(lag=0;lag<lagSize;lag++){
  for(j=0;j<len-lag;j++){
   r[lag]+=(double)x[j] * x[j + lag];
  }
 }
 return r;
}

上記は自己相関関数を求める数式通りのコードですが、より高速にもとめるにはフーリエ変換でパワースペクトルを求めてから逆フーリエ変換を行うことで自己相関関数が得らるので普通はそちらを利用します。
n次ARモデルを構築する際、1次からn次までのARモデルが必要になることがあります。例えば、データの状態を説明するのに、最適な(?)次数を示すAICをもとめる場合などです。 その時、効率的に計算するのがレビンソン法です。n次ARを求める際にn-1次のARを利用して逐次的に計算しています。 次のコードがそれになります。

// LevinsonDurbin ARモデルを構築
// r:自己相関関数
// A:求められたAR係数(戻り値)
// lpcOrder:次数
// 戻り値：ノイズ
double LevinsonDurbin(const double *r,double **A,const int lpcOrder){
 double *Ak=calloc(1,(lpcOrder+1)*sizeof(double));
 int k,n,j;
 double *U,*V;

 Ak[0] = 1.0;
 Ak[1] = - r[1] / r[0];
 double Ek = r[0] + r[1] * Ak[1];
 double lambda = - r[1] / r[0];

 U=calloc(1,sizeof(double)*(lpcOrder+2));
 V=calloc(1,sizeof(double)*(lpcOrder+2));

 for(k=1;k<lpcOrder;k++){
  lambda = 0.0;
  for(j=0;j<k+1;j++){
   lambda -= Ak[ j ] * r[ k + 1 - j ];
  }
  lambda /= Ek;

  U[0]=1.0;
  for(j=1;j<k+1;j++)
   U[j]=Ak[j];
  U[j]=0.0;

  V[0]=0.0;
  for(j=k;j>0;j--)
   V[k-j+1]=Ak[j];
  V[k-j+1]=1.0;
  double *v=vectorScalarMult(V,k+2,lambda);
  Ak=vectorVectorAdd(U,k+2,v);

  Ek *= 1.0 - lambda * lambda;
 }

 *A=&Ak[0];
 return Ek;
}

続いて周波数応答を求めます。

//freqz FIRフィルタの周波数応答を求める
// b a :フィルタ係数
// nB: b配列の要素数
// nA: a配列の要素数
// N: 周波数の分割数（大きい値ほど細かく調べる）
double *freqz(const double *b,const double *a,const int nB,const int nA,const int N){
 double *ret=malloc(sizeof(double)*(N+1));
 int n,i;
 double N2pi=2*M_PI/N;
 for (n = 0; n < N + 1; ++n) {
  double complex z,x,numerator,denominator,t;
  x=-I*N2pi * n;
  z=cexp(x);

  denominator=numerator=0.0+I*0.0;
  for (i = 0; i < nB; ++i) {
   t=b[nB - 1 - i] * cpow(z, i);
   
   numerator +=t;
  }
  for (i = 0; i < nA; ++i) {
   t=a[nA - 1 - i] * cpow(z, i);
   denominator += t;
  
  }
  

  ret[n]=cabs(numerator / denominator);
 }
 return ret;
 
}

全ソースはgithubに上げています。
https://github.com/y-mitsui/lpcSpectrum

C言語でインデクサを実現するm4マクロ

C#のインデクサのようなものをC言語で実現するm4マクロを紹介します。
まだまだ改良の余地がありますが暫定的なものとして公開します。

define(`defObject',`
        setIndexer($2)
        $1 $2=NULL;
        $2=$1Init(shift($@));
');
define(`setIndexer',`define(`$1',
                `if'`else('
                        `$'`#'`,2,'
                        ``$1''`->a['`$'`1'`*'`$'`2'`],
                        ``$1'''
                `)'
        )
');

#include <stdlib.h>

typedef struct{
        int *a;
} testObj;
testObj *testObjInit(void *a,int n){
        struct testObj *result=malloc(sizeof(struct testObj));
        result->a=malloc(sizeof(int)*n);
 return result;
}
int main(void){
        defObject(testObj,objs,1024);    //①
        objs(1,3)+=1;//ここでインデクサのようなアクセスを実現している
 objs->a[2]=2;//普通にアクセスすることも出来る
}

①のdefObjectの第１引数は型名、第２引数は宣言する変数名です。
第３引数以降はtestObjInitに渡すデータです。
defObjectの中で[構造体名]Init関数を呼び出しています。(ここではtestObjInit)
 defObjectの中で呼び出されているsetIndexerマクロがインデクサを実現する機能になります。
setIndexerの中の ``$1''`->a['`$'`1'`*'`$'`2'`], はobjs(1,3)と呼び出されたとき、objs->a[1*3]に置換することを書いています。
 仮にsetIndexerを次のように書いた場合、

define(`setIndexer',`define(`$1',
                `if'`else('
                        `$'`#'`,1,'
                        ``$1''`->a['`$'`1'],
                        ``$1'''
                `)'
        )
');

objs(3)がobjs->a[3]に置換されます.