JAIST Repository
https://dspace.jaist.ac.jp/

Title

変調伝達関数に基づく音声信号処理(1)パワーエンベロ
ープ逆フィルタ処理の原理とその応用について

Author(s)鵜木, 祐史Citation

信号処理, 12(5): 339‑348

Issue Date

2008‑09

Type

Journal Article

Text version

author

URL

http://hdl.handle.net/10119/7754

RightsCopyright (C) 2008 信号処理学会. 鵜木祐史, 信号処理, 12(5), 2008, 339‑348.

Description

Japan Advanced Institute of Science and Technology

変調伝達関数に基づく音声信号処理 (1)—パワーエンベロープ逆フィルタ処理の原理とその応用について—
Speech Signal Processing Based on the Concept of Modulation
Transfer Function (1)— Basis of Power Envelope Inverse Filtering and Its Applications —鵜木祐史
Masashi Unoki

1.

はじめに

私達は，日常，何不自由なく音声を介してコミュニ
ケーションをとっている。しかし，読者はこんな経験をしたことはないだろうか。例えば，お風呂場や教会など音が非常に響く環境（残響環境）や，人で賑わっている雑多な場所，交通量の多い場所といった非常に騒がしい環境（騒音環境）では，静寂な環境に比べて非常に音を聴き取り難く，いつもと同じように簡単に会話をできないと感じたことである。これは，身の回りの音場環境の影響により，音声が歪んだため，音声知覚に重要な情報が欠落したことによるものである1 。
このような音声コミュニケーションの難しさを評価する尺度として，音声明瞭度，単語・文章了解度が利用されている [2]。前者は無意味音節を発声したとき受聴者がその何％を正しく聞き取れたかを，後者は沢山の有意味単語を発声したとき受聴者が正しく聞き取れた
単語数の割合を示すものである。これらの尺度は，音声情報伝達を議論するときに，よく利用されるものであるが，同時に室内音饗学と関係して議論されるとき，
音声レベル，騒音レベル，残響時間等の物理量との関連を見出そうとする検討も古くから行われている [3]。
代表的なものとして，Houtgast と Steeneken によって提唱された変調伝達関数（Modulation Transfer Function: MTF）に基づく音声明瞭度予測理論 [4, 5, 6] がある。これは，音場内において，音声波形の時間的な包絡線情報（以後，エンベロープと呼ぶ）が残響や雑音によって変形することに着目し，100%振幅変調した正弦波を利用して MTF の減衰量から音声伝達指標（Speech
Transmission Index: STI）を予測するものである。この方法は，その後，簡易測定法である RASTI[7] として提案され，現在でも標準的な方法として利用されてい北陸先端科学技術大学院大学 情報科学研究科
923-1292 石川県能美市旭台 1-1
School of Information Science, Japan Advanced Institute of Science and Technology
1-1 Asahidai, Nomi, Ishikawa 923-1292, Japan
E-mail: unoki@jaist.ac.jp

る。STI/RASTI の方法は，理論的に明解であり，実用上多くの利点をもつため，講演会場など室内音響設計[2, 8] にも役立っている。しかしながら，この方法は決して万能であるわけではなく，(1) 音場の時間構造・空間構造を反映していないことや (2) 音源（音声）の物理特性を反映していないことから，音声明瞭度予測に対して適用限界があることが示唆されている [2]。
Houtgast と Steeneken が提唱した音声明瞭度予測
理論は，室内音響を拡散音場と仮定しているため [8]，上記のように，その適用限界があることは間違いない。
しかしながら，室内音響伝達系を入出力の強度情報の関係と残響・雑音に対する MTF を明解に関係づけた
点は，大きな業績であり，評価されることであろう。また，この考えは，他の音声信号処理で残されている諸問題を解決するために利用することもできる。例えば，室内の残響の影響を受けた音声を伝達系を測定せずに回復する方法 [9, 10, 11] や，残響環境下での音声の基
本周波数推定方法 [12] がある。最近では，室内の残響
時間をブラインド推定する方法 [13] や異なる二つの音場空間を考慮した音場再生法 [14] も提案されている。
本論文は，合計３回のシリーズで構成される。これ
らでは，著者が関係した研究分野（残響環境下の音声信号処理）を中心に，MTF を利用した音声信号処理を解説する。本稿では，まず，Houtgast と Steeneken が示した MTF の概念を解説するとともに，その概念に基づいたパワーエンベロープ逆フィルタ法を紹介する。

2.

変調伝達関数（MTF）

変調伝達関数（MTF）の概念は，音声明瞭度を評価するための室内音響の測定として，Houtgast とSteeneken[4, 5, 6] によって導入された。彼らの着眼点
は，室内音響において，残響や単一エコー，妨害雑音の変化によって得られる明瞭度のスコアと MTF から
1 聴覚の優れた能力の一つに音韻修復現象 [1] がある。これは雑音等で欠落した情報から情報を知覚的に補正・修復することで，正しく音声を知覚できるという現象であるが，言語や先験情報に依存しており，高次からの処理が重要と考えられている。同様に雑音・残響環
境下においてもこの能力が発揮されるが，この現象が説明されるような状況とは異なり，音声全体にわたって歪みが存在するため，誤った修復を行い，結果として異聴を引き起こすことが多いと思われる。

2

Journal of Signal Processing, Vol. , No. , 2008

Input intensity

System (Room acoustics)

Output intensity<x2(t)>

1/fm

1/fm

Intensity

Intensity(a)

Reverberation<h2(t)><y2 (t)>

Time

Time

1/2

2
y

1/2

y(t)={e (1+m(fm)sin2πfm(t-θ))} no(t)x(t)={ex2(1+sin2πfmt)} ni(t)(b)

Modulation indexm(fm)Modulation Transfer Function (MTF)
1

Modulation frequency, fm(c)Speech Transmission Index (STI)
Prediction of speech intelligibility

図1

