チャンデバソフトを作ろう!

プログラムを作ろう!

趣味でスピーカーを作るのは楽しい! 同じように趣味のプログラミングも楽しいものです。 私もRaspberry Pi2を買って初めて簡単なサウンドプログラミングをしてみました。 そして思ったことは

サウンドプログラミングって意外と簡単じゃん!みんなやってみようよ!

というわけで今回は3wayのチャンデバソフト(-24dB/OctのLinkwitz–Riley デジタルフィルタ)をきわめていい加減に作ってしまおうという企画です。

環境はLinux、再生ソフトはMPDを想定してます。

Unix哲学

マキロイによるUnix哲学

一つのことを、うまくやれ。

協調して動くプログラムを書け。

標準入出力(テキスト・ストリーム)を扱うプログラムを書け。標準入出力は普遍的インターフェースなのだ。

標準入出力

標準入出力はテキストに使われることが多いのですが、バイナリに対しても使えます。実際にmpdやsoxは標準出力にデータを送ることが出来ます。 またalsa付属のaplayは普通ファイル名を与えて再生しますが、ファイル名を省略した場合には標準入力からのデータを再生します。

標準入出力に対応すればalsaの使い方を知らなくでもサウンドを扱うプログラムを作ることができます。

ただ問題点もあって単に生のデータが入出力されるだけで、データのフォーマットは別に管理しなければなりません。 つまり音声データを取り扱う場合にはサンプリング周波数・チャンネル数・bit深度などは別にコマンドライン引数などで与える必要があります。

データフォーマット

標準入出力に対する音声フォーマットは大抵16bit整数か32bit整数です。一部24bitのデバイスもあります。 数値がch数ぶん並んでいるだけです。

今回も面倒なので32bit入力、32bit出力の固定フォーマットとします。

プログラミング言語

ある程度スピードが必要なのでCかC++で。オブジェクト指向がカッコよさそうなので今回はC++に決定。 私はC++を使うのは始めてで、知識は入門サイトをいくつか覗いてみただけです。 初級入門書以下の知識で実用プログラムを作ろうという無茶な企画です。

2次IIRフィルタの計算

ブロック図省略(ググればでてくるので見てね)

係数b0,b1,b2,a0,a1,a2を掛け算・足し算するだけ。 フィードバックがあるので以前の入力の影響をうける。 状態を持った関数になる。 おおっ!オブジェクトとして表せば良さそう?

2次のIIRフィルタのクラスを作ってみる。

class IIR2  // Second-order IIR filter
{
  private:
    double b0,b1,b2,a1,a2;
    double x0,x1,x2;
  public:
    IIR2(){}
    IIR2(double *coef){
      b0=coef[0];b1=coef[1];b2=coef[2];a1=coef[3];a2=coef[4];
      x0=0;x1=0;x2=0;
    }   
    double f(double s){ 
      x2=x1; x1=x0; //shift xn
      x0=-x1*a1-x2*a2+s;
      return (x0*b0+x1*b1+x2*b2);
    }   
};

2次のIIRフィルタがあればパラメトリックイコライザや2次のバターワースのローパス・ハイパスフィルタが作れる。

4次Linkwitz–Rileyフィルタ

同じ特性の2次Butterworthフィルタを2つ直列に繋ぐと4次のLinkwitz–Rileyフィルタが出来る。

2次のIIRフィルタの係数を受け取りそれを直列2段がけするクラスを作る。 これに2次のバターワースの係数を与えればOK。

class LR4 //series same IIR2 filter
{
  private:
    IIR2 f1,f2;
  public:
    LR4(double *coef){
      f1=IIR2(coef);
      f2=IIR2(coef);
    }
    double f(double s){
      return f2.f(f1.f(s));
    }
};

フィルタ係数の計算

いつも通りGNU Octaveを使う。

関数 butter をつかう。 サンプリング周波数44100、クロス周波数300,2500Hzで計算

fs=44100;
[bl1,al1]=butter(2,300/fs*2) #ローパス
[bh1,ah1]=butter(2,300/fs*2,'high') #ハイパス
[bl2,al2]=butter(2,2500/fs*2) #ローパス
[bh2,ah2]=butter(2,2500/fs*2,'high') #ハイパス

