この記事はWebGL Advent Calendar 2017 7日目の記事です．

はじめに

この記事では，WebGLを利用して音響現象をシミュレートする記事です．
昨今，計算力が必要な分野というとディープラーニング，ブロックチェーンのマイニング，CGあたりが頭に浮かびますが，音響解析にも計算力が必要な分野があります．そこで，WebGLを使って音響シミュレータのひな形を作ってみようと思います．

TL;DR

3Dは扱いません
GPGPU的な話です
CPU版プログラムとGPU版プログラムで同じ演算をさせて，作り方を比較します

対象読者

シミュレータを作ってみたい人
HTML / JavaScriptの知識は持っている人
WebGLに興味はあるけど躊躇している人

取り扱う音響現象

音響現象と言っても，信号処理，建築音響，音声，超音波，聴覚，音響教育など幅広い分野に分かれています．ここでは，シミュレータのひな形を紹介することを念頭に置き，できるだけシンプルな演算として，点音源による音場を可視化するプログラムを作成します．取り扱う理論の選出基準は以下のとおりです．

過去の演算結果を参照しない
理論が複雑ではない

また，点音源による音場を計算する過程で以下の条件を付けています．

減衰は考えない
最初は1つの点音源を考える
次に点音源を増やし，干渉する様子を可視化する

点音源による音場の計算

強さUの点音源が距離r離れた点に作る音圧dPは以下の式で表されます．
ちょっと難しく思えるかもしれませんが，振動によって作られる波は，一次元で考えれば位相が時間とともに変化する正弦波です．
よって，点音源の位相と，点音源から受音点の距離によって波高dPが定まります．ωやρは後で消えるので，ピンと来ない人は無視しておいてください．

dP = j \omega \rho \frac{U}{4 \pi r}e^{-jkr}

複数の点音源によって生成される音場を計算する場合，それぞれの点音源からある点までの距離を元に波高を求め、それらの和を求めることで合成波高Pを算出します．

P = j \omega \rho \sum_{n=1}^N \frac{U}{4 \pi r_n}e^{-jkr_n}

数式の展開

先ほどの式にはeや虚数単位が出てきています。このままでは計算しづらいのでこの部分を、オイラーの公式を用いて展開します。

e^j\theta = \cos \theta + j \sin \theta

これを波高dPの式に代入すると以下のようになります。

dP = j \omega \rho \frac{U}{4 \pi r}\{\cos(-kr) +j\sin(-kr) \}

ここで、cos成分はsinの内部に位相の計算を入れることで消去します。すると、以下の式のようになります。

dP = j \omega \rho \frac{U}{4 \pi r}\{ j\sin(-kr + \phi) \}
= - \omega \rho \frac{U}{4 \pi r}\{ \sin(-kr + \phi) \}

プログラム上の数式

次に、プログラム上で不要な箇所を大胆に省略します。前提条件で、減衰は扱わないこととしていました。これにより、波高の高低がグラフィカルに表されれば良いことになります。どのみち見やすいように振幅を調整するので、前ページの最後の式から振幅に関係する部分をカットします。すると、以下のように単純な数式で表すことができます。

dP = \sin(-kr + \phi)

ここで、波長定数kは以下の式で表されます。

k=\frac{2\pi}{\lambda}

プログラム

ほぼ共通して使用するHTML

titleと読み込むJSファイル名は，次節から示すプログラムに合わせて変更してください．内容は，idがgraph0というcanvasを作っているだけです．

<!DOCTYPE HTML>
<html lang="ja-JP">
<head>
  <meta charset="UTF-8">
  <title>1次元の音波</title>
  <script src="1d.js"></script>
</head>
<body>
  <canvas class="myCanvas" id="graph0" width="512" height="512"></canvas>
</body>
</html>

1次元の振動の伝播

まずは，振動現象を確認するための1次元の音波の伝播シミュレータです．波高は折れ線グラフで描画しています．画面上に説明がありませんが，縦軸が波高で横軸が距離を表しています．

上の方で説明した数式を，ほぼそのままプログラムに落としてあります．1次元ですので，ループはx方向のみとなっています．

演算を行うためのφ（プログラム中ではちょっとした勘違いによりthetaと表記）の値を4.5°進めた上で，requestAnimationFrameを利用して次のフレームを計算・描画させています．波長などのパラメータは，変数に代入している箇所を変更することで調整するものとします．

1d.js

var theta = 0;

