2012/01/07

FIR LPF/HPFの係数プログラム Java版

昨年JavaScriptで書いたものをJavaへ移植してみた。意外とフィルタ係数関係のアクセスが多いので、2012年はデジタルフィルタ関係、もしくは音響関係を少し充実させようかな? 下のFIRなども、原理や用途など後々解説します。個人的には図書館などで資料を少しあさったり、大学の先生のホームページなどを参考にプログラムしているのだが、説明が難解という印象があるので、中学生ぐらいでも分かるように、まとめていこうかと思っている。ちなみに私は元々文系寄りで、情報等の専門教育は受けたことはないです。すべて独学なので間違っていることもあると思います。たまに質問がくるので、ここでお断りしておきます。学生の方は先生に聞いてください。

プログラムの内容はFIRのLPF/HPFで、こちらのページのJavaScript版 (120111窓関数を追加) と基本的には同じ。機能として増やしたのは窓関数をHammingだけでなくHanningも加えたことと、フィルタ係数すべてを配列に格納したところ。実行すると、パラメータの値に応じて、プロンプトにフィルタ係数を出力するようになっているので、自由に書き換えてみてください。そのフィルタ係数をAudacity等に読み込んでみると、下のような波形、つまりインパルス応答になる。
下はLPFの場合。わかる人には馴染み深く、そうでない人はさっぱりの世界。


下はHPF。デルタ関数みたいな地味なかたちだが微妙に波打っている。


参考までにデルタ関数は、1サンプルだけレベルが1になっている。いろいろ便利な関数で、重要なもの。


FIRの係数を求める式

n はタップのn0, n1, n2・・・で、h(n) はタップn番の時の係数hという意味。ωはカットオフ周波数。0~πの範囲で設定する。周波数(サンプリング周波数の1/2)の半分ならπ/2。1/3であればπ/3。上記はLPFであれば、わりとそのまま使える。HPFは、若干いじることで実現できる。下のソースを見ればどう計算しているか分かると思うけど、ちょっとゴチャゴチャしている。変数も減らして、シンプルにきれいに書き直せばよいのだが、とりあえず動くのでそっとしておく。出来た係数を見比べるとLPFとHPFでは中央値以外はプラスマイナスの違いしかない。昨年FIRのプログラムを作るとき、中央値の係数を求める方法が不明だった。参考にした本にはLPFもHPFも中央値が同じだったりして、あれ?という感じだった。中央値が同じになるのは係数が0.5の時だけだと思うのが・・・ 実際作ってみると、不具合でまくりだったので、自分で理屈から考えて値を決めることにした。結果的には音声ファイルを意図どおりに加工できているので、大きな間違いはないと思う。

窓関数
理想的には無限個のタップを使うことなのだが、現実的ではないので、必要に応じた有限個のタップで打ち切る。そうすると、どうしてもギッブズという凸凹が出来てしまう。そこで窓関数を係数に掛けることで、ギッブズを軽減させることができる。そのためFIRでは普通窓関数を使う。図は窓関数を使っていないのが左。カットした部分がボコボコしている。右は窓関数を使って平坦にしている。

窓関数には用途に応じてかなりの種類があるのだが、代表的な窓関数は以下の2つ。

Hamming
周波数の分解能が高い。音声関係では一番使われる窓関数だと思われる。式の表記は区間の指定によって、いろいろあるようだが、今回はプログラムに合わせて以下のように書いてみた。大文字のNはタップの総数で、小文字のnは各タップの係数になる。

下はHammingの波形を出力してAudacityで表示したところ。窓関数の表示方法は、いろいろあるようだが、波形はプラスマイナスが含まれるのが普通なので、このようにしてみた。窓関数を通すと前後が0に近づくように加工される。波形の作り方は、レベルが常時マックスな波形(N数)をHammingに通して作っている。

Hanning
小さな成分を検出しやすい。計測などに利用されている。

波形は前後が完全に0に収束しているのが分かる。

下がソースで、内容が分かりやすいように、係数を求める部分だけ。とりあえず単体で動くようにはしておいた。LPFとHPFは途中まで同じ計算なので、ひとつにまとめようかと思ったけど、見にくくなりそうだったので、別にしておいた。計算式も、なるべく素直にしているので、どうやって求めているか分かると思う。コメントも少し入れておいた。
public class  FIR_LowHigh {
 public static void main(String[] args) {
  
 /* パラメータ 求めたい値に書き換えて! */
  double sFreq= 44100; /* サンプリング周波数 */
  double cFreq = 5000; /* カットオフ周波数 */
  int tapn = 21; /* タップ数 奇数   */
  boolean lpfflag = true ; /* フィルタ切替LPF=true, HPF=false */
  int windowF = 0; /* 窓関数の切替 hamming=0, hanning=1 */
  
 /* 変数 */
  int tapM =((tapn-1)/2); /* タップ中央の位置 */
  double omegaLPF = Math.PI*(((cFreq)/(sFreq))*2);
  double omegaHPF = Math.PI*(1-(((cFreq)/(sFreq))*2));
  double hCenter = 1-(cFreq/(sFreq/2)); /* センターの係数 */
  double tapArray[] = new double [tapn]; /* タップ数 */
  double h = 0; /* 係数 */
  double wHamming = 0; /* ハミング窓初期化 */
  double wHanning = 0; /* ハニング窓初期化 */
  double hWinLPF = 0; /* LPF HPF用 */
  double hWinHPF = 0; /* HPF用 */

 /* LPF */
  if(lpfflag == true){ 
  for(int i=1; i <= tapM; i++){
 h = Math.sin(i*omegaLPF)/(Math.PI*i);
  /* 窓関数の切替 */
  switch(windowF){
  case 0://Hamming
  wHamming = 0.54+0.46*Math.cos((i*Math.PI)/(tapM));
  hWinLPF = h*wHamming;
  break;
  case 1: /* hanning */
  wHanning = 0.5+0.5*Math.cos((i*Math.PI)/(tapM));
  hWinLPF = h*wHanning;
  break;
  }
  tapArray[tapM+i] = hWinLPF; /* 中央から後ろへ */
  tapArray[tapM-i] = hWinLPF; /* 中央から前へ */
  }
  /* 中央値 */
  tapArray[tapM] = 1.0-hCenter; 
  }
 /* HPF */
  else{ 
  for(int i=1; i <= tapM; i++){
  /* 係数を求める */
 h = Math.sin(i*omegaHPF)/(Math.PI*i);
  /* 窓関数の切替 */
  switch(windowF){
  case 0://Hamming
  wHamming = 0.54+0.46*Math.cos((i*Math.PI)/(tapM));
  hWinLPF = h*wHamming;
  break;
  case 1:/* hanning */ 
  wHanning = 0.5+0.5*Math.cos((i*Math.PI)/(tapM));
  hWinLPF = h*wHanning;
  break;
  } 
  hWinHPF = hWinLPF*Math.pow(-1,i);/* HPF */
  /* 配列に入れていく */
  tapArray[tapM+i] = hWinHPF;
  tapArray[tapM-i] = hWinHPF;
  }
  /* 中央値 */
  tapArray[tapM] = hCenter; 
  }
  
 /* タップをプロンプトに出力 */
  for(int i=0; i<tapArray.length; i++){
  System.out.println(tapArray[i]+",");
  }
 }
}



sound programming 目次はこちら