Houtgast と Steeneken によって提案された音声明瞭度予測理論の概要図。(a) 室内音響における入出力信号の強度情報とその伝達関数。(b) 導出された変調伝達関数（MTF）m(fm )。(c) MTF から導出される音声伝達指標（STI）と音声明瞭度の予測
Fig. 1 Schematic illustration of prediction of speech intelligibility proposed by Houtgast and Steeneken: (a)Relationship between intensity of input and output signals and its transfer function, (b) modulation
trasnfer funcion (MTF), and (c) speech transmission index (STI) derived from MTF導出される指標が非常に高い相関関係にあることにある。彼らが想定している室内音響は拡散音場であるが，次のような二つの場合に分類している。(a) 残響に関係する伝達関数は，直接経路の他に，壁や床など様々な反射経路（時間差あり）の総和で構成されること，(b)妨害雑音に関係する伝達関数は，ある点（音源）から室内のある一点までの経路上における妨害量を表すものである。そのため，彼らは，室内を伝達系と見たときの入力，出力の強度変化に着目し，その中の変調度の変化を MTF と定義して，音声明瞭度に関係する指標 STI を予測する体系を作った。本稿では，MTF に特化するため，STI の導出法については触れないが，詳細については，[4, 5, 6] を参照されたい。
まず，彼らは，室内音響における入出力の強度変化を

Input:
Output:Ii (1 + cos(2πfm t))(1)

Io (1 + m(fm ) cos(2πfm (t − θ))) (2)と定義した。ここで，Ii と Io はそれぞれ入出力の強度，fm は変調周波数，θ は位相を表す。入力の強度を
パワーエンベロープとみなしたときに，ちょうど変調度m(fm ) = 1（100%振幅変調）の余弦波をもつものとする。出力は同様の形態を取るが，余弦波の振幅が雑音や残響の影響を受け，ちょうど m(fm ) 倍に変化し，これが振幅変調で見たときの変調度に相当することから，
MTF と呼ばれる。

次に m(fm ) の導出を説明するために，入力信号 x(t)と出力信号 y(t) のパワーエンベロープを

e2x (t) = e2x (1 + cos(2πfm t))e2y (t) = e2y (1 + m(fm ) cos(2πfm t))

(3)
(4)

と定義し（便宜上，θ = 0），次式のように一般化する。

x(t) =
y(t) =ex (t)nx (t)ey (t)ny (t)

(5)
(6)

但し，nx (t) ならびに ny (t) は，白色ガウス雑音の特性を有するランダム変数である。このランダム雑音には次式の特性がある（x と y の添字は省略する）。

< n(t)n(τ ) >= δ(t − τ )(7)

ここで，< · > は集合平均を表す [15]。この特性を活用
し，x(t) と y(t) の 2 乗集合平均を求めると，< x2 (t) > = e2x (t)< y2 (t) > = e2y (t)を得る（導出過程については付録１を参照）。

(8)
(9)

3

2.1

残響環境 h(t) における MTF

1
TR = 0.1 s

0.90

