(57)【要約】【課題】簡便な処理で音響信号の残響成分を調整する。【解決手段】指標算定部５２は、音響信号の強度の移動平均で指標値ＲA(k,m)と指標値ＲB(k,m)とを算定する。指標値ＲB(k,m)は、指標値ＲA(k,m)を上回る追従性で音響信号の強度に追従する。調整値算定部５４は、音響信号の残響成分を抑圧するための調整値Ｇ(k,m)を算定する。具体的には、調整値算定部５４は、指標値ＲA(k,m)と指標値ＲB(k,m)とを分母に包含するとともに指標値ＲB(k,m)を分子に包含する演算式を利用して調整値Ｇ(k,m)を算定する。【選択図】図２
2015004959.tif 000002
【特許請求の範囲】
【請求項１】
音響信号の強度の移動平均で第１指標値を算定する指標算定手段と、
前記音響信号の残響成分を調整するための調整値を算定する手段であって、前記第１指標値と前記第１指標値を上回る追従性で前記音響信号の強度に追従する第２指標値とを分母に包含し、前記第１指標値を上回る追従性で前記音響信号の強度に追従する第３指標値を分子に包含する演算式を利用して前記調整値を算定する調整値算定手段と
を具備する音響処理装置。
【請求項２】
前記第２指標値と前記第３指標値とは共通する
請求項１の音響処理装置。
【請求項３】
加重係数を設定する係数設定手段を具備し、
前記調整値算定手段は、前記第１指標値と、前記係数設定手段が設定した加重係数を乗算した前記第２指標値との加算を分母に包含する前記演算式を利用して前記調整値を算定する
請求項１または請求項２の音響処理装置。
【請求項４】
前記係数設定手段は、前記第１指標値を分母に包含し、前記第２指標値または前記第３指標値を分子に包含する演算式の演算結果に応じて前記加重係数を設定する
請求項３の音響処理装置。
【請求項５】
前記調整値算定手段は、前記第１指標値と前記第２指標値とを分母に包含するとともに前記第３指標値を分子に包含する演算で基礎値を順次に算定し、前記基礎値が第１値を上回る数値に変動した場合に、以降の各基礎値を当該数値に応じて低減することで前記調整値を算定する
請求項１から請求項４の何れかの音響処理装置。
【請求項６】
前記調整値算定手段は、前記第１指標値と前記第２指標値との加算を分母に包含するとともに前記第３指標値を分子に包含する演算で算定される第１基礎値と、前記第１指標値に対する前記第３指標値の相対比である第２基礎値との加重和を前記調整値として算定する
請求項１から請求項５の何れかの音響処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号を処理する技術に関し、特に、音響信号に含まれる残響成分を調整する技術に関する。
【背景技術】
【０００２】
音響信号に含まれる残響成分を抑圧する技術が従来から提案されている。例えば特許文献１には、音響信号に含まれる残響成分を推定する予測フィルタ係数の確率モデルを利用することで残響成分の予測フィルタ係数を推定し、推定後の予測フィルタを利用して残響成分を抑圧する技術が開示されている。また、非特許文献１には、発音源から収音点までの伝達関数の逆フィルタを推定し、推定後の逆フィルタを音響信号に適用することで残響成分を抑圧する技術が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００９−２１２５９９号公報
【非特許文献】
【０００４】
【非特許文献１】K. Furuya, et al."Robust speech dereverberation using multichannel blind deconvolution with spectral subtraction",IEEE Transantions on Audio, Speech, and Language Processing, vol. 15, no. 5, p.1579-1591, 2007
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、特許文献１の予測フィルタ係数や非特許文献１の逆フィルタの高精度な推定には膨大な演算が必要であるという問題がある。以上の事情を考慮して、本発明は、簡便な処理で音響信号の残響成分を調整（抑圧または強調）することを目的とする。
【課題を解決するための手段】
【０００６】
以上の課題を解決するために、本発明の音響処理装置は、音響信号の強度の移動平均で第１指標値（例えば指標値ＲA(k,m)）を算定する指標算定手段と、音響信号の残響成分を調整するための調整値を算定する手段であって、第１指標値と第１指標値を上回る追従性で音響信号の強度に追従する第２指標値（例えば指標値ＲB(k,m)）とを分母に包含し、第１指標値を上回る追従性で音響信号の強度に追従する第３指標値（例えば指標値ＲC(k,m)）を分子に包含する演算式を利用して調整値を算定する調整値算定手段とを具備する。以上の構成では、音響信号の強度の移動平均で調整値が算定されるから、簡便な処理で音響信号の残響成分を調整することが可能である。また、第１指標値を上回る追従性で音響信号の強度に追従する第２指標値を分母に含む演算で調整値が算定されるから、音響信号の残響時間が長い場合でも残響成分を高精度に調整できるという利点がある。
【０００７】
本発明の好適な態様において、第２指標値と第３指標値とは共通する。以上の態様では、第２指標値と第３指標値とが共通するから、第２指標値と第３指標値とを個別に算定する構成と比較して処理負荷が軽減されるという利点がある。ただし、第２指標値と第３指標値とを個別に算定する構成も採用され得る。
【０００８】
本発明の好適な態様に係る音響処理装置は、加重係数を設定する係数設定手段を具備し、調整値算定手段は、第１指標値と、係数設定手段が設定した加重係数を乗算した第２指標値との加算を分母に包含する演算式を利用して調整値を算定する。以上の態様では、加重係数が第２指標値に乗算されるから、演算式の分母における第２指標値の効果を加重係数に応じて調整することが可能である。例えば、残響時間に対して単調増加するように加重係数を算定する構成（例えば第２実施形態）や、前記第１指標値を分母に包含するとともに前記第２指標値または前記第３指標値を分子に包含する演算式の演算結果に応じて前記加重係数を設定する構成（例えば第７実施形態）によれば、残響時間が長い場合でも残響成分を高精度に抑圧できるという前述の作用が強調される。
【０００９】
本発明の好適な態様において、調整値算定手段は、第１指標値と第２指標値とを分母に包含するとともに第３指標値を分子に包含する演算で基礎値を順次に算定し、基礎値が第１値（例えば１）を上回る数値に変動した場合に、以降の各基礎値を当該数値に応じて低減することで調整値を算定する。以上の態様では、基礎値が第１値を上回る場合に以降の各基礎値が低減されるから、音響信号の音響の発音が開始された直後の残響成分（例えば初期反射成分）も高精度に抑圧できるという利点がある。また、基礎値が第１値未満の第２値（例えば閾値ＱTH）を下回ると、調整値の低減を終了する構成によれば、音響信号の強度の過剰な抑圧を防止することが可能である。なお、以上の態様の具体例は、例えば第３実施形態として後述される。
【００１０】
本発明の好適な態様において、調整値算定手段は、第１指標値と第２指標値との加算を分母に包含するとともに第３指標値を分子に包含する演算で算定される第１基礎値と、第１指標値に対する第３指標値の相対比である第２基礎値との加重和を調整値として算定する。以上の態様では、第１基礎値と第２基礎値との加重和が調整値として算定されるから、第１基礎値のみに応じて調整値を算定する構成と比較して残響成分の高精度な調整が期待される。また、音響信号に対する調整値の作用が、音響信号に付与された残響効果の振幅-周波数特性の逆特性に近似するように、第１基礎値および第２基礎値の各々の加重値を算定する構成によれば、残響成分の高精度な調整という効果は格別に顕著である。なお、振幅-周波数特性は、例えば変調伝達関数（例えば変調伝達関数Ψ(ω)）で近似的に表現される。以上の態様の具体例は、例えば第４実施形態として後述される。
【００１１】
残響成分を調整する従来の技術では、残響時間が長い場合に残響成分（特に目的成分の発生から充分な時間が経過した段階の残響成分）の調整精度が低下するという問題がある。以上の事情を考慮して、本発明の他の態様に係る音響処理装置は、音響信号の残響成分の強度を推定する残響推定手段と、音響信号の残響成分を調整するための調整値を算定する手段であって、音響信号の強度と残響推定手段が推定した残響成分の強度とを分母に包含するとともに音響信号の強度を分子に包含する演算式を利用して調整値を算定する調整値算定手段とを具備する。以上の態様では、音響信号の強度に加えて残響成分の強度を分母に包含する演算式を利用して調整値が算定されるから、残響時間が長い場合でも音響信号の残響成分を高精度に調整できるという利点がある。なお、以上の態様の具体例は、例えば第６実施形態として後述される。
【００１２】
以上の各態様に係る音響処理装置は、音響信号の処理に専用されるＤＳＰ（Digital Signal Processor）などのハードウェア（電子回路）によって実現されるほか、ＣＰＵ（Central Processing Unit）等の汎用の演算処理装置とプログラムとの協働によっても実現される。コンピュータが読取可能な記録媒体に格納された形態で提供されてコンピュータにインストールされ得る。記録媒体は、例えば非一過性（non-transitory）の記録媒体であり、ＣＤ-ＲＯＭ等の光学式記録媒体（光ディスク）が好例であるが、半導体記録媒体や磁気記録媒体等の公知の任意の形式の記録媒体を包含し得る。また、例えば、本発明のプログラムは、通信網を介した配信の形態で提供されてコンピュータにインストールされ得る。
【図面の簡単な説明】
【００１３】
【図１】本発明の第１実施形態に係る音響処理装置の構成図である。
【図２】解析処理部の構成図である。
【図３】各指標値と調整値との関係の説明図である。
【図４】演算処理装置の動作のフローチャートである。
【図５】第２実施形態に係る音響処理装置の構成図である。
【図６】第３実施形態に係る音響処理装置の動作のフローチャートである。
【図７】第３実施形態における基礎値の調整の説明図である。
【図８】第６実施形態における解析処理部の構成図である。
【図９】第７実施形態に係る音響処理装置の構成図である。
【図１０】第７実施形態における変数と加重係数との関係の説明図である。
【図１１】第７実施形態に調整値の時間変化と残響時間との関係の説明図である。
【図１２】変形例における音響処理装置の部分的な構成図である。
【発明を実施するための形態】
【００１４】
＜第１実施形態＞
図１は、本発明の第１実施形態に係る音響処理装置１００の構成図である。図１に例示される通り、音響処理装置１００には信号供給装置１２と放音装置１４とが接続される。信号供給装置１２は、音響信号ｘを音響処理装置１００に供給する。音響信号ｘは、発音源から放射された音響に対して音響空間内で反射または散乱した残響成分（初期反射成分および後期残響成分）を付加した音響の時間波形を示すサンプル系列である。例えば、収録音や合成音等の既存の音響に対して事後的に残響効果を付与した音響信号ｘや、残響効果がある音響空間（例えば音響ホール等）内で実際に収録された音響の音響信号ｘが好適に利用される。周囲の音響を収音して音響信号ｘを生成する収音装置や、可搬型または内蔵型の記録媒体から音響信号ｘを取得して音響処理装置１００に供給する再生装置や、通信網から音響信号ｘを受信して音響処理装置１００に供給する通信装置が信号供給装置１２として採用され得る。
【００１５】
第１実施形態の音響処理装置１００は、音響信号ｘの残響成分（後期残響成分）を抑圧した時間領域の音響信号ｙを生成する残響抑圧装置である。放音装置１４（例えばスピーカやヘッドホン）は、音響処理装置１００が生成した音響信号ｙに応じた音波を再生する。なお、音響信号ｙをデジタルからアナログに変換するＤ/Ａ変換器や音響信号ｙを増幅する増幅器等の図示は便宜的に省略した。
【００１６】
図１に例示される通り、音響処理装置１００は、演算処理装置２２と記憶装置２４とを具備するコンピュータシステムで実現される。記憶装置２４は、演算処理装置２２が実行するプログラムや演算処理装置２２が使用する各種のデータを記憶する。半導体記録媒体や磁気記録媒体等の公知の記録媒体または複数種の記録媒体の組合せが記憶装置２４として任意に採用され得る。音響信号ｘを記憶装置２４に記憶した構成（したがって信号供給装置１２は省略される）も好適である。
【００１７】
演算処理装置２２は、記憶装置２４に記憶されたプログラムを実行することで、音響信号ｘから音響信号ｙを生成するための複数の機能（周波数分析部３２，残響調整部３４，波形生成部３６，解析処理部４２）を実現する。なお、演算処理装置２２の各機能を複数の装置に分散した構成や、専用の電子回路（例えばＤＳＰ）が演算処理装置２２の一部の機能を実現する構成も採用され得る。
【００１８】
周波数分析部３２は、音響信号ｘを構成する複数の周波数成分（周波数スペクトル）Ｘ(k,m)を時間軸上の単位期間（フレーム）毎に順次に生成する。記号ｋは、周波数軸上に離散的に設定された複数の周波数（帯域）のうち任意の１個の周波数を指定する変数であり、記号ｍは、時間軸上の任意の１個の単位期間（時間軸上の特定の時点）を指定する変数である。各周波数成分Ｘ(k,m)の算定には、短時間フーリエ変換等の公知の周波数解析が任意に採用され得る。
【００１９】
解析処理部４２は、音響信号ｘの各周波数成分Ｘ(k,m)に応じた調整値Ｇ(k,m)を各周波数成分Ｘ(k,m)について単位期間毎に算定する。第１実施形態の調整値Ｇ(k,m)は、音響信号ｘの残響成分を抑圧するための変数である。概略的には、第ｍ番目の単位期間の音響信号ｘのうち周波数軸上の第ｋ番目の周波数成分Ｘ(k,m)において残響成分が優勢であるほど調整値Ｇ(k,m)は小さい数値に設定されるという傾向がある。
【００２０】
図１の残響調整部３４は、解析処理部４２が算定した各調整値Ｇ(k,m)を音響信号ｘに作用させる。具体的には、残響調整部３４は、各周波数成分Ｘ(k,m)について算定された調整値Ｇ(k,m)を当該周波数成分Ｘ(k,m)に乗算することで音響信号ｙの各周波数成分（周波数スペクトル）Ｙ(k,m)を算定する（Ｙ(k,m)＝Ｇ(k,m)・Ｘ(k,m)）。以上の説明から理解される通り、調整値Ｇ(k,m)は、音響信号ｘの周波数成分Ｘ(k,m)に対するゲイン（スペクトルゲイン）に相当する。
【００２１】
波形生成部３６は、残響調整部３４が算定する各周波数成分Ｙ(k,m)から時間領域の音響信号ｙを生成する。すなわち、波形生成部３６は、各周波数成分Ｙ(k,m)を単位期間毎に短時間逆フーリエ変換で時間領域の信号に変換し、相前後する各単位期間の信号を相互に重複させた状態で加算することで音響信号ｙを生成する。波形生成部３６が生成した音響信号ｙが放音装置１４に供給されて音波として再生される。
【００２２】
図２は、解析処理部４２の構成図である。図２に例示される通り、解析処理部４２は、指標算定部５２と調整値算定部５４とを含んで構成される。指標算定部５２は、音響信号ｘの各周波数成分Ｘ(k,m)に応じた指標値ＲA(k,m)および指標値ＲB(k,m)を単位期間毎に順次に算定する。第１実施形態の指標算定部５２は、音響信号ｘの各周波数成分Ｘ(k,m)の強度Ｐ(k,m)の時系列を平滑化することで各周波数の指標値ＲA(k,m)および指標値ＲB(k,m)を算定する。強度Ｐ(k,m)は、例えば周波数成分Ｘ(k,m)のパワー|Ｘ(k,m)|2である。具体的には、以下の数式(1A)および数式(1B)で表現される通り、音響信号ｘの強度Ｐ(k,m)の指数移動平均が指標値ＲA(k,m)および指標値ＲB(k,m)として算定される。
【数１】
2015004959.tif 000003
【００２３】
数式(1A)の記号αAおよび数式(1B)の記号αBは、指数移動平均の平滑化係数（すなわち、強度Ｐ(k,m)の平滑化の時定数を規定する係数）であり、１未満の正数に設定される。平滑化係数αAは、過去の指標値ＲA(k,m-1)に対する最新（現在）の強度Ｐ(k,m)の加重値に相当し、平滑化係数αBは、過去の指標値ＲB(k,m-1)に対する最新（現在）の強度Ｐ(k,m)の加重値に相当する。指標値ＲB(k,m)の平滑化係数αBは、指標値ＲA(k,m)の平滑化係数αAを上回る（αB＞αA）。したがって、指標値ＲA(k,m)における強度Ｐ(k,m)の平滑化の時定数τAは、指標値ＲB(k,m)における強度Ｐ(k,m)の平滑化の時定数τBを上回る（τA＞τB）。
【００２４】
図３の部分(B)には、音響信号ｘの１個の周波数の周波数成分Ｘ(k,m)から算定される指標値ＲA(k,m)および指標値ＲB(k,m)の時間変化が例示されている。図３の部分(A)のように強度Ｐ(k,m)が指数減衰する室内インパルス応答（ＲＩＲ：Room Impulse Response）を音響信号ｘとして音響処理装置１００に供給した場合の指標値ＲA(k,m)および指標値ＲB(k,m)が図３の部分(B)には図示されている。
【００２５】
図３の部分(B)から理解される通り、指標値ＲA(k,m)および指標値ＲB(k,m)は、音響信号ｘの強度Ｐ(k,m)に追従して経時的に変動する。ただし、前述の通り、指標値ＲA(k,m)の時定数τAは指標値ＲB(k,m)の時定数τBを上回るから、指標値ＲB(k,m)は、指標値ＲA(k,m)と比較して高い追従性で音響信号ｘの強度Ｐ(k,m)の時間変化に追従する。具体的には、図３の部分(B)から把握される通り、室内インパルス応答の開始の時点ｔ0の直後の区間では、指標値ＲB(k,m)が指標値ＲA(k,m)を上回る変化率で急峻に増加する。そして、指標値ＲA(k,m)および指標値ＲB(k,m)は、時間軸上の相異なる時点で極大値に到達し、指標値ＲB(k,m)は指標値ＲA(k,m)を上回る変化率で減少する。
【００２６】
以上のように指標値ＲA(k,m)と指標値ＲB(k,m)とは相異なる変化率で変化するから、指標値ＲA(k,m)と指標値ＲB(k,m)との大小は時間軸上の特定の時点ｔxで反転する。すなわち、時点ｔ0から時点ｔxまでの区間ＳAでは指標値ＲB(k,m)が指標値ＲA(k,m)を上回り、時点ｔx以降の区間ＳBでは指標値ＲA(k,m)が指標値ＲB(k,m)を上回る。区間ＳAは、室内インパルス応答のうち直接音（目的成分）および初期反射音が存在する区間に相当し、区間ＳBは、室内インパルス応答のうち後期残響成分が存在する区間に相当する。
【００２７】
図２の調整値算定部５４は、指標算定部５２が算定した指標値ＲA(k,m)と指標値ＲB(k,m)とに応じた調整値Ｇ(k,m)を各周波数について単位期間毎に算定する。第１実施形態の調整値Ｇ(k,m)の具体例の説明に先立ち、指標値ＲA(k,m)に対する指標値ＲB(k,m)の相対比ＲB(k,m)/ＲA(k,m)を調整値Ｇ(k,m)として利用する場合を第１実施形態の対比例として説明する。なお、相対比ＲB(k,m)/ＲA(k,m)が１を上回る場合、調整値Ｇ(k,m)は１に設定される。
【００２８】
指標値ＲA(k,m)と指標値ＲB(k,m)とが図３の部分(B)のように変動する場合、対比例の調整値Ｇ(k,m)は図３の部分(C)の例示のように経時的に変動する。図３の部分(C)から理解される通り、指標値ＲB(k,m)が指標値ＲA(k,m)を上回る区間ＳA（直接音および初期反射音が存在する区間）では調整値Ｇ(k,m)は最大値１に設定され、指標値ＲA(k,m)が指標値ＲB(k,m)を上回る区間ＳB（後期残響成分が存在する区間）では調整値Ｇ(k,m)は経時的に減少する。したがって、対比例の調整値Ｇ(k,m)を音響信号ｘの各周波数成分Ｘ(k,m)に乗算することで、音響信号ｘの残響成分を抑圧した音響信号ｙを生成することが可能である。
【００２９】
ところで、音響信号ｘの残響時間が比較的に短い場合には指標値ＲB(k,m)が短時間で急峻に減少するから、対比例の調整値Ｇ(k,m)は迅速に充分に小さい数値に減少する。したがって、対比例でも音響信号ｘの残響成分を抑圧することが可能である。しかし、音響信号ｘの残響時間が比較的に長い場合には指標値ＲB(k,m)が長時間にわたり大きい数値に維持されるから、音響信号ｘの残響成分を有効に抑圧し得る充分に小さい数値まで調整値Ｇ(k,m)が減少するには長時間が必要である。したがって、対比例では、音響信号ｘの残響時間が長い場合に残響成分の抑圧精度が低下する可能性がある。
【００３０】
以上の事情を考慮して、第１実施形態の調整値算定部５４は、以下の数式(2)で表現される基礎値Ｑ(k,m)を調整値Ｇ(k,m)として算定する。
【数２】
2015004959.tif 000004
数式(2)から理解される通り、指標値ＲA(k,m)と指標値ＲB(k,m)との加算に対する指標値ＲB(k,m)の相対比である基礎値Ｑ(k,m)が調整値Ｇ(k,m)として算定される（Ｇ(k,m)＝Ｑ(k,m)）。すなわち、調整値算定部５４は、指標値ＲA(k,m)と指標値ＲB(k,m)との双方を分母に包含するとともに指標値ＲB(k,m)を分子に包含する数式(2)を利用して調整値Ｇ(k,m)を算定する。指標値ＲA(k,m)および指標値ＲB(k,m)は非負値であるから、調整値Ｇ(k,m)（基礎値Ｑ(k,m)）は１以下の範囲内で変動する。
【００３１】
前述の通り、音響信号ｘの残響時間が比較的に長い場合には指標値ＲB(k,m)が長時間にわたり大きい数値に維持されるから、数式(2)の分母の指標値ＲB(k,m)は、残響時間が長い場合に、調整値Ｇ(k,m)（基礎値Ｑ(k,m)）を減少させるように作用する。したがって、数式(2)で算定される基礎値Ｑ(k,m)を調整値Ｇ(k,m)として音響信号ｘの各周波数成分Ｘ(k,m)に作用させることで、残響時間が長い場合でも音響信号ｘの残響成分を高精度に抑圧することが可能である。
【００３２】
なお、音響信号ｘのうち残響成分以外の目的成分が優勢な区間（例えば直接音や初期反射音が存在する区間ＳA）では、指標値ＲB(k,m)が指標値ＲA(k,m)と比較して充分に大きい（ＲB(k,m)≫ＲA(k,m)）から、数式(2)で算定される調整値Ｇ(k,m)は１（最大値）に近い数値となる。したがって、数式(2)の分母の指標値ＲB(k,m)の作用で調整値Ｇ(k,m)が減少するとは言え、調整値Ｇ(k,m)が過度に小さい数値に抑制されること（音響信号ｘの目的成分が過剰に抑圧されること）は防止される。すなわち、音響信号ｘの残響成分を高精度に抑圧することが可能である。
【００３３】
図４は、記憶装置２４に記憶されたプログラムに従って演算処理装置２２が実行する残響調整処理のフローチャートである。入力装置（図示略）に対する利用者からの指示（残響調整処理の開始指示）を契機として図４の残響調整処理が開始される。残響調整処理を開始すると、演算処理装置２２は、単位期間毎に以下のステップＳ1からステップＳ5の処理を実行する。まず、周波数分析部３２は、音響信号ｘの各周波数の周波数成分Ｘ(k,m)を算定する（Ｓ1）。
【００３４】
解析処理部４２の指標算定部５２は、周波数分析部３２が算定した各周波数成分Ｘ(k,m)を適用した数式(1A)および数式(1B)の演算で指標値ＲA(k,m)および指標値ＲB(k,m)を算定する（Ｓ2）。調整値算定部５４は、指標算定部５２が算定した指標値ＲA(k,m)および指標値ＲB(k,m)を適用した数式(2)の演算で各周波数の調整値Ｇ(k,m)を算定する（Ｓ3）。
【００３５】
残響調整部３４は、調整値算定部５４がステップＳ3で算定した調整値Ｇ(k,m)を音響信号ｘの各周波数成分Ｘ(k,m)に作用させることで周波数成分Ｙ(k,m)を算定する（Ｓ4）。波形生成部３６は、残響調整部３４が算定した各周波数成分Ｙ(k,m)から単位期間の音響信号ｙを生成する（Ｓ5）。音響信号ｘの全部の単位期間について以上の処理（Ｓ1〜Ｓ5）が反復される（Ｓ6：NO）ことで、音響信号ｘの残響成分を抑圧した音響信号ｙが生成される。
【００３６】
以上に説明した通り、第１実施形態では、音響信号ｘの強度Ｐ(k,m)の移動平均である指標値ＲA(k,m)および指標値ＲB(k,m)に応じて調整値Ｇ(k,m)が算定されるから、特許文献１や非特許文献１の技術と比較して簡便に音響信号ｘの残響成分を抑圧できるという利点がある。
【００３７】
＜第２実施形態＞
本発明の第２実施形態を以下に説明する。なお、以下に例示する各構成において作用や機能が第１実施形態と同様である要素については、第１実施形態の説明で参照した符号を流用して各々の詳細な説明を適宜に省略する。
【００３８】
図５は、第２実施形態における音響処理装置１００の構成図である。図５に例示される通り、第２実施形態の演算処理装置２２は、第１実施形態と同様の要素（周波数分析部３２，残響調整部３４，波形生成部３６，解析処理部４２）に加えて係数設定部４４として機能する。係数設定部４４は、調整値算定部５４が調整値Ｇ(k,m)の算定に適用する係数（以下「加重係数」という）ηを設定する。加重係数ηは正数に設定される。
【００３９】
解析処理部４２の調整値算定部５４は、前掲の数式(2)に代えて、係数設定部４４が設定した加重係数ηを適用した以下の数式(3)の演算で基礎値Ｑ(k,m)を算定する。
【数３】
2015004959.tif 000005
【００４０】
数式(3)から理解される通り、例えば加重係数ηが１を下回る場合、基礎値Ｑ(k,m)は１以上の数値となり得る。そこで、調整値算定部５４は、数式(3)で算定した基礎値Ｑ(k,m)を利用した数式(4)の演算で調整値Ｇ(k,m)を算定する。
【数４】
2015004959.tif 000006
すなわち、調整値Ｇ(k,m)は、１を上限値として基礎値Ｑ(k,m)に設定される。以上に説明した通り、第２実施形態においても指標値ＲA(k,m)および指標値ＲB(k,m)に応じて調整値Ｇ(k,m)が算定されるから、第１実施形態と同様に、音響信号ｘの残響成分を簡便に抑圧することが可能である。なお、数式(4)では調整値Ｇ(k,m)の上限値を１としたが、調整値Ｇ(k,m)の上限値は任意の数値に設定され得る。例えば、いわゆる桁溢れを防止する観点から、１を下回る正数（例えば０.９９）を調整値Ｇ(k,m)の上限値とすることも可能である。
【００４１】
数式(3)から理解される通り、第２実施形態の調整値算定部５４は、指標値ＲA(k,m)と、加重係数ηを乗算した指標値ＲB(k,m)との加算を分母に包含する数式(3)を利用して調整値Ｇ(k,m)を算定する。したがって、加重係数ηが大きいほど調整値Ｇ(k,m)は小さい数値に設定される。すなわち、数式(3)の分母の指標値ＲB(k,m)が調整値Ｇ(k,m)を減少させる作用は、当該指標値ＲB(k,m)に乗算される加重係数ηが大きいほど強調される。
【００４２】
図５の係数設定部４４は、音響信号ｘの残響成分に想定される残響時間（例えば音響信号ｘが発音および収録された音響空間の残響時間）ＴRに応じて加重係数ηを可変に設定する。具体的には、係数設定部４４は、例えば入力装置（図示略）に対する利用者からの指示（残響時間ＴRの指定）に応じて残響時間ＴRを設定し、残響時間ＴRに応じた加重係数ηを設定する。
【００４３】
第２実施形態の係数設定部４４は、残響時間ＴRに対して単調増加するように加重係数ηを設定する。すなわち、残響時間ＴRが長いほど加重係数ηは大きい数値に設定される。例えば、係数設定部４４は、以下の数式(5)の演算で残響時間ＴRに応じた加重係数ηを算定する。
【数５】
2015004959.tif 000007
数式(5)の係数ａは所定の正数に設定され、係数ｂは所定値に設定される。以上の説明から理解される通り、残響時間ＴRが長い（加重係数ηが大きい）ほど、数式(3)の分母の指標値ＲB(k,m)が調整値Ｇ(k,m)を減少させる作用は強調される。したがって、第２実施形態では、残響時間が長い場合でも音響信号ｘの残響成分を高精度に抑圧できるという前述の効果は格別に顕著である。
【００４４】
なお、残響時間ＴRに応じた加重係数ηを算定する演算は数式(5)に限定されない。例えば、以下の数式(6-1)〜(6-3)の演算で係数設定部４４が加重係数ηを算定することも可能である。以下の各数式における係数ａおよび係数ｂは所定値に設定される。
【数６】
2015004959.tif 000008
数式(6-1)の加重係数ηは、残響時間ＴRのγ次関数として表現され、冪指数γは所定の整数（奇偶は不問）に設定される。数式(6-2)の加重係数ηは、残響時間ＴRの指数関数として表現され、数式(6-3)の加重係数ηは、残響時間ＴRの対数関数として表現される。数式(6-3)では、変数｛ａＴR＋ｂ｝が１以上となるように各係数を選定することで加重係数ηを正数の範囲内で設定することが可能である。なお、加重係数ηを、残響時間ＴRに依存しない所定値に設定することも可能である。
【００４５】
また、前述の説明では、利用者が残響時間ＴRを指示する場合を例示したが、音響空間内の音響の収音で音響信号ｘを生成する収音装置を信号供給装置１２として採用した構成では、所定の測定用信号を音響空間内に放音したときに収録される音響信号ｘを解析することで係数設定部４４が残響時間ＴRを実測することも可能である。残響時間ＴRの測定には、例えばインパルス積分法（Schroeder法）やマルチステップ線形予測等の公知の残響解析技術が任意に採用され得る。また、周波数が時間的に連続に変化する時間伸長信号（ＴＳＰ：Time Stretched Pulse）が測定用信号として好適である。
【００４６】
＜第３実施形態＞
前掲の数式(3)の基礎値Ｑ(k,m)が１を上回る区間では調整値Ｇ(k,m)が１に維持されるから、周波数成分Ｘ(k,m)の残響成分（特に初期反射音）は抑圧されず、結果的に抑圧効果が不足する可能性がある。以上の事情を考慮して、第３実施形態では、基礎値Ｑ(k,m)が１を上回る場合に当該基礎値Ｑ(k,m)を低減（正規化）することで、調整値Ｇ(k,m)が１に制限される時間を短縮する。
【００４７】
図６は、第３実施形態の調整値算定部５４が調整値Ｇ(k,m)を算定する処理のフローチャートである。前掲の図４のステップＳ3で図６の処理が実行される。したがって、音響信号ｘの単位期間毎に図６の処理が実行されて、当該単位期間の調整値Ｇ(k,m)が周波数毎に算定される。なお、第１回目の図６の処理が実行される直前に、基礎値Ｑ(k,m)の調整に適用される変数（以下「正規化係数」という）ＧNが初期値１に設定される。
【００４８】
図６の処理を開始すると、調整値算定部５４は、前掲の数式(3)の演算で今回の単位期間の基礎値Ｑ(k,m)を算定する（Ｓ31）。そして、調整値算定部５４は、今回の単位期間の基礎値Ｑ(k,m)が１を上回るか否かを判定する（Ｓ32）。基礎値Ｑ(k,m)が１を上回る場合（Ｓ32：YES）、調整値算定部５４は、現段階の正規化係数ＧNが１であるか否かを判定する（Ｓ33）。基礎値Ｑ(k,m)の調整中には正規化係数ＧNは１を上回る数値に設定される。したがって、ステップＳ33の判定は、基礎値Ｑ(k,m)の正規化を実行している最中であるか否かの判定に相当する。
【００４９】
正規化係数ＧNが１である場合（Ｓ33：YES）、調整値算定部５４は、正規化係数ＧNを、今回の基礎値Ｑ(k,m)（Ｑ(k,m)＞１）に更新する（Ｓ34）。すなわち、以降の各単位期間の基礎値Ｑ(k,m)を正規化するための数値に正規化係数ＧNが設定される。他方、基礎値Ｑ(k,m)が１以下である場合（Ｓ32：NO）や正規化係数ＧNが１である場合（Ｓ33：NO）には、正規化係数ＧNの更新（Ｓ34）は実行されない。
【００５０】
図７は、調整前の基礎値Ｑ(k,m)（実線）と調整後の基礎値Ｑ(k,m)（破線）との関係を例示するグラフである。図７から理解される通り、単位期間毎に順次に算定される複数の基礎値Ｑ(k,m)のうち最初に１を上回る数値に到達した時点ｔ1の基礎値Ｑ(k,m)が（Ｓ32：YES，Ｓ33：YES）、以降の各単位期間の基礎値Ｑ(k,m)を低減（正規化）するための正規化係数ＧNとして利用される。
【００５１】
調整値算定部５４は、ステップＳ31で算定した今回の基礎値Ｑ(k,m)が所定の閾値ＱTHを下回るか否かを判定する（Ｓ35）。基礎値Ｑ(k,m)が閾値ＱTHを下回る場合（Ｓ35：YES）、調整値算定部５４は、正規化係数ＧNを１に初期化する（Ｓ36）。他方、基礎値Ｑ(k,m)が閾値ＱTH以上である場合（Ｓ35：NO）、正規化係数ＧNの初期化（Ｓ36）は実行されない。
【００５２】
調整値算定部５４は、今回の基礎値Ｑ(k,m)を正規化係数ＧNで除算する（Ｓ37）。したがって、図７からも理解される通り、正規化係数ＧNが初期値１に設定された状態（例えば基礎値Ｑ(k,m)が１を上回る以前の状態や基礎値Ｑ(k,m)が閾値ＱTHを下回る数値に低下した直後の状態）では、ステップＳ37の調整の前後で基礎値Ｑ(k,m)は変化しない。他方、正規化係数ＧNが１を上回る数値に設定されている場合、ステップＳ37の調整後の基礎値Ｑ(k,m)は、調整前の基礎値Ｑ(k,m)を下回る数値に調整される。すなわち、図７から理解される通り、基礎値Ｑ(k,m)が最初に１を上回る数値に増加した時点ｔ1から基礎値Ｑ(k,m)が閾値ＱTHを下回る時点ｔ2までの期間では、数式(3)の演算で算定された基礎値Ｑ(k,m)が低減される。調整値算定部５４は、ステップＳ37の実行後の基礎値Ｑ(k,m)について数式(4)の演算を実行することで調整値Ｇ(k,m)を算定する（Ｓ38）。図７では調整値Ｇ(k,m)が鎖線で図示されている。
【００５３】
第３実施形態においても第１実施形態や第２実施形態と同様の効果が実現される。また、第３実施形態では、基礎値Ｑ(k,m)が１を上回る場合に当該基礎値Ｑ(k,m)が低減されるから、音響信号ｘのうち目的成分の直後の残響成分（初期反射音）も高精度に抑圧できるという利点がある。
【００５４】
＜第４実施形態＞
第４実施形態の調整値算定部５４は、以下の数式(7)の演算で調整値Ｇ(k,m)を算定する。
【数７】
2015004959.tif 000009
数式(7)から理解される通り、調整値算定部５４は、基礎値Ｃ(k,m)と基礎値Ｑ(k,m)との加重和を調整値Ｇ(k,m)として算定する。基礎値Ｃ(k,m)は、以下の数式(8)で表現される通り、指標値ＲA(k,m)に対する指標値ＲB(k,m)の相対比であり、前掲の対比例の調整値Ｇ(k,m)に相当する。基礎値Ｑ(k,m)の算定には、例えば第１実施形態の数式(2)や第２実施形態の数式(3)が適用される。
【数８】
2015004959.tif 000010
基礎値Ｃ(k,m)および基礎値Ｑ(k,m)の各々の算定には共通の指標値ＲA(k,m)および指標値ＲB(k,m)が適用される。ただし、基礎値Ｃ(k,m)と基礎値Ｑ(k,m)とで指標値ＲA(k,m)の平滑化係数αAや指標値ＲB(k,m)の平滑化係数αBを相違させることも可能である。
【００５５】
基礎値Ｃ(k,m)の加重値ｗ1と基礎値Ｑ(k,m)の加重値ｗ2とは所定値に設定される。加重値ｗ1および加重値ｗ2の設定方法は任意であるが、調整値Ｇ(k,m)の作用で音響信号ｘの残響成分が有効に抑圧されるように加重値ｗ1および加重値ｗ2を算定する以下の方法が好適である。
【００５６】
音響信号ｘに付与された残響効果の振幅-周波数特性は、変調角周波数ωと残響時間ＴRとを変数とする数式(9)の変調伝達関数Ψ(ω)で近似される。なお、数式(9)の変調伝達関数Ψ(ω)については、例えば、M. Unoki, et. al., "An improved method based on the MTF concept for restoring the power envelope from a reverberant signal", Acoustical science and technology 25(4), p. 232-242にも詳述されている。
【数９】
2015004959.tif 000011
【００５７】
数式(9)の変調伝達関数Ψ(ω)で表現される残響効果を音響信号ｘ（周波数成分Ｘ(k,m)）に対する調整値Ｇ(k,m)の乗算で抑圧するためには、調整値Ｇ(k,m)が変調伝達関数Ψ(ω)の逆特性に近似する必要がある。以上の観点から、調整値Ｇ(k,m)と変調伝達関数Ψ(ω)の逆特性との相違（誤差）を評価するための数式(10)の誤差関数Ｅ(w)を導入する。誤差関数Ｅ(w)は、加重値ｗ1および加重値ｗ2の関数である。なお、数式(10)では、変調角周波数ωを変数とする変調伝達関数Ψ(ω)と対比するために、基礎値Ｃ(k,m)および基礎値Ｑ(k,m)にも変調角周波数ωを導入した。
【数１０】
2015004959.tif 000012
数式(10)から理解される通り、誤差関数Ｅ(w)は、調整値Ｇ(k,m)および変調伝達関数Ψ(ω)の積と１との自乗誤差を変調角周波数ωの全範囲にわたり積算する関数である。したがって、調整値Ｇ(ω)の作用が変調伝達関数Ψ(ω)の逆特性（すなわち残響効果の振幅-周波数特性の逆特性）に近似するほど誤差関数Ｅ(w)は小さい数値となる。第４実施形態の調整値算定部５４は、誤差関数Ｅ(w)が最小となるように加重値ｗ1および加重値ｗ2を算定する。具体的には、加重値ｗ1および加重値ｗ2による誤差関数Ｅ(w)の偏微分をゼロと仮定することで導出される２元１次の連立方程式から加重値ｗ1および加重値ｗ2が算定される。
【００５８】
第４実施形態においても第１実施形態と同様の効果が実現される。また、第４実施形態では、基礎値Ｃ(k,m)と基礎値Ｑ(k,m)との加重和が調整値Ｇ(k,m)として算定されるから、基礎値Ｑ(k,m)のみに応じて調整値Ｇ(k,m)を設定する構成と比較して、多様な音響特性の残響成分を高精度に抑圧できるという効果が実現される。また、第４実施形態では特に、調整値Ｇ(k,m)の作用が残響効果の振幅-周波数特性（変調伝達関数Ψ(ω)）の逆特性に近似するように基礎値Ｃ(k,m)の加重値ｗ1と基礎値Ｑ(k,m)の加重値ｗ2とが設定されるから、残響成分の抑圧精度が向上するという効果は格別に顕著である。
【００５９】
なお、誤差関数Ｅ(w)の内容は適宜に変更される。例えば、前掲の数式(10)では、調整値Ｇ(k,m)の作用が変調伝達関数Ψ(ω)の逆特性に近似するという条件を評価するための誤差関数Ｅ(w)を例示したが、加重値ｗ1と加重値ｗ2との総和が１になるという条件（以下「第１拘束条件」という）を追加した以下の数式(11)の誤差関数Ｅ1(w)を利用することも可能である。
【数１１】
2015004959.tif 000013
数式(11)の係数λ1は所定の正数である。加重値ｗ1と加重値ｗ2との総和が１に近付くほど数式(11)の誤差関数Ｅ1(w)の右辺の第２項は小さい数値となる。すなわち、数式(11)の右辺の第２項が第１拘束条件に相当する。調整値算定部５４は、誤差関数Ｅ1(w)が最小となるように加重値ｗ1および加重値ｗ2を算定する。以上の構成では、加重値ｗ1と加重値ｗ2との総和が１になるという第１拘束条件のもとで加重値ｗ1および加重値ｗ2が算定されるから、例えば基礎値Ｃ(k,m)と基礎値Ｑ(k,m)との双方が最大値１に設定されたような場合でも音響信号ｘと音響信号ｙとの音量差を抑制する（入出力の関係を維持する）ことが可能である。すなわち、音響信号ｙの音量が過剰に増大（ひいては発散）する事態や、音響信号ｙの音量が過剰に低減されて再生音を聴取できない事態等が防止され、結果的に音響信号ｙを安定的に生成できるという利点がある。
【００６０】
また、加重値ｗ1および加重値ｗ2の各々の絶対値が抑制されるという条件（以下「第２拘束条件」という）を追加した以下の数式(12)の誤差関数Ｅ2(w)を利用することも可能である。
【数１２】
2015004959.tif 000014
数式(12)の係数λ2は所定の正数である。各加重値ｗnの自乗ｗn2の総和が小さいほど数式(12)の誤差関数Ｅ2(w)の右辺の第２項は小さい数値となる。すなわち、数式(12)の右辺の第２項が第２拘束条件に相当する。調整値算定部５４は、誤差関数Ｅ2(w)が最小となるように加重値ｗ1および加重値ｗ2を算定する。以上の構成によれば、加重値ｗ1および加重値ｗ2の一方が他方と比較して過度に大きい数値となる事態を防止できる。なお、第１拘束条件および第２拘束条件の双方を加味して加重値ｗ1および加重値ｗ2を算定することも可能である。
【００６１】
基礎値Ｑ(k,m)を、Ｎ個の基礎値Ｑ1(k,m)〜ＱN(k,m)の加重和として算定することも可能である。各基礎値Ｑn(k,m)（ｎ＝１〜Ｎ）は、平滑化係数αA(n)を適用した強度Ｐ(k,m)の指数移動平均で算定された指標値ＲA(k,m)_nと、平滑化係数αB(n)（αB(n)＞αA(n)）を適用した強度Ｐ(k,m)の指数移動平均で算定された指標値ＲB(k,m)_nとを数式(2)または数式(3)に適用することで算定される。平滑化係数αA(n)および平滑化係数αB(n)の一方または双方は基礎値Ｑn(k,m)毎に相違する。数式(7)の基礎値Ｃ(k,m)についても同様に、Ｎ個の基礎値Ｃ1(k,m)〜ＣN(k,m)の加重和として算定され得る。各基礎値Ｑn(k,m)の加重値や各基礎値Ｃn(k,m)の加重値は、前掲の数式(10)の誤差関数Ｅ(w)と同様に、調整値Ｇ(k,m)の作用が残響効果の振幅-周波数特性（変調伝達関数Ψ(ω)）の逆特性に近似するように選定される。
【００６２】
＜第５実施形態＞
以上の各形態ではモノラルの音響信号ｘを想定した。第５実施形態の音響信号ｘは、左チャネルの音響信号ｘLと右チャネルの音響信号ｘRとで構成されるステレオ信号である。解析処理部４２は、音響信号ｘLと音響信号ｘRとを時間領域または周波数領域で加算した音響信号の各周波数成分Ｘ(k,m)を前述の各形態と同様に解析することで各周波数の調整値Ｇ(k,m)を単位期間毎に算定する。周波数分析部３２は、音響信号ｘLの各周波数成分ＸL(k,m)と音響信号ｘRの各周波数成分ＸR(k,m)とを単位期間毎に算定する。残響調整部３４は、周波数成分ＸL(k,m)に調整値Ｇ(k,m)を作用させて周波数成分ＹL(k,m)を生成するとともに周波数成分ＸR(k,m)に調整値Ｇ(k,m)を作用させて周波数成分ＹR(k,m)を生成する。波形生成部３６は、周波数成分ＹL(k,m)から時間領域の音響信号ｙLを生成するとともに周波数成分ＹR(k,m)から時間領域の音響信号ｙRを生成する。すなわち、音響信号ｘLおよび音響信号ｘRの残響成分を抑圧したステレオ形式の音響信号ｙ（ｙL，ｙR）が生成される。第５実施形態においても第１実施形態と同様の効果が実現される。
【００６３】
なお、複数のチャネルで構成される音響信号ｘの残響成分を抑圧する具体的な構成および方法は以上の例示に限定されない。例えば、周波数成分ＸL(k,m)および周波数成分ＸR(k,m)の最大値を周波数成分Ｘ(k,m)として選択したうえで調整値Ｇ(k,m)の算定等を実行することも可能である。また、前掲の数式(1A)および数式(1B)に代えて、以下に例示する各演算で指標値Ｒi(k,m)（ｉ＝Ａ,Ｂ）を算定することも可能である。
【００６４】
以下の数式(13)で表現される通り、音響信号ｘLの周波数成分ＸL(k,m)と音響信号ｘRの周波数成分ＸR(k,m)との相互相関の指数移動平均を指標値Ｒi(k,m)として算定することが可能である。数式(13)の記号＊は複素共役を意味する。
【数１３】
2015004959.tif 000015
【００６５】
以下の数式(14)で表現される通り、音響信号ｘLの各周波数の強度ＰL(k,m)（例えばパワー|ＸL(k,m)|2）と音響信号ｘRの各周波数の強度ＰR(k,m)（例えばパワー|ＸR(k,m)|2）との平均値（単純平均）の指数移動平均を指標値Ｒi(k,m)として算定することが可能である。
【数１４】
2015004959.tif 000016
【００６６】
以下の数式(15)で表現される通り、音響信号ｘLの周波数成分ＸL(k,m)の自己相関と音響信号ｘRの周波数成分ＸR(k,m)の自己相関との乗算値の指数移動平均を指標値Ｒi(k,m)として算定することが可能である。
【数１５】
2015004959.tif 000017
【００６７】
以下の数式(16)や数式(17)を指標値Ｒi(k,m)の算定に利用することも可能である。
【数１６】
2015004959.tif 000018
【数１７】
2015004959.tif 000019
【００６８】
以上の例示から理解される通り、指標値ＲA(k,m)および指標値ＲB(k,m)は、音響信号ｘの強度に追従して時間的に変動する数値として包括的に表現され、音響信号ｘの振幅|Ｘ(k,m)|の任意乗の次元（単位）で表現される物理量の移動平均が典型例として例示される。
【００６９】
なお、以上の説明では、左右２チャネルの音響信号ｘを例示したが、３チャネル以上の音響信号ｘを処理することも可能である。例えば、５チャネル（Ｃ,Ｌ,Ｒ,Ｌs,Ｒs）の音響信号ｘを想定する。中央チャネル（Ｃ）の音響信号ｘCから解析処理部４２が算定した調整値Ｇ(k,m)を残響調整部３４は音響信号ｘCに作用させる。左前方チャネル（Ｌ）の音響信号ｘLと右前方チャネル（Ｒ）の音響信号ｘRとに応じたモノラルの音響信号から解析処理部４２が算定した調整値Ｇ(k,m)を残響調整部３４は音響信号ｘLと音響信号ｘRとに作用させる。また、左後方チャネル（Ｌs）の音響信号ｘLsと右後方チャネル（Ｒs）の音響信号ｘRsとに応じたモノラルの音響信号から解析処理部４２が算定した調整値Ｇ(k,m)を残響調整部３４は音響信号ｘLsと音響信号ｘRsとに作用させる。なお、各チャネルの組合わせは任意である。例えば左前方チャネル（Ｌ）の音響信号ｘLと左後方チャネル（Ｌs）の音響信号ｘLsとから調整値Ｇ(k,m)を算定する構成や、右前方チャネル（Ｒ）の音響信号ｘRと右後方チャネル（Ｒs）の音響信号ｘRsとから調整値Ｇ(k,m)を算定する構成も採用される。
【００７０】
＜第６実施形態＞
第６実施形態の音響処理装置１００は、第１実施形態の解析処理部４２を図８の解析処理部７０に置換した構成である。解析処理部７０は、第１実施形態の解析処理部４２と同様に、各周波数の調整値Ｇ(k,m)を単位期間毎に算定する要素であり、図８に例示される通り、残響推定部７２と調整値算定部７４とを含んで構成される。
【００７１】
残響推定部７２は、音響信号ｘ（各周波数成分Ｘ(k,m)）の残響成分Ｖr(k,m)を推定する。具体的には、残響推定部７２は、音響信号ｘの後期残響成分を近似する以下の数式(18)の演算で、各周波数の残響成分Ｖr(k,m)を単位期間毎に算定する。なお、残響成分を近似する数式(18)については、K. Lebart, et. al., "A New Method Based on Spectral Subtraction for Speech Dereberberation", Acta Acustica, Vol. 87 (2001) p.359-366にも詳述されている。
【数１８】
2015004959.tif 000020
数式(18)の記号Ｔdは、単位期間の個数（したがって自然数）で表現された遅延量を意味する。遅延量Ｔdは、例えば利用者からの指示に応じて可変に設定される。例えば、遅延量Ｔdは、利用者が任意に指定した時間Ｔを各単位期間の時間差（フレームシフト）Ｔsで除算した結果を整数化した数値（Ｔd＝round(Ｔ/Ｔs)）に設定される。例えば音響信号ｘのサンプリング周波数が８ｋHzや１６kHzである場合、遅延量Ｔdは単位期間の２個分から３個分程度の数値に設定される。また、数式(18)の時間Ｔxは、例えば利用者が指定した時間Ｔに応じて設定される。具体的には、時間Ｔに応じた遅延量Ｔdと各単位期間の時間差Ｔsとの積が数式(18)の時間Ｔxとして好適である（Ｔx＝Ｔd・Ｔs）。記号Δは、以下の数式(19)で表現される通り、残響時間ＴRに応じた変数である。
【数１９】
2015004959.tif 000021
【００７２】
また、数式(18)の記号Ｖx(k,m)は、音響信号ｘのパワー|Ｘ(k,m)|2の移動平均（以下の例示では指数移動平均）であり、以下の数式(20)で定義される。数式(20)の記号β1は、指数移動平均の平滑化係数であり、例えば０.９程度の数値（１未満の正数）に設定される。
【数２０】
2015004959.tif 000022
【００７３】
第６実施形態の調整値算定部７４は、以下の数式(21)の演算で各周波数の調整値Ｇ(k,m)を単位期間毎に算定する。
【数２１】
2015004959.tif 000023
数式(21)は、ウィナーフィルタ（Wiener Filter）の理論を利用したフィルタ（Quasi- Wiener Filter）であり、音響信号ｘに作用することで残響成分Ｖr(k,m)を抑圧する。以上の説明から理解される通り、第６実施形態の調整値算定部７４は、音響信号ｘの強度（以上の例示では振幅|Ｘ(k,m)|）と音響信号ｘから推定される残響成分Ｖr(k,m)とを分母に包含するとともに音響信号ｘの強度を分子に包含する演算式（数式(21)）を利用して調整値Ｇ(k,m)を算定する。調整値算定部７４が数式(21)の演算で生成した調整値Ｇ(k,m)を、残響調整部３４は、第１実施形態と同様に音響信号ｘ（各周波数成分Ｘ(k,m)）に作用させる。
【００７４】
数式(21)の残響成分Ｖr(k,m)は、数式(18)および数式(19)から理解される通り、残響時間ＴRの関数である。具体的には、残響時間ＴRが長いほど残響成分Ｖr(k,m)は長時間にわたり大きい数値に維持される。すなわち、数式(21)の残響成分Ｖr(k,m)は、第１実施形態の指標値ＲB(k,m)と同様に、残響時間ＴRが長い場合に調整値Ｇ(k,m)を減少させるように作用する。したがって、第６実施形態によれば、第１実施形態と同様に、残響時間が長い場合でも音響信号ｘの残響成分を高精度に抑圧できるという利点がある。
【００７５】
なお、第６実施形態の調整値Ｇ(k,m)は、事後ＳＮ比σpos（σpos＝(|Ｘ(k,m)|2／Ｖr(k,m))−１）を利用して以下の数式(22)のように変形され得る。なお、事後ＳＮ比σposの利用については、数式(18)の説明で参照した文献（Lebart）にも開示されている。
【数２２】
2015004959.tif 000024
数式(22)の事後ＳＮ比σposは、以下の数式(23)で表現される通り、事前ＳＮ比σpriで近似的に表現される。
【数２３】
2015004959.tif 000025
数式(23)の記号β2は、１未満の正数に設定される。数式(23)の記号max[0, ]は、括弧内の数値を非負値に制限する演算を意味する。
【００７６】
以上の近似を数式(22)に適用すると、事前ＳＮ比σpriに応じて調整値Ｇ(k,m)を算定するための以下の数式(24)が導出される。したがって、調整値算定部７４は、数式(21)に代えて数式(24)の演算で調整値Ｇ(k,m)を算定することも可能である。
【数２４】
2015004959.tif 000026
【００７７】
なお、以上の説明では、ウィナーフィルタの理論を残響成分Ｖr(k,m)の抑圧に適用した場合を例示したが、事前ＳＮ比σpriを利用したウィナーフィルタの理論を、MMSE-STSA法やMMSE対数振幅法等の公知の雑音抑圧技術に拡張することで、音響信号ｘから残響成分Ｖr(k,m)を抑圧することも可能である。
【００７８】
＜第７実施形態＞
第２実施形態の係数設定部４４は、調整値Ｇ(k,m)の算定に適用する加重係数ηを、残響時間ＴRに応じて可変に設定する構成としたが、第７実施形態の係数設定部４４は、指標値ＲA(k,m)および指標値ＲB(k,m)に応じて加重係数ηを可変に設定する。なお、第７実施形態においても第２実施形態と同様にモノラルの音響信号ｘを想定した。
【００７９】
図９は、第７実施形態における音響処理装置１００の構成図である。図９に例示される通り、第７実施形態の音響処理装置１００は、第２実施形態と同様の要素を包含する。第７実施形態の解析処理部４２は、第２実施形態と同様に、指標算定部５２と調整値算定部５４とを含んで構成される。指標算定部５２は、数式(1A)の演算により指標値ＲA(k,m)を算定するとともに数式(1B)の演算により指標値ＲB(k,m)を算定する。調整値算定部５４は、係数設定部４４が設定した加重係数ηを適用した数式(3)および数式(4)の演算で調整値Ｇ(k,m)を算定する。
【００８０】
第７実施形態の係数設定部４４は、指標算定部５２が算定した指標値ＲA(k,m)および指標値ＲB(k,m)を数式(25)に適用して変数ｇ(k,m)を算定し、変数ｇ(k,m)に応じた加重係数ηを設定する。変数ｇ(k,m)は、指標値ＲA(k,m)に対する指標値ＲB(k,m)の相対比であって、前述の対比例における調整値Ｇ(k,m)に相当する。
【数２５】
2015004959.tif 000027
【００８１】
図１０は、変数ｇ(k,m)と加重係数ηとの関係を例示するグラフである。係数設定部４４は変数ｇ(k,m)を関数Ｗ(ｇ(k,m))に適用することで加重係数ηを設定する（η＝Ｗ(ｇ(k,m))）。関数Ｗ(ｇ(k,m))は、図１０に例示される非線形関数である。図１０から理解される通り、加重係数ηは、変数ｇ(k,m)が０以上かつ１未満の場合は変数ｇ(k,m)に連動して増加（単調増加）する一方、変数ｇ(k,m)が負数または１以上である場合には０に設定される。図１０に示されるように、変数ｇ(k,m)が０に近接するとき（音響信号ｘの残響時間が比較的短い場合）、加重係数ηは小さい数値に設定されるから、調整値Ｇ(k,m)は大きい数値に設定される。他方、変数ｇ(k,m)が１に近接するとき（音響信号ｘの残響時間が比較的に長い場合）、加重係数ηは大きい数値（所定の正数ａに近い数値）に設定されるから、調整値Ｇ(k,m)は小さい数値に設定される。
【００８２】
図１１は、加重係数ηを適用して算定された調整値Ｇ(k,m)の時間変化と残響時間との関係を例示するグラフである。前述の対比例における調整値Ｇ(k,m)（Ｇ(k,m)＝ＲB(k,m)／ＲA(k,m)）の時間変化が鎖線で併記されている。図１１の部分(A)は、音響信号ｘの残響時間が比較的に短い場合を例示し、図１１の部分(B)は、音響信号ｘの残響時間が比較的に長い場合を例示する。前述の各形態と同様に、図１１における区間ＳAは、直接音および初期反射音が存在する区間に相当し、区間ＳBは後期残響成分が存在する区間に相当する。
【００８３】
音響信号ｘの残響時間が比較的に長い場合には、指標値ＲB(k,m)は長時間にわたり大きい数値に維持されて、変数ｇ(k,m)は１に近い数値に維持される。加重係数ηは、変数ｇ(k,m)の値が１に近接するときには、比較的大きい数値に設定される。その結果、調整値Ｇ(k,m)は長時間にわたって小さい数値に設定される。図１１(B)を参照して調整値Ｇ(k,m)と対比例とを比較すると、時点ｔｘ以降では、調整値Ｇ(k,m)は対比例に対して充分に低減されていることが判る。以上の構成によれば、音響信号ｘの残響時間が比較的に長い場合でも、調整値Ｇ(k,m)を有効に抑制することができる。
【００８４】
以上に説明した通り、第７実施形態では、音響信号ｘの残響時間が長い場合には、係数算定部４４によって加重係数ηが大きい数値に設定されるから、第２実施形態と同様に、音響信号ｘの残響成分を高精度に抑圧するように調整値Ｇ(変数ｇ(k,m))を設定することが可能である。第７実施形態の構成によれば、第２実施形態と比較して、残響時間ＴRの推定が不要であることから残響時間ＴRの誤差に影響されずに残響成分を抑圧できるという効果を奏する。
【００８５】
＜変形例＞
以上に例示した各形態は多様に変形され得る。具体的な変形の態様を以下に例示する。以下の例示から任意に選択された２以上の態様は適宜に併合され得る。
【００８６】
（１）前述の各形態における周波数分析部３２と残響調整部３４と波形生成部３６とは、音響信号ｘを時間領域で処理する図１２の要素（周波数分析部６２，残響調整部６４，波形生成部６６）に置換され得る。なお、解析処理部４２の構成および動作は第１実施形態と同様である。
【００８７】
周波数分析部６２は、帯域分割部６２２と包絡抽出部６２４とを含んで構成される。帯域分割部６２２は、信号供給装置１２から供給される音響信号ｘを、相異なる周波数帯域に対応する複数の帯域成分ｘ(k)（ｘ(1),ｘ(2),……）に時間領域で分解する。例えば、通過帯域が相違する複数の帯域通過フィルタで構成されるフィルタバンクが帯域分割部６２２として利用される。包絡抽出部６２４は、複数の帯域成分ｘ(k)の各々を包絡成分ｘE(k)と残余成分ｘR(k)とに分解する。包絡成分ｘE(k)は、帯域成分ｘ(k)の時間波形の包絡線に相当する成分であり、残余成分ｘR(k)は、帯域成分ｘ(k)から包絡成分ｘE(k)を除外した成分である。包絡成分ｘE(k)の抽出には、例えばヒルベルト変換等の公知の信号処理技術が任意に採用される。
【００８８】
残響調整部６４は、周波数分析部６２（包絡抽出部６２４）が生成した各包絡成分ｘE(k)に、解析処理部４２が生成した調整値Ｇ(k,m)を作用させる。具体的には、残響調整部６４は、包絡成分ｘE(k)に調整値Ｇ(k,m)を乗算することで包絡成分ｙE(k)を生成する。包絡成分ｘE(k)の時間軸上の各サンプルには、解析処理部４２が当該サンプルの時点について算定した調整値Ｇ(k,m)が乗算される。
【００８９】
波形生成部６６は、第１合成部６６２と第２合成部６６４とを含んで構成される。第１合成部６６２は、残響調整部６４が各周波数帯域について生成した包絡成分ｙE(k)と、当該周波数帯域の残余成分ｘR(k)とを合成（例えば乗算や加算）することで帯域成分ｙ(k)を生成する。以上の説明から理解される通り、帯域成分ｙ(k)は、帯域成分ｘ(k)から残響成分を抑圧した音響成分である。第２合成部６６４は、第１合成部６６２が生成した複数の帯域成分ｙ(k)を合成（例えば加算）することで音響信号ｙを生成する。図１２の構成でも、前述の形態と同様の効果が実現される。
【００９０】
（２）第１実施形態では、前掲の数式(2)の基礎値Ｑ(k,m)を調整値Ｇ(k,m)として利用したが、数式(2)の基礎値を利用した以下の数式(26)の演算で調整値Ｇ(k,m)を算定することも可能である。
【数２６】
2015004959.tif 000028
数式(26)の係数νは数値１を上回る所定値に設定される。数式(26)から理解される通り、調整値Ｇ(k,m)は、数値１を上限値として係数νと基礎値Ｑ(k,m)との乗算値に設定される。数式(2)の基礎値Ｑ(k,m)を調整値Ｇ(k,m)として利用する構成では、調整値Ｇ(k,m)が１には殆ど到達しないから、音響信号ｙの強度が音響信号ｘと比較して不足する可能性がある。他方、数式(26)で調整値Ｇ(k,m)を算定する構成では、調整値Ｇ(k,m)が１に設定される頻度が充分に確保されるから、音響信号ｘと音響信号ｙとの強度差を低減することが可能である。
【００９１】
（３）音響信号ｘの目的成分の発音開始から相当の時間が経過して強度Ｐ(k,m)が充分に減衰した区間（以下「後期減衰区間」という）では、指標値ＲA(k,m)または指標値ＲB(k,m)が変動することで基礎値Ｑ(k,m)（調整値Ｇ(k,m)）が不安定となり、音響信号ｙの音量が揺動して再生音の音質が低下する可能性がある。以上の傾向を考慮すると、後期減衰区間にて音響信号ｙの音量の変動を抑制する構成が好適である。
【００９２】
具体的には以下の数式(27)で表現される通り、調整値算定部５４は、指標値ＲA(k,m)が所定の閾値ＲTHを下回る場合（すなわち後期減衰区間内である場合）に、調整値Ｇ(k,m)を、基礎値Ｑ(k,m)とは無関係の所定値Ｇminに固定し、指標値ＲA(k,m)が閾値ＲTHを上回る場合に基礎値Ｑ(k,m)を調整値Ｇ(k,m)として採択する。以上の構成では、指標値ＲA(k,m)が閾値ＲTHを下回る後期減衰区間内において調整値Ｇ(k,m)が所定値Ｇminに維持されるから、後期減衰区間内で調整値Ｇ(k,m)の変動に起因した音響信号ｙの音量の変動を抑制することが可能である。
【数２７】
2015004959.tif 000029
【００９３】
なお、数式(27)では、指標値ＲA(k,m)が閾値ＲTHを下回る場合に調整値Ｇ(k,m)を所定値Ｇminに設定したが、指標値ＲA(k,m)が閾値ＲTHを下回るという条件に加えて、基礎値Ｑ(k,m)が所定の閾値を上回る数値に変動した場合（すなわち調整値Ｇ(k,m)が不安定に変動し得る場合）に調整値Ｇ(k,m)を所定値Ｇminに設定することも可能である。すなわち、指標値ＲA(k,m)が閾値ＲTHを下回る場合でも、基礎値Ｑ(k,m)が安定している場合には基礎値Ｑ(k,m)が調整値Ｇ(k,m)として採択される。また、以上の説明では、後期減衰区間の内外の判定に指標値ＲA(k,m)を利用したが、指標値ＲA(k,m)を上回る追従性で音響信号ｘの強度Ｐ(k,m)に追従する指標値（例えば強度Ｐ(k,m)の移動平均）が閾値ＲTHを下回るか否かに応じて後期減衰区間の内外を判定することも可能である。
【００９４】
（４）前述の各形態では、音響信号ｘのパワー|Ｘ(k,m)|2を強度Ｐ(k,m)として平滑化することで指標値ＲA(k,m)および指標値ＲB(k,m)を算定したが、指標算定部５２による平滑化の対象はパワー|Ｘ(k,m)|2に限定されない。例えば、音響信号ｘの振幅|Ｘ(k,m)|や振幅|Ｘ(k,m)|の４乗|Ｘ(k,m)|4を音響信号ｘの強度Ｐ(k,m)として指標値ＲA(k,m)および指標値ＲB(k,m)を算定することも可能である。すなわち、指標値Ｒi(k,m)は、任意の冪指数ρを適用した以下の数式(28)でも表現され得る。また、第６実施形態で例示した数式(20)についても同様に、音響信号ｘの振幅|Ｘ(k,m)|の冪指数（音響信号ｘの強度の次元）は任意である。
【数２８】
2015004959.tif 000030
【００９５】
（５）前述の各形態では、音響信号ｘの強度Ｐ(k,m)の指数移動平均を指標値ＲA(k,m)および指標値ＲB(k,m)として算定したが、指標値ＲA(k,m)および指標値ＲB(k,m)の算定方法は適宜に変更される。例えば、音響信号ｘの強度Ｐ(k,m)の単純移動平均（あるいは加重移動平均）を指標値ＲA(k,m)および指標値ＲB(k,m)として算定することも可能である。
【００９６】
具体的には、指標算定部５２は、ＭA個の単位期間にわたる強度Ｐ(k,m)の単純移動平均を指標値ＲA(k,m)として算定し、ＭB個の単位期間にわたる強度Ｐ(k,m)の単純移動平均を指標値ＲB(k,m)として算定する。平均個数ＭAは、前述の指数移動平均に適用される平滑化係数αAに対応し、平均個数ＭBは、指数移動平均に適用される平滑化係数αBに対応する。すなわち、平均個数ＭAが平均個数ＭBを上回る数値に設定されることで、前述の形態と同様に、指標値ＲA(k,m)は、指標値ＲB(k,m)と比較して低い追従性で音響信号ｘの強度Ｐ(k,m)の変動に追従する。以上の説明から理解される通り、指数移動平均に適用される平滑化係数（αA，αB）に加えて単純移動平均の平均個数（ＭA，ＭB）も、移動平均に適用される移動平均係数の概念に包含される。
【００９７】
（６）前述の各形態では、基礎値Ｑ(k,m)の演算式（数式(2)）の分母および分子の双方に共通の指標値ＲB(k,m)が包含される場合を例示したが、以下の数式(29)で基礎値Ｑ(k,m)（調整値Ｇ(k,m)）を算定することも可能である。数式(3)も同様に変形され得る。
【数２９】
2015004959.tif 000031
数式(29)の指標値ＲC(k,m)は、例えば、平滑化係数αCを適用した音響信号ｘの強度Ｐ(k,m)の指数移動平均である。すなわち、調整値算定部５４は、指標値ＲA(k,m)と指標値ＲB(k,m)とを分母に包含するとともに指標値ＲC(k,m)を分子に包含する演算式を利用して調整値Ｇ(k,m)を算定する。指標値ＲB(k,m)および指標値ＲC(k,m)の各々は、指標値ＲA(k,m)を上回る追従性で音響信号ｘの強度Ｐ(k,m)に追従する（例えばαB＞αA，αC＞αA）。すなわち、指標値ＲA(k,m)の平滑化の時定数τAは、指標値ＲB(k,m)の平滑化の時定数τBおよび指標値ＲC(k,m)の平滑化の時定数τCを上回る。時定数τBと時定数τCとの大小（平滑化係数αBと平滑化係数αCとの大小）は不問である。
【００９８】
前述の第１実施形態から第５実施形態は、数式(29)の指標値ＲB(k,m)と指標値ＲC(k,m)とを共通させた構成に相当する（ＲB(k,m)＝ＲC(k,m)）。他方、指標値ＲB(k,m)と指標値ＲC(k,m)とを個別に算定する（平滑化係数αBと平滑化係数αCとを相違させる）ことも可能である。前述の各形態のように指標値ＲB(k,m)と指標値ＲC(k,m)とを共通させた構成によれば、指標値ＲB(k,m)と指標値ＲC(k,m)とを個別に算定する構成と比較して指標算定部５２の処理負荷が軽減されるという利点がある。なお、第７実施形態において、ＲB(k,m)／ＲA(k,m)（第２指標値／第１指標値）の代わりに、ＲC(k,m)／ＲA(k,m)（第３指標値／第１指標値）を変数ｇ(k,m)として当該変数ｇ(k,m)に応じて加重係数ηを算定することも可能である。
【００９９】
（７）前述の各形態では、音響信号ｘの強度Ｐ(k,m)の移動平均を指標値ＲB(k,m)として算定したが、音響信号ｘの強度Ｐ(k,m)を指標値ＲB(k,m)として利用することも可能である。すなわち、指標値ＲB(k,m)の算定について強度Ｐ(k,m)の移動平均は省略され得る。同様に、数式(29)の指標値ＲC(k,m)として音響信号ｘの強度Ｐ(k,m)を利用することも可能である。以上の説明から理解される通り、指標値ＲB(k,m)および指標値ＲC(k,m)は、指標値ＲA(k,m)を上回る追従性で音響信号ｘの強度Ｐ(k,m)の時間変化に追従する数値として包括的に表現される。
【０１００】
（８）基礎値Ｑ(k,m)（調整値Ｇ(k,m)）の算定方法は前掲の数式(2)や数式(3)に限定されない。任意の冪指数κを適用して数式(2)や数式(3)の右辺を冪乗することで基礎値Ｑ(k,m)を算定することも可能である。具体的には、数式(2)を変形（一般化）した以下の数式(30A)や数式(3)を変形した以下の数式(30B)が利用され得る。
【数３０】
2015004959.tif 000032
【０１０１】
（９）前述の各形態では、音響信号ｘの残響成分を抑圧する調整値Ｇ(k,m)を例示したが、音響信号ｘの残響成分を強調（抽出）する場合にも本発明は適用される。例えば、基礎値Ｑ(k,m)を所定値（例えば１）から減算することで調整値Ｇ(k,m)（Ｇ(k,m)＝１−Ｑ(k,m)）を算定する構成では、残響成分を強調した音響信号ｙを生成することが可能である。以上の説明から理解される通り、調整値算定部５４は、音響信号ｘの残響成分を調整（抑圧または強調）するための調整値を算定する要素として包括される。
【０１０２】
（１０）第２実施形態および第７実施形態では、加重係数ηを分母に包含する数式(3)の基礎値Ｑ(k,m)を利用した数式(4)の演算で調整値Ｇ(k,m)を算定したが、加重係数ηを適用した調整値Ｇ(k,m)の算定方法は以上の例示に限定されない。例えば、以下の数式(31)に例示される通り、分母に加えて分子にも加重係数ηを含む演算により基礎値Ｑ(k,m)を算定することも可能である。数式(31)の基礎値Ｑ(k,m)を適用した調整値Ｇ(k,m)の算定には前掲の数式(4)が第２実施形態と同様に利用される。加重係数ηの算定には、第２実施形態（数式(5)，(6-1)〜(6-3)）や第７実施形態（η＝Ｗ(g(k,m))）と同様の方法が適用される。以上の構成においても、残響時間が長い場合に残響成分を高精度に抑圧できるという前述の作用が強調される。
【数３１】
2015004959.tif 000033
【０１０３】
（１１）音響信号ｘの残響成分を抑圧する前述の各形態を、公知の音声強調技術やエコー除去技術と併用することも可能である。例えば、音声強調技術やエコー除去技術を利用した処理後の音響信号ｘを対象として、前述の各形態に例示した残響成分の抑圧が実行される。残響成分の抑圧後の音響信号ｙは、例えば音声認識や話者認識等の各種の音響処理の対象としても好適である。
【０１０４】
（１２）携帯電話機等の端末装置と通信するサーバ装置（典型的にはウェブサーバ）で音響処理装置１００を実現することも可能である。例えば、音響処理装置１００は、端末装置から受信した音響信号ｘから音響信号ｙを生成して端末装置に送信する。なお、音響信号ｘの各周波数成分Ｘ(k,m)が端末装置から送信される構成（例えば端末装置が周波数分析部３２を具備する構成）では音響処理装置１００から周波数分析部３２が省略され、残響成分の調整後の各周波数成分Ｙ(k,m)を音響処理装置１００から端末装置に送信する構成（例えば端末装置が波形生成部３６を具備する構成）では音響処理装置１００から波形生成部３６が省略される。また、端末装置が残響調整部３４を具備する構成では、音響処理装置１００から残響調整部３４が省略され、解析処理部４２が生成した調整値Ｇ(k,m)が音響処理装置１００から端末装置に提供される。以上の説明から理解される通り、本発明は、音響信号ｘの残響成分を調整するための調整値Ｇ(k,m)を算定する装置としても実現され得る。
【０１０５】
（１３）音響空間内での反射や散乱に起因した狭義の残響成分に加えて、例えば楽器の演奏音等の響き成分（共鳴成分）も残響成分に含意される。具体的には、ピアノ等の鍵盤楽器の響板による共鳴成分やバイオリン等の弦楽器の共鳴成分（胴鳴り，箱鳴り）の調整にも本発明を適用することが可能である。すなわち、本発明の残響成分は、経時的に減衰する成分（減衰成分）を意味する。
【符号の説明】
【０１０６】
１００……音響処理装置、１２……信号供給装置、１４……放音装置、２２……演算処理装置、２４……記憶装置、３２，６２……周波数分析部、３４，６４……残響調整部、３６，６６……波形生成部、４２，７０……解析処理部、４４……係数設定部、５２……指標算定部、５４，７４……調整値算定部、７２……残響推定部。
【図１】
2015004959.tif 000034
【図２】
2015004959.tif 000035
【図３】
2015004959.tif 000036
【図４】
2015004959.tif 000037
【図５】
2015004959.tif 000038
【図６】
2015004959.tif 000039
【図７】
2015004959.tif 000040
【図８】
2015004959.tif 000041
【図９】
2015004959.tif 000042
【図１０】
2015004959.tif 000043
【図１１】
2015004959.tif 000044
【図１２】
2015004959.tif 000045