function calc() {
  let canvas = document.getElementById( 'graph0' );
  let ctx = canvas.getContext( "2d" );

  let lambda = 8.6; // 波長は8.6mm （音で言えば40kHz）
  let c_size = 512; // canvasのサイズ
  let w_number = 4; // フィールド内の波の数

  ctx.clearRect( 0, 0, 512, 512 );

  let k = 2.0 * Math.PI / lambda;

  ctx.beginPath();
  ctx.moveTo( 0, 200 );
  for( var x=0; x<c_size; x++ ) {
    let r = x * ( lambda * w_number ) / c_size; // 距離
    let amp = Math.sin( - k * r + Math.PI * 2.0 / 360.0 * theta );
    ctx.lineTo( x, 100 * amp + 200 );
  }
  ctx.stroke();

  theta += 4.5; // 1フレームで位相を4.5°進める
  requestAnimationFrame( calc );
}

window.addEventListener( 'load', function() {
  calc();
});

2次元の音波の伝播（音源が1つの場合）

続いて，2次元に拡張したものです．ループがX方向とY方向の2重になっています．表示は，ある空間を上から眺めているものです．波高は色の明るさで表しています．（黒＝波高が負，明るい色＝波高が正）

1次元のプログラムと比べると，距離の計算部分が若干異なりますが，波高（変数amp）の計算自体は同一です．その後，ampの値（-1〜+1）を，0〜255の値に変換し，pixelDataに代入しています．

one_source.js

let theta = 0;

function calc() {
  let canvas = document.getElementById('graph0');
  let ctx = canvas.getContext( "2d" );

  let lambda = 8.6; // 波長は8.6mm （音で言えば40kHz）
  let c_size = 512; // canvasのサイズ
  let w_number = 16; // フィールド内の波の数

  ctx.clearRect( 0, 0, 512, 512);

  let imageData = ctx.createImageData( 512, 512 );
  let pixelData = imageData.data;

  let k = 2.0 * Math.PI / lambda;
  let m_size = ( lambda * w_number ) / c_size; // 1ピクセルのサイズ

  for( let y=0; y<c_size; y++ ) {
    for( let x=0; x<c_size; x++ ) {
      let px = c_size / 2.0 - x;
      let py = c_size / 2.0 - y;
      let r = Math.sqrt( ( px * m_size * px * m_size ) + ( py * m_size * py * m_size ) );
      let amp = Math.sin( - k * r + Math.PI * 2.0 / 360.0 * theta );
      let wh = Math.floor( 127 + 126.0 * amp );
      if( wh<0 || 255<wh ) console.log( "？？？" );
      pixelData[ y * c_size * 4 + x * 4 + 0 ] = 0;	//R
      pixelData[ y * c_size * 4 + x * 4 + 1 ] = wh;	//G
      pixelData[ y * c_size * 4 + x * 4 + 2 ] = wh;	//B
      pixelData[ y * c_size * 4 + x * 4 + 3 ] = 255;	//a
    }
  }
  ctx.putImageData( imageData, 0, 0 );

  theta += 4.5; // 1フレームで位相を4.5°進める
  requestAnimationFrame( calc );
}

window.addEventListener( 'load', function() {
  calc();
});

WebGLで計算を行う際の基礎知識

WebGLでは頂点を計算するVertex Shaderと，色情報を計算するFragment Shaderを作成します．WebGLを利用するシミュレータのプログラムは，大まかに言うと以下の順序で流れていきます．

Vertex Shaderの記述とコンパイル
Fragment Shaderの記述とコンパイル
プログラムオブジェクトの生成（リンク）
Shaderに渡す変数の設定
板の貼付け・描画
次のフレームの処理

コンパイルやリンクの部分は決まり文句と思って構いません．

ちょっと嬉しく感じるのが，canvas要素の描画コンテキストを作成する部分で，燦然と輝く"webgl"の眩しさです．（あくまで個人の感想です）

  let canvas = document.getElementById( 'graph0' );
  let gl = canvas.getContext( "webgl" );

変数の受け渡し

シミュレータを作る上で外せないのが，JSからShaderに変数を渡す部分です．詳しいところは次節のプログラムを参照して欲しいのですが，getAttribLocationでShader内の変数名を指定し，単独の数値の場合（ここではtheta_nやdinstance_n）はvertexAttrib1fを使って受け渡します．該当箇所を抜粋します．

  var theta_n = gl.getAttribLocation(program, "theta_n");
  gl.vertexAttrib1f( theta_n, theta );