両辺を 2 乗した集合平均 < · > を求めると，次式を得る（詳細な導出過程については付録２を参照）。
]
[
e2
β
(11)
< y2 (t) >= e2y (t) = x 1 + cos(2πfm t)
α
α
∫∞
∫∞
但し，α = 0 h2 (τ )dτ ，β = 0 h2 (τ ) exp(−jωm τ )dτである。ここで，式 (4) との対応を見れば，m(fm ) =β/α がわかり，次式に示す複素表現の MTF を得る。
∫∞ 2
h (t) exp(−jωt)dt
m(ω) = 0 ∫ ∞ 2(12)h (t)dt0
これは h2 (t) に対する一種の Fourier 変換を表している。
次に，式 (12) のインパルス応答 h(t) を次式のように定義する。これは，室内音響特性の統計的近似として良く知られた Schroeder のインパルス応答 [16] である。
()
6.9t
h(t) = a exp −nh (t) = eh (t)nh (t), t ≥ 0
TR(13)
但し，a は振幅項，TR は T60（h(t) のパワーがそのピーク
からちょうど 60 dB 減衰するとるまでに要する時間）で定義される残響時間である [8, 16, 17]。このとき，MTF
は，式 (12) に式 (13) を代入することで得られる。
[()2 ]−1/2TRm(fm ) = |m(2πfm )| = 1 + 2πfm
13.8
(14)
図 2 は，変調周波数 fm（時間的なパワーエンベロー
プにおける主要な周波数）の関数とした式 (14) の MTF
m(fm ) の特性図を示す。これらの曲線は，TR = 0.1, 0.3,
0.5, 1.0, 2.0 s と ω = 2πfm としたときのものである。
m(fm ) は，100%振幅変調したときの fm に対する変調度と解釈することもできる。これらの曲線は，TR と fmに依存しており，その概形からわかるように，残響環境の MTF は一種の低域通過フィルタ特性を示している。

2.2

雑音環境として，定常な加法性雑音 w(t) を考える。

(15)

ここで，図 1 における入出力関係は

y(t) = x(t) + w(t)

0.7
0.6

TR = 0.3 s

0.5
0.402
0.4

TR = 0.5 s

0.3
0.2

T =1s

0.1

TR = 2 s

0

R

0

5

(16)

である。前節にならい，両辺の 2 乗集合平均を取ると< y2 (t) > = < x2 (t) > + < w2 (t) >e2y (t) = e2x (t) + e2w (t)
(
)
)
(
e2x
2
2
ex + ew1+=cos(2πfm t)
e2x + e2w

1015
Modulation frequency, fm (Hz)20

図2

様々な残響時間（TR = 0.1, 0.3, 0.5, 1.0, 2.0 s）
における変調伝達関数 m(fm ) の特性
Fig. 2 Theoretical curves representing the MTF,
m(fm ), for various conditions with TR = 0.1,
0.3, 0.5, 1.0, and 2.0 sを得る。但し，e2w (t) は定常と仮定され，e2w (t) = e2w =∫1 T 2
T 0 ew (τ )dτ とする。ここで T は信号長である。このとき，式 (3) と式 (4) の関係を眺めると，余弦波にかかる m(fm ) は fm に依存せず，次式のように一定な減衰を表すことがわかる。

m(fm ) =e2x
e2x + e2w=

1
1 + 10−(S/N )/10

(17)

但し，S/N は 10 log10 e2x /e2w (dB) である。例えば，e2x (t)と e2w (t) がほぼ同じ強度であるとき（S/N = 0 dB），
MTF は m(fm ) = 0.5 となる。

2.3 雑音残響環境（w(t) と h(t)）における MTF
残響と雑音の両方を考慮した場合の変調伝達関数を考える。前節にならい，図 1 の入出力関係を考えると∫ ∞h(τ )x(t − τ )dτ + w(t)y(t) =0e2y (t) =雑音環境 w(t) における MTF

w(t) = ew (t)nw (t)0.8

MTF, m(fm)

残響環境における出力信号 y(t) は入力信号 x(t) と室内インパルス応答 h(t) の畳み込み演算から得られる。
∫ ∞y(t) =h(τ )x(t − τ )dτ
(10)

[
]
β
e2x
1 + cos(2πfm t) + e2w (t)ααを得る。式 (3) と式 (4) の関係を眺め，前節の導出過程を踏まえると，余弦波にかかる m(fm ) は，次式のように残響・雑音の影響をそれぞれ受けることがわかる。

1e2x×m(fm ) = √()
e2x + e2w
TR 2
1 + 2πfm 13.8

(18)以上が，Houtgast と Steeneken によって導入された
MTF の概念である。MTF に着目すれば，残響は変調度に関する低域通過フィルタの効果，雑音は変調度のゲ
インの減衰効果をもたらすことがわかる。本論文では，残響環境下での音声信号処理に焦点を当てるため，以後，残響環境での MTF の概説に特化することにする。

4

Journal of Signal Processing, Vol. , No. , 20083.

パワーエンベロープ逆フィルタ処理

3.1 処理方法
室内音響における入出力関係式（x(t)，h(t)，y(t)）と MTF m(fm ) の対応を再度眺める。２節ではパワー
エンベロープ（e2x (t) と e2y (t)）を単一周波数成分からなるものと仮定したが，これを限定しない形で式 (10) を展開，整理すると，e2x (t)，e2h (t)，e2y (t) の間には，

⟨ 2 ⟩
y (t)

⟨{∫
∞
=
∫ ∞
=
−∞

−∞

}2 ⟩
x(τ )h(t − τ )dτ

e2x (τ )e2h (t − τ )dτ = e2y (t) (19)の関係があることがわかる（詳細な導出過程は付録３を参照）。同様に，式 (14) と式 (19) の MTF の関係を眺めると m(fm ) は e2h (t) に関するパワーを正規化した一種の Fourier 変換であることもわかる。
以上を踏まえると，MTF に基づく定式化においては，信号だけでなくパワーエンベロープに関しても線形性が成り立っていることがわかる。このことから，仮に残響の効果（e2h (t)）を事前に知ることができれば，この特性が最小位相特性を有するときに限り，e2y (t) から原信号のパワーエンベロープ e2x (t) を求めることができる。これは逆畳み込みの関係であることから，パワー
エンベロープ逆フィルタ処理と呼ばれる [9]。
次に，これらを計算機で取り扱うために，標本化定理に基づき連続時間信号を離散時間信号に変換する。ここでは，サンプリング周波数 fs = 20 kHz，サンプル
点 n とし，e2x [n]，e2h [n]，e2y [n]，x[n]，h[n]，y[n] と表記する。これら三つのパワーエンベロープの z 変換をそれぞれ，Ex (z)，Eh (z)，Ey (z) とする。ここで，インパルス応答のパワーエンベロープの z 変換を考える。
式 (13) のパワーエンベロープ e2h (t) は単純な指数減衰関数であり，t < 0 では e2h (t) = 0 と仮定していることから，最小位相特性を有することがわかる。この特性を踏まえ，e2h [n] の z 変換 Eh (z) を求めると

Eh (z) =(a2)
1 − exp − T13.8
z −1
·f
R s

(20)

を得る（導出過程については付録４を参照）2 。式 (19)の畳み込みの関係から，e2x (t) の変調スペクトル Ex (z)は e2y (t) の変調スペクトル Ey (z) と MTF の逆数（逆フィルタ，IMTF）1/Eh (z) の積で求められる。
{
(
)
}
13.8
Ey (z)−1
1 − exp −z(21)
Ex (z) =a2
TR · fs
この逆フィルタは，1 次の IIR フィルタであり，e2y [n]に適用することで，最終的に e2x [n] を得る [9]。
図 3 は，MTF に基づくパワーエンベロープ逆フィルタ処理を説明するものである。図 3 (a) は，e2x (t) =2 パワーエンベロープ畳み込みの他，エンベロープを 2 乗しない形での展開も議論されている [18, 19]。

0.5(1 + sin(2πfm t)) を示す。ここで，fm は 10 Hz であり，変調度は１である。図 3(b) は，式 (5) を利用して得られた原信号 x(t) = ex (t) · nx (n) を示す。図 3 (c) は，
TR = 0.5 s としたときの式 (13) の e2h (t) を示す。図 3(d) は，残響インパルス応答 h(t) = eh (t) · nh (t) を示す。
図 3(e) と図 3(f) は，式 (19) のように e2x (t) と e2h (t) の畳み込みから得られた e2y (t) と，x(t) と h(t) の畳み込みから得られた y(t) を示す。図の左側（(a), (c), (e)）はパ
ワーエンベロープを，図の右側（(b), (d), (f)）は各信号波形を示す。この図では，変調度が 1.0（図 3(a)）から0.404 (図 3(a) の点線と図 3(e) の点線の間の偏差の最大値) に減衰していることがわかる。MTF は，fm と TRの関数とした変調度 [6] を示すため，式 (14) に fm = 10
Hz と TR = 0.5 s を代入することで，m(fm ) = 0.402 を得る。図 3 (g) の実線は，TR = 0.5 として，式 (21) を利用した際の，e2y (t)（図 3 (e)）から求めた ê2x (t) を示す。この結果から，パワーエンベロープ逆フィルタが，パワーエンベロープの大きさと形状に関して，正確に回復可能であることを示している。
3.2 応用上の問題点
上記で概説したパワーエンベロープ逆フィルタの原理は，広林ら [9] によって最初に提案されたが，次に述べるような原理上の問題点が多々残されていた。
図 4 は，パワーエンベロープ逆フィルタ処理のブロッ
クダイアグラムを示す。このモデルの概念において，観測された残響信号 y(t) から完全に正確にパワーエンベ
ロープ e2y (t) を抽出でき，かつ室内インパルス応答のパラメータ（振幅項 a と残響時間 TR ）を事前に知っているという条件であれば，図 4 の方法は観測された残響信号 y(t) から原信号のパワーエンベロープを正確に回復することができる。しかし，e2y (t) が y(t) から正確に抽出されなかった場合，この方法は式 (21) の逆フィル
タ法を利用して y(t) から e2x (t) を正確に回復できない。
更に，e2y (t) が y(t) から完全に正確に抽出されたとしても，二つのパラメータ（a と TR ）の値が不正確であれ
ば，この方法は正しく ê2x (t) を回復できない。
例えば，図 3 (g) の実線とは対照的に，逆フィルタが不適切なパラメータ値（例えば，TR = 0.3 s や TR = 1.0s）で適用された場合，図 3 (g) の破線，一点鎖線に示すように，e2x (t) は正確に回復されない。ここで，前者の場合，“過小回復” であり，後者の場合，“過剰回復”となる。特に，TR = 1.0 s（図 3 (g) の破線）のとき，
ê2x (t) は e2y (t) から過剰に処理され，これは負の値をもつことになる（物理的にパワーが負になることは有り得ない）。ここで，êx (t)2 の改善度は，TR = 0.3 s（3(g) の点線）のときよりも高くなり，TR = 0.5 s での回復より低くなる。これは，式 (19) が積分処理であるのに対し，式 (21) が微分処理であることに起因している。
以上のことから，パワーエンベロープ逆フィルタ処理を利用するには，次の二つの問題点を解決しなければならない。(1) パワーエンベロープをどのようにして正確に抽出するか，(2) 室内インパルス応答のパラメー
タ（振幅項 a と残響時間 TR ）をどのように決定するか。
これらの解決法について，次節で概説する。

5

(a)

(b)
4
2
x(t)

2

ex (t)

1
0.5

0
−2

0
−3

(c)

x 10

(d)
0.2
TR =0.5 (s)

h(t)

2

eh(t)

4
2

0
−0.2

0
(e)

(f)
4
2
y(t)

2

ey (t)

1
0.5

0
−2

0

0

(g)

0.2

2

^ ex (t)0.4

0.6Time (s)21TR=0.3 (s)T =0.5 (s)R
T =1.0 (s)

0

R

0

0.2

0.4

0.6

Time (s)図 3

MTF の概念に基づいた伝達系のパワーエンベロープ間の関係。(b) 入力信号 x(t) と (a) そのパワーエンベ
ロープ e2x (t)，(d) 室内音響インパルス応答 h(t) と (c) そのパワーエンベロープ e2h (t)，(e) e2x (t) ∗ e2h (t) から得られたパワーエンベロープ e2y (t)，(f) x(t) ∗ h(t) から得られた出力信号 y(t)，(g) 回復されたパワーエンベ
ロープ ê2x (t)Fig. 3 Example of the relationship between the power envelopes of a system based on the MTF concept: (a)power envelope e2x (t) of (b) original signal x(t), (c) power envelope e2h (t) of (d) impulse response h(t),(e) power envelope ey (t)2 derived from e2x (t) ∗ e2h (t), (f) reverberant signal y(t) derived from x(t) ∗ h(t),and (g) restored power envelope ê2x (t)

3.3

改良法

3.3.1 パワーエンベロープ抽出法
パワーエンベロープ抽出法に関しては，信号伝送処理で扱われるような振幅変調信号に対する復調処理を考えることで，適切な解法を得ることができる。しかし，MTF に基づくパワーエンベロープ逆フィルタ処理
では，キャリアが正弦波ではなく白色雑音であるため，一般的な復調方式（例えば，半波整流法 [20] など）では対応できない。これに対し，著者の研究グループに
よって二つの抽出法が提案された [10]。一つは，式 (9)の集合平均を利用する方法である3 。もう一つの方法は，Hilbert 変換の関係 [21] を利用したものである。両手法
とも同程度の抽出精度を有するが，計算が簡単で計算量が低いことから，本稿では後者を紹介する。
3 式 (9) では，キャリアが白色ガウス過程の確率変数としている
ため，実際の観測信号 y(t) から集合平均を求めることができない。
白色ガウス過程に基づいて生起された白色雑音は互いに無相関であることから，相異なる白色雑音を掛け合わせたものも白色雑音であると考える。提案法は，この仮定に基づき，有限個の白色雑音を観測信号に掛け合わせて得た擬似的な観測集合 y(t) の集合平均を求めることで，e2y (t) を得るものである。

Hilbert 変換器は，90 度位相器と呼ばれ，偶関数・奇関数の変換や解析信号からの瞬時振幅の算出によく利用される [21]。観測信号 y(t) のキャリアが偶関数あるいは奇関数で構成されていると仮定すれば，Hilbert
変換を利用して瞬時振幅を求め，これに低域通過フィ
ルタ処理（LPF）をすることでエンベロープを抽出することができる。この方法を次に示す。
[]
2ê2y (t) = LPF |y(t) + j · Hilbert(y(t))|(22)
ここで，LPF のカットオフ周波数は 20 Hz とした。これは音声の主要な変調周波数が 1〜16 Hz にあるという報告に基づいたためである [10].

3.3.2 インパルス応答のパラメータ推定
パワーエンベロープ逆フィルタ法では，残響時間 TRと振幅項 a が逆フィルタを構成する主要素となるため，残響回復の前にこれらの正確な値を知らなければならない。著者の研究グループでは，次のような考えに沿って，パラメータ推定法を提案した [10]。
まず，図 2 に示すように，パワーエンベロープの変調周波数 fm が単一の場合を考える。先に述べたように，

6

Journal of Signal Processing, Vol. , No. , 2008

Reverberant signal

Power envelope
extractiony(t)e 2y (t)A

Recovered power-envelope

Power-envelope
inverse-filteringC^2
e x (t)B
Parameter (a, TR)estimation

図 4 パワーエンベロープ逆フィルタ処理のブロックダイアグラム
Fig. 4 Block diagram of the power envelope inverse ﬁltering method

fm と m(fm ) の関係は，残響による e2x (t) の変調度の減衰量を表しているのだから，1/m(fm ) 倍（IMTF）だ
け e2y (t) の変調度を補正することが，残響回復となる。
そのため，完全な回復量を示す TR を探し出すことで，MTF を形成したときの TR を知ることができる。この算出は，fm が単一周波数のときに容易であるが，複数周波数のときは難しい。そこで，IMTF の算出を変調
周波数領域ではなく，時間領域で行うことを考える。図2(g) に示したように，パワーエンベロープの過小・過剰処理の結果から，残響を付与したときの残響時間と同じ TR の値を利用してパワーエンベロープを回復すると，ê2x (t) が e2x (t) と一致することがわかる。
パワーエンベロープは少なくとも一つの 0 の値を取るディップないし無音区間をもつため，ここでは，e2x (t)の変調度が 1 であると仮定する。上記の発見に基づき，オリジナルと回復されたパワーエンベロープの間の一致条件が，残響によって減少した変調度を回復することであると定義する。この条件は，ちょうど TR を推定する条件となっており，次式のように，回復されたパワーエンベロープ ê2x (t) の最大のディップが 0 であると
ころ，あるいは ê2x (t) の負の面積がちょうど 0 であるところを検出することで推定することができる [10]。
(
)
∫ T
¯
)¯
( 2
¯
¯
T̂R = max
arg minmin êx,T (t), 0 dt
TR,min ≤TR ≤TR,max

R

0

(23)

但し，T は y(t) の信号長，ê2x,TR (t) は TR を変数として得られたパワーエンベロープである。“max(arg min{·})”
は，これらの候補 ê2x,TR (t) の中から，負の面積がちょうど 0 になるときの TR の最大値を求めることを意味する。そのため，この式は，回復されたパワーエンベ
ロープがちょうど負の値をもつ前のところで制約を受けていることを意味する。ここで，TR,min と TR,max はそれぞれ TR の取り得る範囲の下限値と上限値を表す。
例えば，図 2(g) に TR = 0.3, 0.5, 1.0 s のときに回復した ê2x,TR (t) の候補を示す。ここで，TR,min = 0.0
s，TR,max = 3.0 s とした。 TR = 0.3 s のときの候補
ê2x,TR (t) はちょうど過小回復の状態にあり，TR = 1.0
s のときの候補 ê2x,TR (t) は過剰回復の状態にある。式(23) を利用して TR を推定すると，最適値 T̂R = 0.5 s
が三つの候補の中で最適であることがわかる。
次に振幅項 a は，室内インパルス応答の増幅度に関係する。式 (13) と式 (21) では， 同じ値を利用したが，実際，室内における残響効果は信号の増幅よりも反射による遅延の重ね合わせの効果が強い。そこで，e2h (t)の伝達特性と MTF の整合をとるために次式のようにパワーによる正規化を行うことにした。
√ ∫
(
)
T
13.8t
exp −
dt
(24)
â = 1/
T̂R
0
実環境における応用では，この値を推定しなければい
けないかもしれないが，室内インパルス応答のパワーエ
ンベロープから最適化手法により求めることもできる。

4.

評価シミュレーション

ここでは，著者の研究グループによって改良された
パワーエンベロープ逆フィルタ法 [10] を評価する。そこで，次に示す三つのパワーエンベロープと白色雑音の積で構成される原信号 x(t) を利用する。

1. 正弦波: e2x (t) = 1 − cos(2πF t)1 ∑sin(2πkF0 t + θk )K
K

2. 調波複合音: e2x (t) = 1 +k=1

3. 帯域制限雑音: e2x (t) = LPF[nw (t)]但し，F = 10 Hz, F0 = 1 Hz, K = 20, θk はランダム位
相，LPF[·] のカットオフ周波数は 20 Hz とした。ここでは，パワーエンベロープ抽出精度とパワーエンベロー
プの回復精度の両方を同時に評価する。評価シミュレー
ションでは，三つのパワーエンベロープに対し，100 個のキャリア（白色雑音）を乗じて得た音源信号 x(t) および，５種類の残響時間（TR = 0.1, 0.3, 0.5, 1.0, 2.0 s）
に対して，それぞれ 100 種類の室内インパルス応答 h(t)を畳み込んで得られた合計 1, 500（= 3 × 5 × 100）個の残響信号 y(t) を用意する。評価尺度は，パワーエン
ベロープに対する (1) 相関値と (2) SNR（S をオリジナ
ル，N をオリジナルと推定したものとの差）とした。
Corr(e2x , ê2x )=∫T ( 2√{∫T)(ex (t) − e2x (t)0)

ê2x (t) − ê2x (t) dt} {∫(e2x (t) − e2x (t))2 dt
0T(ê2x (t) − ê2x (t))2 dt
0

}

(25)
SNR(e2x , ê2x )=∫T

10 log10 ∫ T
00(e2x (t))2 dt(e2x (t) − ê2x (t))2 dt(dB)

但し，e2x (t) は e2x (t) の平均値を表す。

(26)

7

Correlation

1

0.95

0.9
(a)

Sinusoid
Harmonics
Random(b)

Sinusoid
Harmonics
RandomSNR (dB)

0.85
30

20

10

0

0

0.5

1
Reverberation time, TR (s)1.52

図 5 パワーエンベロープの抽出精度：(a) 相関値，(b)SNR
Fig. 5 Extraction accuracy of the power envelope:(a) correlation, (b) SNR

4.1

パワーエンベロープ抽出精度の評価

はじめに，パワーエンベロープ抽出法の精度について評価する（図 4 の A 点での評価）。図 5 に，x(t) な
いし，y(t) から抽出された三つのパワーエンベロープの抽出精度を示す。これは式 (22) を利用して得られた結果であり，図 5(a) は相関値を，図 5(b) は SNR を示す。e2x (t) の抽出精度については，オリジナルとの差を，
e2y (t) の抽出精度については，式 (19) から得られたパ
ワーエンベロープとの差を評価している。TR = 0 s のときは，e2x (t) の評価に相当する。この結果から，残響時間にほとんど影響を受けずに，高い精度でパワーエ
ンベロープを y(t) から抽出できることが確認された。

4.2

残響時間 TR の推定精度の評価

次に，式 (23) を利用した残響時間の推定精度を評価する（図 4 の B 点での評価）。図 6(a) は，各エンベロー
プに対して推定された残響時間 T̂R を示す。各点は，各100 個の観測信号から残響時間を推定したときの平均値を示し，エラーバーはそのときの標準偏差を示す。図中の点線は，オリジナルの残響時間 TR を示す。残響時
間が 0.5 s 付近までは，推定値がオリジナルの TR とほ
ぼ等しくなっているが，それ以降は徐々に差が生じる傾向にある。この影響については，次節で説明する。

4.3

パワーエンベロープ回復精度の評価

最後に，図 5 と図 6(a) の結果を踏まえ，式 (21) から求められた ê2x (t) の回復精度を評価する（図 4 の C
点での評価）。ここでは，回復精度として，相関値の改善度 Corr(e2x , ê2x ) − Corr(e2x , e2y ) と，SNR の改善度
SNR(e2x , ê2x ) − SNR(e2x , e2y ) を評価した。図 3(a) と図3(e) に示すように，各エンベロープの変調度の関係（MTF) は，TR が増加するにつれ減衰する。同様に，

