OCamlで信号処理

はじめに

　OCamlで音声データにフィルタを掛けてリアルタイム再生するのを試します（多分こういうのはCかC++あたりでやるのが良いと思いますが）。

スペック

　スペックが高くないマシンでもリアルタイム再生出来るのを確認したいのでRaspberry Piで試します。

• マシン: Raspberry Pi (Model B)
• OS: RASPBIAN Debian Wheezy (September 2014)
• USB Audio I/F: BEHRINGER UCA202

使用するフィルタ

　今回はReaktorで作った次のローパスフィルタをOCamlで実装します。Reaktorは半額セールのときに買うのがおすすめです。

　各オブジェクトの入出力に数値以外の文字列が書かれているものがありますが、文字列が一致しているところは線で繋がっていると考えて良いです。文字列で接続を表現することによりスパゲッティにならないようにしています。"SR.R"はサンプリングレート(44100や48000)を表しています。

信号の種類について

　信号の種類は大きく分けて3つあると思います。

種類	説明
オーディオ信号	オーディオのサンプル毎に更新される信号
イベント信号	パラメータ値が変更されたときに更新される信号
定数信号	接続を変更しない限りは一定の値をとる信号

　信号の種類で結線に色付けしてみます。赤:オーディオ信号、紫:イベント信号、青:定数信号とします。

　今回は再生途中でパラメータの変更は無しとしたいので次のようにイベント信号も定数信号として考えることにします。

実装の方針

　オブジェクト毎にカレント値を持たせ、AudioInから開始して繋がるオブジェクトを次々と更新することによりAudioOutの値を算出する方法が考えられますが、その場合はカレント値をmutable(破壊的代入が可能)な変数にする方法しか思いつかなかったのであまりやりたくありません。今回はAudioOutを開始として値算出に必要な経路を探索してみます。

AuidoOutを親とするツリー構造

　経路はAuidoOutを親とするツリー構造のようなもので表現します。探索時の注意ですが**z^-1（単位遅延演算子）**が見つかったらその先は探索せずに止めます。z^-1は前回値を保持していますのでAuidoOut値の算出にはそれを使用します。

z^-1を親とするツリー構造

　AuidoOut値を算出するためには、z^-1もサンプル毎に漏れなく更新する必要がありますのでz^-1を親とするツリー構造も考えます。探索した結果、自分自身に到達することもあります。ツリー構造と言っておきながら途中で合流してるようにも見えますね。あまり気にしないことにします。

ツリーを表現する

　次のパラメータ付の再帰ヴァリアントでツリーを表現します。

抜粋

(* ツリー構造 *)
type tree =
| AudioIn            (* オーディオ信号入力 *)
| ZIn of int         (* 遅延オーディオ信号入力 *)
| Const of float     (* 定数 *)
| SampleRate         (* サンプリングレート *)
| Add of tree * tree (* 加算 *)
| Sub of tree * tree (* 減算 *)
| Mul of tree * tree (* 乗算 *)
| Div of tree * tree (* 除算 *)
| Tan of tree        (* 三角関数のtan *)
| AudioOut of tree   (* オーディオ信号出力 *)
| ZOut of tree       (* 遅延オーディオ信号出力 *)

　AudioOutとZOutを親としてツリー構造を表現します。単位遅延演算子は入力のZInと出力のZOutの2つに分けます。ZInには識別のためのindexを保持し（indexは適当にidでソートして昇順で付けます）、Constには定数を保持しておきます。数学関数は今回使用するTanのみを追加しています。

接続情報の作成

　再帰ヴァリアントでツリー構造を表現する前に元となる接続情報を作成します。各オブジェクトを識別するために適当に番号を付けます。

JSONで接続情報を記述する

　接続情報をハードコードするのもアレなのでJSONファイルを手作りします。各オブジェクトの演算と接続を記述していきます。

test.json