配列を渡す場合はgetAttribLocationで変数名を指定した後に，enableVertexAttribArrayで変数領域へのアクセスを許可し，vertexAttribPointerを使って要素数や型を指定しつつ値を受け渡します．（と言いつつ，このプログラムでは定数を渡しているので，削除可能です．すみません，消し忘れました）

  var vertex = gl.getAttribLocation(program, "vertex");
  gl.enableVertexAttribArray(vertex);
  gl.vertexAttribPointer(vertex, 2, gl.FLOAT, false, 0, 0);

変数は，Vertex Shaderが受け取ります．attributeというキーワードが付いている変数宣言の箇所がありますが，これがJSから変数を受け取ることを表しています．

また，変数宣言時にvaryingが付いているものは，Vertex ShaderからFragment Shaderに変数を渡すための書き方です．つまり，このVertex Shaderはvertexとtheta_nとdistance_nをJSから受け取り，thetaとdistanceという変数の内容をFragment Shaderに渡しています．（受け渡しはしていますが，今のところ使用していない変数も含まれています）

    "attribute vec2 vertex;",
    "attribute float theta_n;",
    "attribute float distance_n;",
    "varying float theta;",
    "varying float distance;",

Fragment Shaderにはvaryingというキーワードを伴う変数宣言がありますが，これらがProgramable Shaderから送られてきた変数を受け取る部分です．

    "varying float theta;",
    "varying float distance;",

Fragment Shaderでは，パラメータに従って演算を行い，各ピクセルの色を決定します．ここで，JS版のプログラムと異なるのが，「ピクセルごとの繰り返し処理を書かず，各ピクセルで行う計算のみを記述すれば良い」ことです．描画するピクセルの位置は自動的にgl_FragCoordに代入されています．

2次元の音波の伝播（音源が1つの場合）のWebGL版

続いて，GPUで音波の伝播を演算・描画するソースコードを示します．今回のプログラムでは，3D空間上に1枚の板を貼り付けて，その板のピクセルごとの色を計算させるというやり方を採用しています．よって，計算自体はFragment Shader内で全て行っています．本来であれば各パラメータはJSから受け渡したいところですが，サンプルコードを作っていて力尽きました．

できるだけJS版と対応するように記述しているので，比較してみてください．計算そのものにはあまり関係しませんが，JS版だと色を0〜255で指定していましたが，WebGL版だと0〜1で指定します．

one_source_gpu.js

let theta = 0;

function calc() {


  let lambda = 8.6; // 波長は8.6mm （音で言えば40kHz）
  let c_size = 512; // canvasのサイズ
  let w_number = 16; // フィールド内の波の数

  gl.clearColor( 1.0, 1.0, 1.0, 1.0 );
  gl.clear(gl.COLOR_BUFFER_BIT);

  let imageData = gl.createBuffer();
  gl.bindBuffer( gl.ARRAY_BUFFER, imageData );

  // バーテックス(頂点)シェーダー
  var vSource = [
    "precision mediump float;",
    "attribute vec2 vertex;",
    "attribute float theta_n;",
    "attribute float distance_n;",
    "varying float theta;",
    "varying float distance;",
    "void main(void) {",
    "gl_Position = vec4(vertex, 0.0, 1.0);",
    "theta = theta_n;",
    "distance = distance_n;",
    "}"
  ].join("\n");

  var vShader = gl.createShader(gl.VERTEX_SHADER);
  gl.shaderSource(vShader, vSource);
  gl.compileShader(vShader);
  gl.getShaderParameter(vShader, gl.COMPILE_STATUS);

  // フグメントシェーダー
  var fSource = [
    "precision mediump float;",
    "varying float theta;",
    "varying float distance;",
    "void main(void) {",
    "const float PI = 3.1415926535897932384626433832795;",
    "const float lambda = 8.6;",
    "const float c_size = 512.0;",
    "const float w_number = 16.0;",

    "const float k = 2.0 * PI / lambda;",
    "const float m_size = ( lambda * w_number ) / c_size;",

    "float px = c_size / 2.0 - gl_FragCoord.x;",
    "float py = c_size / 2.0 - gl_FragCoord.y;",
    "float r = sqrt( ( px * m_size * px * m_size ) + ( py * m_size * py * m_size ) );",
    "float amp = sin( -k * r + PI * 2.0 / 360.0 * theta );",
    "gl_FragColor = vec4( 0, (amp+1.0)/2.0, (amp+1.0)/2.0, 1.0 );",
    "}"
  ].join("\n");

  var fShader = gl.createShader(gl.FRAGMENT_SHADER);
  gl.shaderSource(fShader, fSource);
  gl.compileShader(fShader);
  gl.getShaderParameter(fShader, gl.COMPILE_STATUS);


  // プログラムオブジェクトの生成
  var program = gl.createProgram();
  gl.attachShader(program, vShader);
  gl.attachShader(program, fShader);
  gl.linkProgram(program);
  gl.getProgramParameter(program, gl.LINK_STATUS);
  gl.useProgram(program);


  // シェーダー側の変数をjs側から設定する
  var vertex = gl.getAttribLocation(program, "vertex");
  gl.enableVertexAttribArray(vertex);
  gl.vertexAttribPointer(vertex, 2, gl.FLOAT, false, 0, 0);

  var theta_n = gl.getAttribLocation(program, "theta_n");
  gl.vertexAttrib1f( theta_n, theta );

  var distance_n = gl.getAttribLocation(program,"distance_n");
  gl.vertexAttrib1f( distance_n, 1 );

  // 4点の座標をセット
  var vertices = [
    -1.0, 1.0,
    1.0, 1.0,
    -1.0, -1.0,
    1.0, -1.0
  ];

  // 描画する
  gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(vertices), gl.DYNAMIC_DRAW);
  gl.drawArrays(gl.TRIANGLE_STRIP, 0, vertices.length/2);

  theta += 4.5; // 1フレームで位相を4.5°進める
  requestAnimationFrame( calc );
}