Corr(e2x , e2y ) と SNR(e2x , e2y ) も，TR の増加とともに減衰する。そのため，ê2x (t) が正しく回復されたのであれ
ば，相関値と SNR の改善度は，両方とも正の値を示すはずである。もしいずれかの値が負の値を示すのであれば，ê2x (t) が適切に回復されたとはいえない。
図 6 に，パワーエンベロープ逆フィルタ処理による回復精度を示す。図 6(b), (d), (f) は，それぞれ，三つのパワーエンベロープに対する相関値の改善度を，図6(c), (e), (g) は SNR の改善度を示す。図中の実線は，式 (23) を利用して TR を推定して ê2x (t) を求めたときの結果を，破線は TR を既知として ê2x (t) を求めたときの結果を示す。各点とエラーバーは，それぞれ結果の平均値と標準偏差を示す。図 6(b)-(g) の改善度は，すべての条件で正の値を示していることから，パワーエン
ベロープ逆フィルタ法が，残響によって歪んだパワー
エンベロープを非常に効果的に回復しているといえる。
図中の実線と破線の改善度の違いについて検討する。
図 6 (b)-(c) の TR が 0.3〜1.0 s のときの結果と図 6(d)(g) の TR が 0.5〜2.0 s のときの結果を除き，二つの改善度には大きな違いが見られない。逆に，ここであげた違いが見られる場所は，図 6(a) に示す TR の推定精度が低下しているところと関係している。残響時間が比較的短い場合（TR < 0.5），推定値 T̂R は非常に点線に近い値となっており，相関値と SNR の改善度も非常によい値になっている。しかし，TR が 0.5 を超えたあたりから，その推定精度は低下し，SNR の改善度も正の値を示しているが若干低下する傾向にある。相関値の改善度は逆に増加するか同じ値を保持していることから，TR 推定の違いが回復されたパワーエンベロープの形状には大きな影響を与えていないということがわ
かる。SNR の改善度が低下しているのは，TR の過小推定が振幅項 â に影響を与えているためと考えられる。
比較的長い残響時間（TR > 0.5）では，TR を既知とした場合よりも，TR を推定したほうがよい改善度を示している。これは，図 6(a) の TR 推定の結果が悪いのではなく，図 6(b)-(g) に示す回復精度のレベルで最良値を取るような制約を式 (23) が表しているといえる。
上記について，事例を通して検討する。図 7 に，TR =1.0 s のときの正弦波ならびに帯域制限雑音のエンベロー
プに対する回復結果の一例を示す。図 7(a) と図 7(c) の
e2y (t) は，y(t) から抽出されたものである（破線）。図中の１点鎖線と実線は，それぞれ，e2x (t) と式 (19) を利用して得られた e2y (t) を示す。図 7(b) と図 7(d) で示される負の面積は，パワーエンベロープのディップの最大値付近で検出される（図 7(b) の場合，0.1 s 付近；図7(d) の場合，0.2 s 付近）。ここで，T̂R は図 7(b) のとき
0.73 s，図 7(d) のとき 0.78 s であった。図 7(a) の e2x (t)と e2y (t) の間の相関値と SNR は 0.15 と 4.44 dB であったのに対し，図 7(b) の e2x (t) と ê2x (t) の間の相関値とSNR はそれぞれ 0.98 と 15.09 dB であった。これらの結果から，相関値と SNR の改善度は，それぞれ，0.83 と10.65 dB であった。対照的に，TR を既知としてパワー
エンベロープを回復したとき，その改善度は，相関値で0.82（= 0.97 − 0.15），SNR で 7.59（= 12.03 − 4.44）
dB であった。そのため，TR を既知として利用するよりも推定した TR を利用したほうが，適切な改善度が得られることがわかる。このときの効果は，相関値で8Journal of Signal Processing, Vol. , No. , 2008