{
  "objects": [
    {"id": "obj1", "operator": "audioin"},
    {"id": "obj2", "operator": "const", "value": 400.0},
    {"id": "obj3", "operator": "const", "value": 0.707107},
    {"id": "obj4", "operator": "audioout", "inputs": ["obj26"]},
    {"id": "obj5", "operator": "z^-1", "inputs": ["obj24"]},
    {"id": "obj6", "operator": "z^-1", "inputs": ["obj27"]},
    {"id": "obj7", "operator": "const", "value": 3.14159},
    {"id": "obj8", "operator": "samplerate"},
    {"id": "obj9", "operator": "div", "inputs": ["obj7", "obj8"]},
    {"id": "obj10", "operator": "mul", "inputs": ["obj2", "obj9"]},
    {"id": "obj11", "operator": "tan", "inputs": ["obj10"]},
    {"id": "obj12", "operator": "const", "value": 1.0},
    {"id": "obj13", "operator": "div", "inputs": ["obj12", "obj3"]},
    {"id": "obj14", "operator": "add", "inputs": ["obj11", "obj13"]},
    {"id": "obj15", "operator": "mul", "inputs": ["obj11", "obj14"]},
    {"id": "obj16", "operator": "const", "value": 1.0},
    {"id": "obj17", "operator": "add", "inputs": ["obj16", "obj15"]},
    {"id": "obj18", "operator": "div", "inputs": ["obj11", "obj17"]},
    {"id": "obj19", "operator": "mul", "inputs": ["obj5", "obj14"]},
    {"id": "obj20", "operator": "add", "inputs": ["obj6", "obj19"]},
    {"id": "obj21", "operator": "sub", "inputs": ["obj1", "obj20"]},
    {"id": "obj22", "operator": "mul", "inputs": ["obj21", "obj18"]},
    {"id": "obj23", "operator": "add", "inputs": ["obj5", "obj22"]},
    {"id": "obj24", "operator": "add", "inputs": ["obj23", "obj22"]},
    {"id": "obj25", "operator": "mul", "inputs": ["obj23", "obj11"]},
    {"id": "obj26", "operator": "add", "inputs": ["obj25", "obj6"]},
    {"id": "obj27", "operator": "add", "inputs": ["obj26", "obj25"]}
  ]
}

　可読性は良くないですが我慢します。今回はイベント信号は定数信号として考えますので"obj2"と"obj3"は定数オブジェクトとし、適当な値を設定しました。またオーディオ信号入力とオーディオ信号出力の数は1個ずつに制限します。

operator	説明
"audioin"	オーディオ信号入力
"delay"	遅延オーディオ信号入力
"const"	定数
"samplerate"	サンプリングレート
"add"	加算
"sub"	減算
"mul"	乗算
"div"	除算
"tan"	三角関数のtan
"audioout"	オーディオ信号出力
"z^-1"	単位遅延演算子

JSONを読み込む

　JSONファイルを読み込んでツリー構造を作成していきます。

抜粋

(* ツリー作成 *)
let create_tree ~j:json terminal_id = (
    let rec create_tree2 id = (
        let obj = find_object_by_id ~j:json id in
        match (get_operator obj) with
        | "audioin"     -> AudioIn
        | "const"       -> Const(get_value obj)
        | "z^-1"        -> ZIn(z_index_of ~j:json obj)
        | "samplerate"  -> SampleRate
        | _ as operator ->
            (match (operator, get_inputs obj) with
            | ("add", x::y::[]) -> Add(create_tree2 x, create_tree2 y)
            | ("sub", x::y::[]) -> Sub(create_tree2 x, create_tree2 y)
            | ("mul", x::y::[]) -> Mul(create_tree2 x, create_tree2 y)
            | ("div", x::y::[]) -> Div(create_tree2 x, create_tree2 y)
            | ("tan", x::[])    -> Tan(create_tree2 x)
            | _                 -> raise NotSupportedJSONFormat
            )
    ) in match (get_operator terminal_id, get_inputs terminal_id) with
         | ("audioout", x::[]) -> AudioOut(create_tree2 x)
         | ("z^-1", x::[])     -> ZOut(create_tree2 x)
         | _                   -> raise NotSupportedJSONFormat
)

　パターンマッチング便利ですね。上記のコードには直接登場しませんがJSONファイル読み込みにyojsonとcoreライブラリを使用しています。ツリー構造を作成する前に無限ループにならないことをチェックする必要がありますが、今回は頑張りませんでした。