window.addEventListener( 'load', function() {
  calc();
});

複数の点音源による音波の干渉シミュレータ

ここまでのプログラムでは計算量が少なく，JS版であっても最近のPCのCPUであれば十分実用的な速度で動作します．つまり，WebGLを利用する恩恵が受けられないテーマでした．そこで，点音源の数を増やして，音波が干渉する様子をシミュレートしましょう．

まずはJS版です．これまでと異なるのは，8つの点音源を想定し，ある点の波高を計算する際に点音源からの距離を計算→波高を計算→波高の和を計算→平均を取り波高としています．

ここまで計算量が増えると，最新のCPUを用いてもフレームレートが下がってきます．

many_source.js

let theta = 0;

function calc() {
  let canvas = document.getElementById('graph0');
  let ctx = canvas.getContext( "2d" );

  let lambda = 8.6; // 波長は8.6mm （音で言えば40kHz）
  let c_size = 512; // canvasのサイズ（単位：ピクセル）
  let w_number = 16; // フィールド内の波の数
  let s_number = 8; // 音源は8つ
  let interval = 16; // 音源間隔（単位：ピクセル）

  ctx.clearRect( 0, 0, 512, 512);

  let imageData = ctx.createImageData( 512, 512 );
  let pixelData = imageData.data;

  let k = 2.0 * Math.PI / lambda;
  let m_size = ( lambda * w_number ) / c_size;

  for( let y=0; y<c_size; y++ ) {
    for( let x=0; x<c_size; x++ ) {
      let sx = -interval * s_number/2.0 + interval/2.0;
      let amp = 0;
      for( let n=0; n<s_number; n++ ) {
        let px = c_size / 2.0 - x - sx;
        let py = c_size / 2.0 - y;
        let r = Math.sqrt( ( px * m_size * px * m_size ) + ( py * m_size * py * m_size ) );
        amp += Math.sin( - k * r + Math.PI * 2.0 / 360.0 * theta );
        sx += interval;
      }
      let wh = Math.floor( 127 + 126.0 * amp / s_number );
      if( wh<0 || 255<wh ) console.log( "？？？" );
      pixelData[ y * c_size * 4 + x * 4 + 0 ] = 0;	//R
      pixelData[ y * c_size * 4 + x * 4 + 1 ] = wh;	//G
      pixelData[ y * c_size * 4 + x * 4 + 2 ] = wh;	//B
      pixelData[ y * c_size * 4 + x * 4 + 3 ] = 255;	//a
    }
  }
  ctx.putImageData( imageData, 0, 0 );

  theta += 4.5; // 1フレームで位相を4.5°進める
  requestAnimationFrame( calc );
}

window.addEventListener( 'load', function() {
  calc();
});

複数の点音源による音波の干渉シミュレータ（WebGL版）

同じ計算をGPUで行ってみます．基本的な流れはこれまで説明してきたとおりです．GPUは単純な計算の繰り返しに向いているため，PCのみでなく，タブレットやスマートフォンでもほとんどフレームレートが一定（60FPS）のままアニメーションが進みます．

many_source_gpu.js

let theta = 0;