Improved correlation

1
2
(a)
1.8

1.4
1.2

0.8
0.6
0.4

Sinusoid
Harmonic
Random

0.2
0

0

0.5

1
1.5
Reverberation time, TR (s)Estimated TR
Known TR

0.2

Estimated TR
Known TR

(e)
10
5
0
0

0.5

1
1.5
Reverberation time, TR (s)

2

1
(b)

0.8
0.6
0.4
Estimated TR
Known TR

0.2

Improved correlation

Improved correlation

0.4

2

1

(f)
0.8
0.6
0.4
Estimated TR
Known TR0.2
015

Estimated TR
Known TR

(c)
10
5
0
0

0.5

1
1.5
Reverberation time, TR (s)2

Improved SNR (dB)015
Improved SNR (dB)

0.6

0
15

1

Improved SNR (dB)

Estimated TR (s)

1.6

(d)
0.8

Estimated TR
Known TR

(g)
10
5
0
0

0.5

1
1.5
Reverberation time, TR (s)2

図6

評価シミュレーションの結果。(a) 残響時間の推定値（点線は理想化された残響時間）。正弦波（(b) と (c)），
調波複合音（(d) と (e)），帯域制限雑音（(f) と (g)）のパワーエンベロープに対する回復精度の比較（相関値の改善度と SNR の改善度）
Fig. 6 Results of simulations for evaluation: (a)Estimated reverberation time. The dotted line shows the idealized
reverberation time. Comparisons with the envelope restoration accuracy for sinusoidal power envelope ((b)
improved correlation and (c) improved SNR), a harmonic power envelope ((d) improved correlation and(e) improved SNR), and a band-limited noise power envelope: (f) improved correlation and (g) improved
SNR