出力信号の算出

　AudioOutとZOutを更新するための関数を作成していきます。定数信号と演算は固定ですのでオーディオ入力信号a_inと遅延信号z_in1, z_in2,...を引数とする関数が作成出来れば良さそうです。

サンプル毎の処理イメージ

a_out  = f_a (a_in, z_in1, z_in2, ...);
z_out1 = f_z1(a_in, z_in1, z_in2, ...);
z_out2 = f_z2(a_in, z_in1, z_in2, ...);
...

　定数信号を適用しながら演算を合成した関数f_a, f_z1, f_z2,...を作成してサンプル毎にオーディオ入力信号・遅延信号を適用する方法を考えたのですが、上手く出来る気がしなかったので次のような制限を付け加えることにします。

*　　関数 f_a, f_z1, f_z2,... はそれぞれ a_in, z_in1, z_in2,... の線形結合で表現出来るものに制限する*

　この制限を付け加えれば、関数 f_a, f_z1, f_z2,... はベクトル(a_in, z_in1, z_in2,...)と適当なベクトルとの内積で表現出来ます。そうすると演算の合成は$\mathbb{R}^n$（$n$はオーディオ入力信号数＋遅延信号数）上のベクトルの合成と等価と見なすことが出来そうです。

信号の表現

　オーディオ入力信号・遅延信号は単位ベクトルだと思うことにします。

　信号タイプは次のように定義します。

抜粋

(* 信号タイプ *)
type signal =
| AudioSignal of (float list) (* オーディオ入力信号・遅延信号 *)
| ConstSignal of float        (* 定数信号 *)

演算の合成

　再帰処理で信号ベクトルを合成します。線形結合で表現出来るものに制限しますので、非線形な演算等は対象外とします。

抜粋

(* ベクトル作成 *)
let create_vector ~z_count:z_count terminal_obj = (
    let rec create_vector2 out = (
        match out with
        | AudioIn      -> AudioSignal(1.0::(zero_vector z_count))
        | ZIn(index)   -> AudioSignal(0.0::(unit_vector z_count index))
        | Const(value) -> ConstSignal(value)
        | SampleRate   -> ConstSignal(44100.0)
        | Add(x, y)    ->
            (match (create_vector2 x, create_vector2 y) with
            | (AudioSignal(a), AudioSignal(b)) -> AudioSignal(List.map2_exn ~f:(+.) a b)
            | (ConstSignal(a), ConstSignal(b)) -> ConstSignal(a +. b)
            | (AudioSignal(_), ConstSignal(_)) -> raise NotSupportedJSONFormat
            | (ConstSignal(_), AudioSignal(_)) -> raise NotSupportedJSONFormat
            )
        | Sub(x, y) ->
            (match (create_vector2 x, create_vector2 y) with
            | (AudioSignal(a), AudioSignal(b)) -> AudioSignal(List.map2_exn ~f:(-.) a b)
            | (ConstSignal(a), ConstSignal(b)) -> ConstSignal(a -. b)
            | (AudioSignal(_), ConstSignal(_)) -> raise NotSupportedJSONFormat
            | (ConstSignal(_), AudioSignal(_)) -> raise NotSupportedJSONFormat
            )
        | Mul(x, y) ->
            (match (create_vector2 x, create_vector2 y) with
            | (AudioSignal(_), AudioSignal(_)) -> raise NotSupportedJSONFormat
            | (ConstSignal(a), ConstSignal(b)) -> ConstSignal(a *. b)
            | (AudioSignal(a), ConstSignal(b)) -> AudioSignal(List.map ~f:(fun x->x*.b) a)
            | (ConstSignal(a), AudioSignal(b)) -> AudioSignal(List.map ~f:(fun x->a*.x) b)
            )
        | Div(x, y) ->
            (match (create_vector2 x, create_vector2 y) with
            | (AudioSignal(_), AudioSignal(_)) -> raise NotSupportedJSONFormat
            | (ConstSignal(a), ConstSignal(b)) -> ConstSignal(a /. b)
            | (AudioSignal(a), ConstSignal(b)) -> AudioSignal(List.map ~f:(fun x->x/.b) a)
            | (ConstSignal(_), AudioSignal(_)) -> raise NotSupportedJSONFormat
            )
        | Tan(x) ->
            (match (create_vector2 x) with
            | AudioSignal(_) -> raise NotSupportedJSONFormat
            | ConstSignal(a) -> ConstSignal(tan a)
            )
        | AudioOut(x) -> create_vector2 x
        | ZOut(x)     -> create_vector2 x
    ) in match (create_vector2 terminal_obj) with
         | AudioSignal(a) -> a
         | ConstSignal(_) -> raise NotSupportedJSONFormat
)