function calc() {
  let canvas = document.getElementById( 'graph0' );
  let gl = canvas.getContext( "webgl" );

  let lambda = 8.6; // 波長は8.6mm （音で言えば40kHz）
  let c_size = 512; // canvasのサイズ
  let w_number = 16; // フィールド内の波の数

  gl.clearColor( 1.0, 1.0, 1.0, 1.0 );
  gl.clear(gl.COLOR_BUFFER_BIT);

  let imageData = gl.createBuffer();
  gl.bindBuffer( gl.ARRAY_BUFFER, imageData );

  // バーテックス(頂点)シェーダー
  var vSource = [
    "precision mediump float;",
    "attribute vec2 vertex;",
    "attribute float theta_n;",
    "attribute float distance_n;",
    "varying float theta;",
    "varying float distance;",
    "void main(void) {",
    "gl_Position = vec4(vertex, 0.0, 1.0);",
    "theta = theta_n;",
    "distance = distance_n;",
    "}"
  ].join("\n");

  var vShader = gl.createShader(gl.VERTEX_SHADER);
  gl.shaderSource(vShader, vSource);
  gl.compileShader(vShader);
  gl.getShaderParameter(vShader, gl.COMPILE_STATUS);

  // フグメントシェーダー
  var fSource = [
    "precision mediump float;",
    "varying float theta;",
    "varying float distance;",
    "void main(void) {",
    "const float PI = 3.1415926535897932384626433832795;",
    "const float lambda = 8.6;",
    "const float c_size = 512.0;",
    "const float w_number = 16.0;",
    "const float s_number = 8.0;",
    "const float interval = 16.0;",

    "const float k = 2.0 * PI / lambda;",
    "const float m_size = ( lambda * w_number ) / c_size;",
    "float sx = -interval * s_number/2.0 + interval/2.0;",
    "float amp = 0.0;",
    "for( int n=0; n<int(s_number); n++ ) {",
    "float px = c_size / 2.0 - gl_FragCoord.x - sx;",
    "float py = c_size / 2.0 - gl_FragCoord.y;",
    "float r = sqrt( ( px * m_size * px * m_size ) + ( py * m_size * py * m_size ) );",
    "amp += sin( -k * r + PI * 2.0 / 360.0 * theta );",
    "sx += interval;",
    "}",
    "amp /= 8.0;",
    "gl_FragColor = vec4( 0, (amp+1.0)/2.0, (amp+1.0)/2.0, 1.0 );",
    "}"
  ].join("\n");

  var fShader = gl.createShader(gl.FRAGMENT_SHADER);
  gl.shaderSource(fShader, fSource);
  gl.compileShader(fShader);
  gl.getShaderParameter(fShader, gl.COMPILE_STATUS);


  // プログラムオブジェクトの生成
  var program = gl.createProgram();
  gl.attachShader(program, vShader);
  gl.attachShader(program, fShader);
  gl.linkProgram(program);
  gl.getProgramParameter(program, gl.LINK_STATUS);
  gl.useProgram(program);


  // シェーダー側の変数をjs側から設定する
  var vertex = gl.getAttribLocation(program, "vertex");
  gl.enableVertexAttribArray(vertex);
  gl.vertexAttribPointer(vertex, 2, gl.FLOAT, false, 0, 0);

  var theta_n = gl.getAttribLocation(program, "theta_n");
  gl.vertexAttrib1f( theta_n, theta );

  var distance_n = gl.getAttribLocation(program,"distance_n");
  gl.vertexAttrib1f( distance_n, 1 );

  // 4点の座標をセット
  var vertices = [
    -1.0, 1.0,
    1.0, 1.0,
    -1.0, -1.0,
    1.0, -1.0
  ];

  // 描画する
  gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(vertices), gl.DYNAMIC_DRAW);
  gl.drawArrays(gl.TRIANGLE_STRIP, 0, vertices.length/2);

  theta += 4.5; // 1フレームで位相を4.5°進める
  requestAnimationFrame( calc );
}

window.addEventListener( 'load', function() {
  calc();
});

あとがき

WebGLを用いると，環境構築の手間無しにGPUを利用した演算パワーを手に入れることができます．また，そのパフォーマンスも十分なものと言えます．

これまで，高速な演算環境が必要であり，GPUに興味はあるけれど手を出しづらいと感じていた方のお力となれば幸いです．

プログラム中で波長を8.6mmとしているのは，研究遂行上必要だったので，この記事のもととなっているプログラム中でハードコートしていました．よく考えたら波長と描画範囲と音源の間隔の関係を比で表せば良いので，気が向いたら修正してアップロードし直します．

参考

松田晃一, WebGL+HTML5 3DCGプログラミング入門, カットシステム, 2012年