0.01，SNR で 3.06 dB であった。
図 7(d) に，帯域制限雑音をパワーエンベロープとした場合に対する同様の評価を示す。ここでは，TR を既知とした場合，相関値で 0.57 （= 0.97 − 0.40），SNR
で 6.44 dB（= 12.8 − 6.36） dB であったのに対し，TRを推定してエンベロープを回復した場合，相関値で 0.58（= 0.98 − 0.40），SNR で 8.04 dB（= 14.4 − 6.36）であった。この結果， TR を推定したことによる効果は，相関値で 0.01，SNR で 1.60 dB であった。
5.

おわりに

本稿では，変調伝達関数に基づく音声信号処理（全３回シリーズ）の第１稿として，Hougast と Steeneken
によって導入された変調伝達関数の考えと，変調伝達関数に基づいたパワーエンベロープ逆フィルタ処理の原理を概説した。また，この原理を実際の信号処理に適用するときの問題点と改善法を紹介した。次回（第2 稿）では，本稿で紹介したパワーエンベロープ逆フィルタ処理を，残響音声回復に応用した手法を概説する。
謝辞

本論文で紹介する研究は，科学研究費補助金若手研
究 B（No.14780267），若手研究 A（No. 18680017），
萌芽研究（No. 17650048），科学技術振興調整費（若手研究支援プログラム），矢崎科学技術振興記念財団（特定研究助成）ならびに総務省 戦略的情報通信研究開発推進制度（課題番号 071705001）の援助を受けて行われた。研究協力者である，本学 赤木正人教授，Lu
Xugang 博士，本学修了生の古川正和君，酒田恵吾君，
戸井真智君，柴野洋平君，細呂木谷敏弘君，平松壮太
君，本学在学生の山崎 悠君に心より感謝する。

