【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