【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]に置換されます.