92Power envelope(a)Original e (t)x21

Original e (t)y2Extracted e (t)y0.5Restored power envelope

0

2

(b)

Original e (t)x1

Estimated T

R

Known T

R

0.5

0

2

Power envelope(c)Original e (t)x21

Original e (t)y2Extracted e (t)y0.5Restored power envelope

0

2

(d)

Original e (t)x1

Estimated T

R

Known T

R

0.5

0
0

0.1

0.2

0.3

0.4
Time (s)

0.5

0.6

0.7

0.8

図 7 結果例：(a) 正弦波信号のパワーエンベロープと(b) TR = 1.0 s のときに回復されたパワーエン
ベロープ，(c) 帯域制限雑音のパワーエンベロー
プと (d) TR = 1.0 s のときに回復されたパワー
エンベロープFig. 7 Sample results: (a) extracted power envelopeand (b) restored power envelope for sinusoidal
with TR = 1.0 s; (c) extracted power envelopeand (d) restored power envelope for a bandlimited noise with TR = 1.0 s

参考文献[1] A.S. Bregman: Auditory Scene Analysis: The Perceptual
Organization of Sound. MIT Press, Cambridge, MA, 1990.
[2] 戸井田 義徳: 小特集–音声の明瞭度と認識率–, 空間内における音声情報伝達, 音響誌, Vol. 51, No. 4, pp. 312–316, 1995.
[3] N. R. French and J. C. Steinberg: Factors governing theintelligibility of speech sounds, J. Acoust. Soc. Am., Vol.
19, pp. 90–119, 1947.
[4] T. Houtgast and H. J. M. Steeneken: The modulationtransfer function in room acoustics as a predictor of speech
intelligibility, Acustica, Vol. 28, pp. 66–73, 1973.
[5] T. Houtgast, H. J. M. Steeneken, and R. Plomp: Predicting speech intelligibility in rooms from the Modulation

Transfer Function. I. General Room Acoustics, Acustica,
Vol. 46, pp. 60–72, 1980.
[6] T. Houtgast and H. J. M. Steeneken: A review of the
MTF concept in room acoustics and its use for estimating
speech intelligibility in auditoria,” J. Acoust. Soc. Am.,
Vol. 77, pp. 1069–1077, 1985.
[7] 小椋靖夫，浜田晴夫，三浦種敏: 音場における音声伝達品質のための MTF と STI, 音響誌, Vol. 40, pp. 181–191, 1984.
[8] H. Kuttruﬀ: Room Acoustics, Taylor & Francis, fourth
edition, London, 2000.
[9] 広林茂樹, 野村博昭, 小池恒彦, 東山 三樹夫: パワーエンベロー
プ伝達関数の逆フィルタ処理による残響音声の回復, 信学論 A,
Vol. J81-A, No. 10, pp. 1323–1330, 1998.
[10] M. Unoki, M. Furukawa, K. Sakata and M. Akagi: An
improved method based on the MTF concept for restoring
the power envelope from a reverberant signal, Acoustical
Science and Technology, Vol. 25, No. 4, pp. 232–242, 2004.
[11] M. Unoki, K. Sakata, M. Furukawa, and M. Akagi: A
speech dereverberation method based on the MTF concept in power envelope, Acoustical Science and Technology, Vol. 25, No. 4, pp. 243–254, 2004.
[12] M. Unoki and T. Hosorogiya: Estimation of fundamental
frequency of reverberant speech by utilizing complex cepstrum analysis, J. Signal Processing, Vol. 12, No. 1, pp.
31–44, Jan. 2008.
[13] S. Hiramatsu and M. Unoki: A Study on the Blind Estimation of Reverberation Time in Room Acoustics, J.
Signal Processing, Vol. 12, No. 4, pp. 323–326, July 2008.
[14] S. Irie and S. Hirobayashi: Simulation of an acoustic system using power envelope inverse ﬁltering, Proc.
ICASSP2008, pp. 325–328, Las Vegas, 2008.
[15] A. Papouris: Probability: Random Variables, andStochastic Processes, Third Edition, MacGraw-Hill, Inc.,New York, 1991.
[16] M. R. Schroeder: Modulation Transfer Functions: Deﬁnition and Measurement, Acustica, Vol. 49, pp. 179–182,
1981.
[17] M. R. Schroeder: New Method of Measuring Reverberation Time, J. Acoust. Soc. Am., Vol. 37, No. 6, pp. 409,
1965.
[18] 広林茂樹, 寺島洋行，山淵龍夫: 残響音場における音響信号のエンベロープ推定法の評価, シミュレーション, Vol. 22, No. 3,
pp. 68–75 (208–215), 2003.
[19] X. Lu, M. Unoki, and M. Akagi: A study on restoration methods for temporal envelope of reverberant speechbased on the modulation transfer function concept, Proc.
JCA2007, P-2-11, Sendai, Japan, 2007.
[20] F. J. Taylor: Principles of Signals and Systems, MacGrawHill, Inc., New York, 1994.
[21] A. V. Oppenheim and R. W. Schafer: Digital Signal Processing, Prentice-Hall, Inc., London, 1975.