今回は面倒なので計算した係数を直接ソースにコピペする。

とりあえず完成!

32bit2ch入力で32bit8ch出力、そのうち6chを使って3wayを構成。

こんな短いプログラムにライセンスをつける意味はないけどとりあえずGPLにしてみる。

amaiir3way.cpp

// Copyright (C) 2015 Amanogawa Audio Labo
//
// This program is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free Software
// Foundation; either version 3 of the License, or (at your option) any later
// version.
//
// This program is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License along with
// this program; if not, see <http://www.gnu.org/licenses/>.

// amaiir3way.cpp version 0.1

#include<cstdio>
#include<cstdlib>
using namespace std;
const int N=8192; //buffersize = N * ch (int_32)

// b0,b1,b2,a1,a2   (a0=1)  
double C44k_300l[]={4.43273022494514e-04,8.86546044989029e-04,4.43273022494514e-04,-1.939570207351671,0.941343299441649 };
double C44k_300h[]={0.970228376698330,-1.940456753396659,0.970228376698330,-1.939570207351671,0.941343299441649};
double C44k_2500l[]={0.0251761145544014,0.0503522291088028,0.0251761145544014,-1.503695341299221,0.604399799516827};
double C44k_2500h[]={0.777023785204012,-1.554047570408024,0.777023785204012,-1.503695341299221,0.604399799516827};

class IIR2  // Second-order IIR filter
{
  private:
    double b0,b1,b2,a1,a2;
    double x0,x1,x2;
  public:
    IIR2(){}
    IIR2(double *coef){
      b0=coef[0];b1=coef[1];b2=coef[2];a1=coef[3];a2=coef[4];
      x0=0;x1=0;x2=0;
    }
    double f(double s){
      x2=x1; x1=x0; //shift xn
      x0=-x1*a1-x2*a2+s;
      return (x0*b0+x1*b1+x2*b2);
    }
};

class LR4 //series same IIR2 filter
{
  private:
    IIR2 f1,f2;
  public:
    LR4(double *coef){
      f1=IIR2(coef);
      f2=IIR2(coef);
    }
    double f(double s){
      return f2.f(f1.f(s));
    }
};

int main()
{
  int32_t buf_in[N][2];
  int32_t buf_out[N][8]={};
  int i,ret_n;

  LR4 lowl(C44k_300l);
  LR4 lowr(C44k_300l);
  LR4 midl1(C44k_300h);
  LR4 midr1(C44k_300h);
  LR4 midl2(C44k_2500l);
  LR4 midr2(C44k_2500l);
  LR4 highl(C44k_2500h);
  LR4 highr(C44k_2500h);

  while(true){
    ret_n=fread(buf_in,sizeof(int32_t)*2,N, stdin);
    for(i=0; i<ret_n; i++){
      buf_out[i][0]=lowl.f(buf_in[i][0]);
      buf_out[i][1]=lowr.f(buf_in[i][1]);
      buf_out[i][2]=midl2.f(midl1.f(buf_in[i][0]));
      buf_out[i][3]=midr2.f(midr1.f(buf_in[i][1]));
      buf_out[i][4]=highl.f(buf_in[i][0]);
      buf_out[i][5]=highr.f(buf_in[i][1]);
    }
    fwrite(buf_out, sizeof(int32_t)*8, ret_n, stdout);
    if (ret_n != N)
      return 0;
  }
}

コンパイル

開発環境をインストール

sudo apt-get install build-essential

g++ -O2 amaiir3way.cpp -o amaiir3way

sudo cp amaiir3way /usr/bin/

使い方

mpd.confを編集

audio_output { 
  type  "pipe" 
  name  "amaiir" 
  command  "amaiir3way | aplay -r 44100 -f S32_LE -c 8 -D hw:1,0" 
  format  "44100:32:2" 
} 

出力デバイスは8chのを選んでね。(HDMIなど)

計算速度的には192000でも問題無いです(出力デバイスさえ対応していれば)

おわりに

わずか100行以下で実用になるチャンデバソフトが出来上がってしまいました。

俺の方がマシなプログラムを作れるぜ!というかたは是非作って見てください。

2015年8月22日