5
インパルス応答の測定
建築音響の分野，あるいは騒音制御の分野において，インパルス応答を測定して種々の評価を行なうこ とが多い。これは下図に示すように，インパルス信号がすべての周波数範囲において均等なエネルギー を持ち，これを放射して得られた応答が，原理的には音源と受音点の間の音響的情報のすべてを含むか らである。ここでは，最近良く用いられる測定手法の原理について学ぼう。なお，出展，参考文献は以 下に示すとおりである。 ജ ജ ࠗࡦࡄ࡞ࠬାภ ࠗࡦࡄ࡞ࠬࠬࡐࡦࠬ Time Time Frequency Le ve l Frequency Le ve l ᵄᢙಽᨆ ᵄᢙಽᨆ ቶౝߩ㧘ᱷ㗀ߥߤߩ ޟᔕ╵ޠ߇⇥ߺㄟ߹ࠇࠆ 㖸Ḯ߆ࠄฃ㖸ὐ߹ߢߩ㧘ߔߴߡߩ㖸㗀⊛ᖱႎࠍߺ߁ࠆ߽ߩ Fig. 5.1: インパルス応答を用いた測定の概念参考文献
[1] 橘 秀樹，矢野博夫，環境騒音・建築音響の測定，コロナ社，2004. [2] P. A. Nelson, S. J. Elliott, Active Noise Control, Academic Press, 1992.[3] M. R. Schroeder, Integrated-impulse method measuring sound decay without using impulses, J. Acoust. Soc. Am. 66 (2), pp. 497–500, 1979.
[4] D. D. Rife, Transfer-function measurement with maximum-length sequences, J. Audio Eng. Soc. 37 (6), pp. 419–444, 1989. [5] 城戸健一,ディジタルフィルタの原理(その4),日本音響学会誌, 43 (7), pp. 538–543, 1987. [6] 鈴木陽一,浅野 太,曽根敏夫,音響系の伝達関数の模擬をめぐって(その2),日本音響学会誌, 45 (1), pp. 44–50, 1989. [7] http://tosa.mri.co.jp/sounddb/tsp/tsp design.html [8] 佐藤史明,室内音響インパルス応答の測定技術,日本音響学会誌, 58 (10), pp. 669–676, 2002. [9] 佐藤史明, Swept-Sine法に基づく音響伝搬測定,日本音響学会誌, 63 (6), pp. 322–327, 2007. [10] 藤本卓也,低域バンドでのSN改善を目的としたTSP信号に関する検討,日本音響学会講演論文集, pp. 433-434, 1999,9.
5.1 インパルス信号の生成
インパルス信号は本当にすべての周波数成分を含むのか，それを簡単なプログラムを使って確かめて みよう。三角関数 cos(ωt) を，ω を変化させながら加算していこう。どうなるか？ 例えば，周波数を 1 Hz から 100 Hz まで変化させ，時刻 t を −0.5 秒から +0.5 秒までプロットして みよう。なお，cos sum.m というプログラムを利用してよい。 1) 加算する周波数を増やしたり減らしたりすると，加算した波形がどのようになるのか，調べて みよう。 2) 計算する時間範囲を変化させてみて，どのような波形が得られるのか確かめよう。 上記のプログラムを実行して分かるように，時刻 t = 0 での位相が揃った，同振幅で色々な周波数の 正弦（余弦）波を加算することで，インパルス信号に近いものが生成されることがわかる。つまり，こ れを非常に細かく加算する，つまり積分することで，インパルス信号を生成することが可能となる。こ のことから，インパルス信号は非常に広い範囲の周波数帯域を満遍なく含むことが理解できるであろう。 このインパルス信号を音場に放射した場合の応答を求める方法はいくつか存在する。最も直感的な方 法は，本当にインパルス的な信号を放射することである。しかし，一般的にはこれは難しい。スパーク 放電や風船の破裂音を使う方法も試みられているが，スペクトルの幅が十分でなかったり，あるいはス ピーカを使う測定では，その動特性や動作範囲の限界などで十分なパワーを供給できなくなることなど が理由である。5.2 ランダム信号を用いた測定法
まずは，一様乱数やガウス乱数を用いて測定する方法を紹介する。信号の継続時間はインパルス信号 に比べて桁違いに長いので，S/N の良い，つまり精度の良い測定が期待できるはずである。 下図の様に，システムへの入力信号を x(t)，出力信号を y(t) と表したとき，両者は以下の畳み込み の関係で結び付けられるとしよう。 y(t) = h(t) ∗ x(t) = ∞ −∞h(t)x(t − τ )dτ (1) システムのインパルス応答 h(t) を求めることがここでの問題である。x(t)
h(t)
y(t)
ജ ജ 㖸႐ߥߤߩࠪࠬ࠹ࡓ Fig. 5.2: 入力 x(t), システム h(t), 及び出力 y(t) の関係。5.2.1 直接相関法 入力信号 x(t) の自己相関関数および出力信号 y(t) との相互相関関数はそれぞれ以下のように表さ れる。 Rx(τ ) = E[x(t)x(t + τ )] = limT →∞ 1 T T 0 x(t)x(t + τ )dt (2) Rxy(τ ) = E[x(t)y(t + τ )] = lim T →∞ 1 T T 0 x(t)y(t + τ )dt (3) ただし，E[x(t)x(t + τ )] は x(t)x(t + τ ) の期待値（時間平均）を求めるという意味である。 ここで，Rxy(τ ) を以下のように変形する。 Rxy(τ ) = E x(t) ∞ −∞h(τ1)x(t + τ − τ1)dτ1 = ∞ −∞h(τ1)E [x(t)x(t + τ − τ1)dt] dτ1 = ∞ −∞h(τ1)Rx(τ − τ1)dτ1 = h(τ ) ∗ Rx(τ ) (4) つまり，入力信号 x(t) と出力信号 y(t) の相互相関関数 Rxy(τ ) は，入力信号の自己相関関数 Rx(τ ) とシステムのインパルス応答 h(τ ) の畳み込みで表される。 ここで，入力信号として自己相関関数がデルタ関数であるような信号を用いれば，自己相関関数を求 めることでそのままシステムのインパルス応答が求められることになる。つまり，Rx(τ ) = δ(τ ) である 関数を用いればよいわけである。この目的のために，ホワイトノイズなどを用いることができる。 この方法による測定で，外部から無相関のノイズが混入した場合のことを考察する。ノイズ込みの出 力波形を y(t) = y(t) + n(t) とすると， Rxy(τ ) = E[x(t)y(t + τ )] = E[x(t){y(t + τ ) + n(t + τ )}] = E[x(t)y(t + τ )] + E[x(t)n(t + τ )] (5) 一般的に右辺第２項はゼロと考えられるために， Rxy(τ ) = E[x(t)y(t + τ )] = Rxy(τ ) となり，原理的にこの方法はノイズの影響を受けないことが分かる。
5.2.2 クロススペクトル法 自己相関関数 Rx，相互相関関数 Rxy をフーリエ変換することで，それぞれパワースペクトル密度関 数，クロススペクトル密度が得られる。これを式で書くと，以下のようになる。 Sx(f ) = ∞ −∞Rx(τ )e −jωτdτ (6) Sxy(f ) = ∞ −∞Rxy(τ )e −jωτdτ (7) (8) 更に，時間領域の畳み込みの関係が周波数領域では積となる合成積則より，Eq. (4) は以下と等価である。 Sxy(f ) = H(f ) · Sx(f ) (9) ここに H(f ) はインパルス応答の周波数領域表現（伝達関数）である。したがって， H(f ) = Sxy(f ) Sx(f ) (10) を求め，これを逆フーリエ変換することでインパルス応答が求められることになる。ここで，パワース ペクトル密度と，信号のフーリエ変換の関係を思い出すと， Sx(f ) = E[X∗(f ) X(f )] (11) Sxy(f ) = E[X∗(f ) Y (f )] (12) ということになり，つまりは観測した信号 x(t)，y(t) か，実際にはその離散時間表現 x(n), y(n) のフー リエ変換を求めて，上記の演算を行なうことでインパルス応答の推定が行なえることになる。フーリエ 変換を用いることが出来る，ということは演算時間短縮に関する大きなメリットである。
5.2.3 Maximum Length Sequence: MLS 法
直接相関法の計算で時間がかかるのはもちろん相関を求める部分である。ここに，白色性の擬似ラン ダム信号である M 系列信号 (maximum length sequence signal) を用いる場合，高速アダマール変換を 用いることで相互相関関数の計算が効率化され，きわめて高速にインパルス応答を求めることが可能に なる。この方法は建築音響関係の測定に多く用いられているが，音場の時変性の影響などを大きく受け ることが明らかになっている。 M系列信号とは長さ 1 L = 2N− 1 の周期を持つ 2 値（−1 と 1）系列である。ちなみに，mls.m と いうファイルに，M 系列ノイズを生成するアルゴリズムが記述されている。例えば， >> xm=mls(5);
と入力することで，25− 1 = 31 という長さの系列を生成する。 1) 適当な次数の M 系列信号を生成してみよう。 2) この信号のパワースペクトルや，自己相関関数などを調べてみよう。 5.2.4 おまけ：乱暴な方法 前述のように，入力信号と出力信号には以下の畳み込みの関係がある。 y(t) = h(t) ∗ x(t) この関係を，周波数領域に変換して考えてみる。この場合，畳み込みは単純な掛け算に置き換えること ができ，以下のような関係が求められる。 Y (f ) = H(f ) · X(f ) (13) 大文字は，それぞれの信号の周波数領域表現（例えば DFT した結果）とする。この式を見れば，例えば H(f ) = Y (f ) X(f ) (14) を求めればインパルス応答の周波数領域表現（伝達関数）が直接的に求められそうである。 しかし，例えば入力信号にノイズを用いた場合には，ランダムな信号によって生じた１回の結果から 確定的なインパルス応答を推測することはあまりにも乱暴であろう。このために，多数回の演算の平均 を求める操作が必要になりそうである。その他にも，いずれかの周波数成分に零点があれば，この演算 は成立しなくなる。乱暴な方法と呼ぶ所以である。
1) impulse data.matというファイルをロードしよう。ir, x, y, fs という 4 変数がロードさ れる。ir はサンプリング周波数を 1,000 Hz にダウンサンプルした，あるホールのインパルス 応答であり，長さは 2,000 サンプル（2 秒）である。入力信号 x は，最大値を 1 に基準化した 継続時間は 30 秒のガウシアンノイズであり，また出力信号 y は y = fftfilt(ir, x); によっ て得られた信号である。x, y という素材信号を使って，インパルス応答を推定して，正解であ る ir との比較を行おう。。 2) imp fft.mというスクリプトファイルは，上記の乱暴な方法を行なうものである。これを実行 してみよう。 今回の課題 1 imp fft.mを参考に，上記の直接相関法とクロススペクトル法でインパルス応答を測定する方法に ついて各自考案し，真のインパルス応答と，結果として得られたインパルス応答とを一枚の Fig. に 描いて，repo8-1.fig （直接相関法），repo8-2.fig （クロススペクトル法）というファイル名で提出す ること。なお，いずれの図でも subplot(2,1,1) に真のインパルス応答 (ir) を，subplot(2,1,2) に計算で得られたインパルス応答を示してください。また，どのような方策を行なったのか，メー ルの本文に必ず記述してください。
5.3 Swept-Sine 法
直接的にインパルス信号を用いずにインパルス応答を算出する方法のもうひとつの代表格がここで紹 介する Swept-Sine 法である。これまでに，インパルスを時間的に引き伸ばしたものとして TSP: Time Stretched Pulseと呼ばれることも多かったが，実際には用いる引き伸ばし時間が長く，もはやパルスと は考えにくい，ということで近年は Swept-Sine という呼び方に統一されつつある。 5.3.1 Swept-Sine 信号の生成 一般的によく知られている Swept-Sine 信号の生成方法は，以下の式に示すとおり，周波数領域で設 計して，それを逆フーリエ変換するというものである。 S(k) = ⎧ ⎪ ⎨ ⎪ ⎩ exp −j4mπk2 N2 , 0 ≤ k ≤N 2 S∗(N − k), N 2 < k < N (15) ここで，m = N/4, N = 2n, n は任意に設定する FFT の次数である。なお，この信号から生成される 音源信号と畳み込むことでインパルスを生じる，逆の特性を持つ信号は，以下の特性を逆フーリエ変換 することで求められる。 S−1(k) = ⎧ ⎪ ⎨ ⎪ ⎩ exp j4mπk2 N2 , 0 ≤ k ≤N 2 S−1∗(N − k), N 2 < k < N (16) tsp designという関数を用意している。これは上記の n を与えると音源信号 s とその逆特性を持った 信号 g を与えるものである。使用例は以下の通り。 >> [s, g] = tsp design(10); これで，210= 1024の長さの信号が得られる。 1) tsp designという関数を使って，Swept-Sine 信号の様子を観測しよう。 2) sと g は，逆の特性になっている（インパルスを作り出せる）のか，確認しよう。 このような信号を用いてインパルス応答を測定するわけであるが，その方法には大別して 2 種類ある。 円状畳み込みをを用いた測定方法と直線状畳み込みを用いる方法である。それぞれの特徴やメリット・ デメリットなどはウェブページ 7)に詳述されているが，ここではその内容を概説する。5.3.2 円状畳み込みの原理を使ったインパルス応答測定法 この方法の特徴は以下のようにまとめられる。 ◦ 用いる信号の長さ N と，取り込むべき信号長 L は等しい。(L = N) また，計測するインパルス 応答の長さ I との関係は，L > I であれば良い。 ◦ このため，あらかじめインパルス応答のおおよその長さを知っておく必要がある。 以下，具体的な測定手順を示す。 1. Fig. 5.3 (a)に示すインパルス応答を求めることとする。 2. このインパルス応答の長さは 2000 であり，これよりも長く 2 のべき乗である長さ 2048 の信号を 作る。具体的には Fig. 5.3 (b) のようになる。 3. 同期加算を行なうために，これを隙間無くならべる（Fig.5.3 (c))。今回は 3 回の加算を仮定して いるが，そのために (d) に示すように 4 回分の応答を取り込む。これは 3 回目の信号の響き部分 だけを取り込むためである。 4. 4つの区間の信号を加算することで，周期的応答を 3 回加算したことになる（Fig.5.3 (e))。 5. (f)のように音源信号と逆の特性を持つ信号を用意する。 6. 応答信号 (e) に逆特性を作用させる。具体的には，(e) の信号，(f) の信号共に N 点の FFT を行 い，結果を要素ごとに掛け算する。 7. さらにその結果を逆 FFT することで，(g) のようにインパルス応答を推定することができた。
(a) ߎߩࠗࡦࡄ࡞ࠬࠬࡐࡦࠬࠍ᷹ቯߒߚ (b) ߘߩߚߦ↪ᗧߒߚSwept-Sineାภ (c) หᦼട▚ߩߚߦାภࠍਗߴࠆ (d) 㖸႐ߩᔕ╵߇⇥ߺㄟ߹ࠇߚ⚿ᨐ (e) 㖸Ḯାภ㐳ߢಾࠅߒߡട▚ᐔဋߒߚ⚿ᨐ (f) ㅒ․ᕈࠍਈ߃ࠆାภ 200 400 600 800 1000 1200 1400 1600 1800 2000 500 1000 1500 2000 500 1000 1500 2000 2000 4000 6000 8000 2000 4000 6000 8000 500 1000 1500 2000 500 1000 1500 2000 (g) ផቯߐࠇߚࠗࡦࡄ࡞ࠬᔕ╵ Fig. 5.3: 円状畳み込みを用いたインパルス応答測定方法
5.3.3 直線状畳み込みの原理を使ったインパルス応答測定法 この方法の特徴は以下のようにまとめられる。 ◦ 用いる信号の長さ N と，計測するインパルス応答の長さ I との間に，特に大きな制限は無い。イ ンパルス応答よりも短い信号で測定可能である。 ◦ このため，短めの信号で実際に音場の響きを聞きながら測定できるというメリットがある。 以下，具体的な測定手順を示す。 1. Fig. 5.4 (a)に示すインパルス応答を求めることとする。 2. このインパルス応答の長さは 2000 であり，これよりも短く，2 のべき乗である長さ 1024 の信号 を作る。具体的には Fig. 5.4 (b) のようになる。 3. 同期加算を行なうために，これを適当な間隔をあけてならべる（Fig.5.4 (c))。今回は 3 回の加算 を仮定している。あける間隔は，インパルス応答長よりも長く設定する。この図の例では，イン パルス応答長+300 サンプルを仮定している。 4. 応答を取り込む際には，N + I − 1 よりも長い間隔で取り込むこと。 5. 3つの区間の信号を，先頭をそろえて加算する（Fig.5.4 (e))。 6. (f)のように音源信号と逆の特性を持つ信号を用意する。 7. 応答信号 (e) に逆特性を作用させる。直線状畳み込みを行なう関数としては，conv() あるいは fftfiltを用いると良い。 8. (g)のようにインパルス応答を推定することができた。 今回の課題
impulse dataという mat ファイルに格納されているインパルスレスポンスと，関数 tsp design を用いて，各自 Fig. 5.3, 5.4 を再現しなさい。見易さを工夫して，それぞれ repo9-1.fig, repo9-2.fig として提出してください。
(a) ߎߩࠗࡦࡄ࡞ࠬࠬࡐࡦࠬࠍ᷹ቯߒߚ (b) ߘߩߚߦ↪ᗧߒߚSwept-Sineାภ (c) หᦼട▚ߩߚߦାภࠍਗߴࠆ (d) 㖸႐ߩᔕ╵߇⇥ߺㄟ߹ࠇߚ⚿ᨐ (e) 㖸Ḯାภ㐳+ αߢಾࠅߒߡട▚ᐔဋߒߚ⚿ᨐ (f) ㅒ․ᕈࠍਈ߃ࠆାภ 500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 200 400 600 800 1000 2000 4000 6000 8000 2000 4000 6000 8000 1000 2000 3000 200 400 600 800 1000 (g) ផቯߐࠇߚࠗࡦࡄ࡞ࠬᔕ╵ ାภ㐳 ขㄟ㐳 Fig. 5.4: 直線状畳み込みを用いたインパルス応答測定方法
5.3.4 Swept-Sine 信号のバリエーション Swept-Sine信号は，その原理が比較的明快に与えられているため，バリエーションもいくつか提案さ れている。その中で，藤本によって提案された Pink TSP 信号10)を紹介しておく。これは低音域ほど スイープ音の掃引速度が遅く，振幅特性として Pink 特性を持ったものである。高音域での S/N を犠牲 にすることになるが，広帯域のバランスをとる場合には有効な手段として知られている。 S(k) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 k = 0 exp(jak log k)√ k , 0 < k ≤ N 2 S∗(N − k), N 2 < k < N (17) ここで，aN 2 log N 2 = 2mπ となるように設定し，また m = N/4, N = 2n, n は任意に設定する FFT の次数である。なお，この信号から生成される音源信号と畳み込むことでインパルスを生じる，逆の特 性を持つ信号は，以下の特性を逆フーリエ変換することで求められる。 S−1(k) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 k = 0 √
k exp (−jak log k) , 0 ≤ k ≤ N 2 S−1∗(N − k), N 2 < k < N (18) 1) tsp designという関数を参考に，上記の Pink-Swept-Sine 信号を生成しよう。 2) spectrogramなどを使って，周波数スウィープの様子を観測しよう。 3) sと g は，逆の特性になっている（インパルスを作り出せる）のか，確認しよう。 5.3.5 インパルス応答の活用例 最後に，得られたインパルス応答を加工して用いる例として，残響曲線を描く手法を修得しよう。室 内音響学においても述べたように，室内において得られたインパルス応答から，エネルギーの減衰曲線 を描くことができる。その関係式は以下の通りである。 p2 (t) = ∞ t h 2 (t) (19) ここに，p2(t) は室内における音圧自乗値，つまりエネルギー減衰過程のアンサンブル平均で，求めた い減衰曲線そのものである。この関係式は，インパルス応答 h(t) の自乗値を用いて，ある時間 t 以降 の積分値が得られれば，それが減衰曲線であることを示している。これは例えば，以下のコマンドで実 行することが出来る。
実際には，インパルス応答を所望の帯域フィルタに通して，その応答に上記のエネルギー積分を行なっ て，周波数ごとの減衰曲線を求めることが多い。この減衰曲線に直線を当てはめて，60 dB 減衰するま での時間を求めれば，残響時間が得られることになる。 1) real irというファイルに，サンプリング周波数 48000 Hz で測定した 2 秒間のインパルス応答 が ir という変数名で格納されている。これを load し，プロットしてみよう。 2) fir1などを使ってバンドパスフィルタを作り，これを通すことで 500 Hz 中心のオクターブバ ンドのインパルス応答を求めよう。 3) 上記の応答から，fliplr と cumsum を用いて減衰曲線を求めよう。 4) 残響時間を求めてみよう（ヒント： polyfit という関数で，最小自乗直線近似が行なえます）。