付録１: 式 (8) と式 (9) の導出

v(t) をランダム変数，f (v, t) を v(t) の確率密度関数とする。このとき，集合平均は次式で定義される [15]。
∫ ∞< v(t) >=vf (v, t)dv
(27)
−∞

10

Journal of Signal Processing, Vol. , No. , 2008

式 (7) の関係に注意して，式 (8) の左辺を求めると< x2 (t) > ===−∞< e2x (t)n2x (t) >e2x (t) < n2x (t) >e2x (t)付録２: 式 (11) の導出
式 (10) の < y2 (t) > は x(t) と h(t) の畳み込みを 2乗した集合平均と等しいため，次式となる。
⟨{∫} ⟩< y(t)2 >=−∞ex (τ1 )eh (t − τ1 )∫ ∞−∞eh (τ2 )eh (t − τ2 )×< nx (τ1 )nx (τ2 ) >×< nh (t − τ1 )nh (t − τ2 ) >dτ1 dτ2を得る。但し，< n2x (t) >= δ(0) = 1 である。式 (9) の< y2 (t) > についても同様の導出過程から得られる。

∞

∫ ∞
=

2

h(τ )x(t − τ )dτ

(28)を得る。但し，nx (t) と nh (t) は，平均 0，分散 1 を有するランダム関数であり，式 (7) の関係をもつ。ここで，τ = τ1 = τ2 と式 (9) を代入すると，次式を得る。

< y2 (t) >
∫ ∞
∫ ∞
=
ex (τ1 )eh (t − τ1 )ex (τ2 )eh (t − τ2 )−∞−∞×δ(τ2 − τ1 ) dτ1 dτ2
2

∫ ∞
=
=

e2x (τ )e2h (t − τ )dτ
−∞e2x (t) ∗ e2h (t) = e2y (t)これを整理すると
付録４: 式 (21) の導出< y2 (t) >∫ ∞=h2 (τ )e2x (t − τ ) < nx (t − τ )nx (t − τ ) > dτ∫0 ∞=h2 (τ )(e2x (1 + cos(ωm (t − τ )))dτ
0
[∫ ∞
]
∫ ∞
= e2x
h2 (τ )dτ +h2 (τ ) cos(ωm (t − τ )))dτ
0
[∫0 ∞
{∫ ∞
}]
h2 (τ )dτ + R
h2 (τ )ejωm (t−τ ) dτ
= e2x
0
[∫0 ∞
= e2x
h2 (τ )dτ
0
{
}]
∫ ∞
j2πfm t2−j2πfm τ+R eh (τ )e
dτ
0

=

[
]
β
1 + cos(2πfm t)
α
α

e2x

∫∞
を得る。但し，R は実数部を意味し，α = 0 h2 (τ )dτ ，∫∞ 2
β = 0 h (τ ) exp(−jωm τ )dτ である。ここで，式 (4)との対応を見れば，m(fm ) = β/α から，式 (12) を得る。

パワーエンベロープ e2h [n] の z 変換は

Eh (z) =∞∑e2h [n]z −nn=−∞である。但し，t = n/fs であり，t < 0 で e2h (t) = 0 であるため，n の下限は 0 である。これを整理すると，
)
(
∞
∑
13.8n
2
z −n
Eh (z) =a exp −
fs · TRn=0= a2 γ n z −n
= a2 + a2 γz −1 + a2 γ 2 z −2 + a2 γ 3 z −3 + · · ·a2=
1 − γz −1a2()=
1 − exp − f13.8
z −1
s ·TR
を得る。但し，γ = exp(−13.8/fs · TR ) である。

付録３: 式 (19) の導出
2

式 (19) の < y (t) > は x(t) と h(t) の畳み込みを 2乗した集合平均と等しいため，次式となる。
⟨{∫} ⟩< y(t)2 >=

∞

−∞

2

x(τ )h(t − τ )dτ

ここで，式 (5) と式 (13) を代入し，整理すると< y2 (t) >⟨∫ ∞=ex (τ1 )nx (τ1 )eh (t − τ1 )nh (t − τ1 )dτ1
−∞
⟩
∫ ∞
×
ex (τ2 )nx (τ2 )eh (t − τ2 )nh (t − τ2 )dτ2
−∞

鵜木 祐史
1994 年職業能力開発大学校情報工学科卒。1996 年北陸先端科
学技術大学院大学情報科学研究科博士前期課程修了。1999 年同大情報科学研
究科博士後期課程修了。博士（情報科学）。同年 ATR 人間情報通信研究所第一研究室客員研究員，2000 年英国ケンブリッジ大学生理学部 CNBH 客員研究員，2001 年北陸先端科学技術大学院大学情報科学研究科助手を経て，2005 年
同大助教授に奉職（2007 年に准教授）。
現在に至る。1998 年〜2001 年の間，日本学術振興会特別研究員（DC2, PD の2 期）を兼任。主に，聴覚機能のモデル化とそれに基づく信号処理
ならびに音声信号処理（残響音声回復，骨導音声回復，音響電子透かし）の研究に従事。日本音響学会佐藤論文賞（1999 年度）ならびに山下太郎学術奨励賞（2005 年度）受賞。信号処理学会，日本電子
情報通信学会，日本音響学会，アメリカ音響学会，IEEE，ISCA，EURASHIP 各会員。