リアルタイム再生

　最後にWAVファイルの音声データにフィルタを掛けてリアルタイム再生します。WAVファイル読み込みは自前で実装しましたが、運良くALSAのAPIが叩けるライブラリがありました。では再帰で回してきます。

抜粋

let a_vector, z_matrix = get_matrix "test.json"
and wav_handle = In_channel.create ~binary:true "test.wav" in
let wav_frames = get_frames wav_handle     (* フレーム数取得 *)
and alsa_pcm = open_default_pcm            (* オーディオデバイスを開く *)
and alsa_buffer = String.create (2*2*1024) (* PCMデータ転送用バッファ *)
and wav_buffer = Buffer.create (2*2*1024)  (* WAVデータ読み込み用バッファ *) in

let rec audio_process remain_frames z_in_ch1 z_in_ch2 = (
    (* WAVファイルからデータを1ブロック読み込み *)
    Buffer.clear wav_buffer;
    let block_size = if remain_frames > 1024 then 1024 else remain_frames in
    Buffer.add_channel wav_buffer wav_handle (block_size*4);

    (* ブロック処理 *)
    let rec audio_process_block alsa_frame z_in_ch1 z_in_ch2 = (
        if alsa_frame < block_size then (
            (* 1フレーム取得 *)
            let a_in_ch1, a_in_ch2 = read_wav_data wav_buffer alsa_frame in
            let in_ch1 = a_in_ch1::z_in_ch1
            and in_ch2 = a_in_ch2::z_in_ch2 in

            (* ch1更新 *)
            let a_out_ch1
                = List.fold2_exn ~f:(fun a x y->a+.x*.y) ~init:0.0 in_ch1 a_vector
            and z_out_ch1
                = List.map ~f:(List.fold2_exn ~f:(fun a x y->a+.x*.y) ~init:1.0e-100 in_ch1) z_matrix

            (* ch2更新 *)
            and a_out_ch2
                = List.fold2_exn ~f:(fun a x y->a+.x*.y) ~init:0.0 in_ch2 a_vector
            and z_out_ch2
                = List.map ~f:(List.fold2_exn ~f:(fun a x y->a+.x*.y) ~init:1.0e-100 in_ch2) z_matrix in

            (* PCMデータ転送用バッファ書き込み *)
            write_alsa_buffer alsa_buffer alsa_frame a_out_ch1 a_out_ch2;

            (* next frame *)
            audio_process_block (alsa_frame + 1) z_out_ch1 z_out_ch2
        )
        else z_in_ch1, z_in_ch2 (* end block *)
    ) in
    let z_out_ch1, z_out_ch2 = audio_process_block 0 z_in_ch1 z_in_ch2 in (* start block *)

    (* PCMデータ転送 *)
    ignore (Alsa.Pcm.writei alsa_pcm alsa_buffer 0 block_size);

    if remain_frames > 1024 then audio_process (remain_frames - 1024) z_out_ch1 z_out_ch2 (* next block *)
) in
let z_in_init = zero_vector (List.length z_matrix) in
audio_process wav_frames z_in_init z_in_init; (* start *)

　1ブロックあたり1024フレームとしてブロック単位でALSAにPCMデータを転送します。1.0e-100という数が登場しますが非正規化数が発生するのをソフトウェア側で防止するためのものです。あとPCMデータ転送で渡しているバッファは初回に作成したものに破壊的代入を行って使い回していますが、そこはご了承ください。

• ソースコード