JP6191238B2 - 音響処理装置および音響処理方法 - Google Patents

Description

本発明は、音響信号を処理する技術に関し、特に、音響信号に含まれる残響成分を調整する技術に関する。

音響信号に含まれる残響成分を抑圧する技術が従来から提案されている。例えば特許文献１には、音響信号に含まれる残響成分を推定する予測フィルタ係数の確率モデルを利用することで残響成分の予測フィルタ係数を推定し、推定後の予測フィルタを利用して残響成分を抑圧する技術が開示されている。また、非特許文献１には、発音源から収音点までの伝達関数の逆フィルタを推定し、推定後の逆フィルタを音響信号に適用することで残響成分を抑圧する技術が開示されている。

特開２００９−２１２５９９号公報

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

しかし、特許文献１の予測フィルタ係数や非特許文献１の逆フィルタの高精度な推定には膨大な演算が必要であるという問題がある。以上の事情を考慮して、本発明は、簡便な処理で音響信号の残響成分を調整（抑圧または強調）することを目的とする。

以上の課題を解決するために、本発明の音響処理装置は、残響時間に応じた第１移動平均係数（例えば平滑化係数α₁，係数ζ₁（ζ_1n），平均個数Ｍ₁）を算定する係数算定手段と、係数算定手段が算定した第１移動平均係数を適用した音響信号の強度の移動平均で第１指標値（例えば指標値Ｒ₁(k,m)）を算定する指標算定手段と、音響信号の残響成分を調整するための調整値を、指標算定手段が算定した第１指標値と、当該第１指標値を上回る追従性で音響信号の強度に追従する第２指標値（例えば指標値Ｒ₂(k,m)）とに応じて算定する調整値算定手段とを具備する。例えば、調整値算定手段は、第１指標値を分母に包含するとともに第２指標値を分子に包含する演算により調整値を算定する。以上の構成では、音響信号の強度の移動平均で調整値が算定されるから、簡便な処理で音響信号の残響成分を調整することが可能である。また、残響時間に応じて第１移動平均係数が可変に算定されるから、音響信号の残響成分を高精度に抑圧可能な調整値を算定できるという利点もある。

本発明の好適な態様において、係数算定手段は、残響時間と、音響信号のうち強度が算定される各単位期間の時間差（例えば時間差ＴS）とに応じて、第１移動平均係数を算定する。以上の態様では、残響時間と各単位期間の時間差とに応じて第１移動平均係数が算定されるから、各単位期間の時間差が変更された場合でも、音響信号の残響成分を高精度に抑圧可能な調整値を算定できる利点がある。

本発明の好適な態様において、指標算定手段は、第１移動平均係数を適用した音響信号の強度の移動平均で第１指標値を算定するとともに、第１移動平均係数とは相違する第２移動平均係数（例えば平滑化係数α₂，係数ζ₂（ζ_2n），平均個数Ｍ₂）を適用した音響信号の強度の移動平均で第２指標値を算定し、係数算定手段は、第１移動平均係数に対応する第１調整成分（例えば調整成分Ｇ₁(ω)）と第２移動平均係数に対応する第２調整成分（例えば調整成分Ｇ₂(ω)）とを含む調整値のうちの第１調整成分が、音響信号に付与された残響効果の振幅-周波数特性の逆特性に近似するように、第１移動平均係数を算定する。以上の構成では、調整値の第１調整成分が残響効果の振幅-周波数特性の逆特性に近似するように第１移動平均係数が算定されるから、音響信号の残響成分を高精度に抑圧可能な調整値を算定できるという効果は格別に顕著である。

ただし、第１調整成分が残響成分の振幅-周波数特性の逆特性に近似するという条件だけでは、変調周波数の高域側の音響成分が過剰に強調される可能性がある。以上の事情を考慮すると、第２調整成分が、第１調整成分と音響信号に付与された残響効果の振幅-周波数特性との積の逆特性に近似するように、第２移動平均係数を算定する構成が好適である。以上の構成によれば、第１調整成分の振幅-周波数特性に起因した高域側の過剰な強調が抑制されるから、音響信号の残響成分を高精度に抑圧可能な調整値を算定できるという効果は格別に顕著である。なお、残響成分の振幅-周波数特性は、例えば変調伝達関数（例えば変調伝達関数Ψ(ω)）で近似的に表現される。

以上の各態様に係る音響処理装置は、音響信号の処理に専用されるＤＳＰ（Digital Signal Processor）などのハードウェア（電子回路）によって実現されるほか、ＣＰＵ（Central Processing Unit）等の汎用の演算処理装置とプログラムとの協働によっても実現される。コンピュータが読取可能な記録媒体に格納された形態で提供されてコンピュータにインストールされ得る。記録媒体は、例えば非一過性（non-transitory）の記録媒体であり、ＣＤ-ＲＯＭ等の光学式記録媒体（光ディスク）が好例であるが、半導体記録媒体や磁気記録媒体等の公知の任意の形式の記録媒体を包含し得る。また、例えば、本発明のプログラムは、通信網を介した配信の形態で提供されてコンピュータにインストールされ得る。

本発明の好適な態様に係る音響処理装置の構成図である。 解析処理部の構成図である。 各指標値と調整値との関係の説明図である。 変調伝達関数のグラフである。 調整成分Ｇ₁(ω)の振幅-周波数特性である。 調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積の振幅-周波数特性である。 調整成分Ｇ₂(ω)の振幅-周波数特性である。 調整値Ｇ(ω)と変調伝達関数Ψ(ω)との積の振幅-周波数特性である。 演算処理装置の動作のフローチャートである。 変形例における音響処理装置の部分的な構成図である。

図１は、本発明の好適な実施形態に係る音響処理装置１００の構成図である。図１に示すように、音響処理装置１００には信号供給装置１２と放音装置１４とが接続される。信号供給装置１２は、音響信号ｘを音響処理装置１００に供給する。音響信号ｘは、発音源から放射された音響に対して音響空間内で反射または散乱した残響成分（初期反射成分および後期残響成分）を付加した音響の時間波形を示すサンプル系列である。例えば、収録音や合成音等の既存の音響に対して事後的に残響効果を付与した音響信号ｘや、残響効果がある音響空間（例えば音響ホール等）内で実際に収録された音響の音響信号ｘが好適に利用される。周囲の音響を収音して音響信号ｘを生成する収音装置や、可搬型または内蔵型の記録媒体から音響信号ｘを取得して音響処理装置１００に供給する再生装置や、通信網から音響信号ｘを受信して音響処理装置１００に供給する通信装置が信号供給装置１２として採用され得る。

音響処理装置１００は、音響信号ｘの残響成分（後期残響成分）を抑圧した時間領域の音響信号ｙを生成する残響抑圧装置である。放音装置１４（例えばスピーカやヘッドホン）は、音響処理装置１００が生成した音響信号ｙに応じた音波を再生する。なお、音響信号ｙをデジタルからアナログに変換するＤ/Ａ変換器や音響信号ｙを増幅する増幅器等の図示は便宜的に省略した。

図１に示すように、音響処理装置１００は、演算処理装置２２と記憶装置２４とを具備するコンピュータシステムで実現される。記憶装置２４は、演算処理装置２２が実行するプログラムや演算処理装置２２が使用する各種のデータを記憶する。半導体記録媒体や磁気記録媒体等の公知の記録媒体や複数種の記録媒体の組合せが記憶装置２４として任意に採用され得る。音響信号ｘを記憶装置２４に記憶した構成（したがって信号供給装置１２は省略される）も好適である。

演算処理装置２２は、記憶装置２４に記憶されたプログラムを実行することで、音響信号ｘから音響信号ｙを生成するための複数の機能（周波数分析部３２，残響調整部３４，波形生成部３６，変数設定部４２，解析処理部４４）を実現する。なお、演算処理装置２２の各機能を複数の装置に分散した構成や、専用の電子回路（例えばＤＳＰ）が演算処理装置２２の一部の機能を実現する構成も採用され得る。

周波数分析部３２は、音響信号ｘを構成する複数の周波数成分（周波数スペクトル）Ｘ(k,m)を時間軸上の単位期間（フレーム）毎に順次に生成する。記号ｋは、周波数軸上に離散的に設定された複数の周波数（帯域）のうち任意の１個の周波数を指定する変数であり、記号ｍは、時間軸上の任意の１個の単位期間（時間軸上の特定の時点）を指定する変数である。時間軸上で相前後する各単位期間の時間差（フレームシフト）ＴSは単位期間の時間長を下回る。したがって、相前後する各単位期間は時間軸上で相互に重複する。各周波数成分Ｘ(k,m)の算定には、短時間フーリエ変換等の公知の周波数解析が任意に採用され得る。

変数設定部４２は、音響信号ｘの処理に適用される各変数の数値を設定する。具体的には、本実施形態の変数設定部４２は、周波数分析部３２による各周波数成分Ｘ(k,m)の生成（短時間フーリエ変換）に適用される各単位期間の時間差ＴSと、音響信号ｘの残響成分に想定される残響時間（例えば音響信号ｘが発音および収録された音響空間の残響時間）ＴRとを可変に設定する。具体的には、変数設定部４２は、例えば入力装置（図示略）に対する利用者からの指示に応じて時間差ＴSを設定する。

また、変数設定部４２は、例えば入力装置に対する利用者からの指示（残響時間ＴRの指定）に応じて残響時間ＴRを設定する。なお、音響空間内の音響の収音で音響信号ｘを生成する収音装置を信号供給装置１２として採用した構成では、所定の測定用信号を音響空間内に放音したときに収録される音響信号ｘを解析することで変数設定部４２が残響時間ＴRを実測することも可能である。残響時間ＴRの測定には、例えばインパルス積分法（Schroeder法）やマルチステップ線形予測等の公知の残響解析技術が任意に採用され得る。また、周波数が時間的に連続に変化する時間伸長信号（TSP：Time Stretched Pulse）が測定用信号として好適である。

解析処理部４４は、音響信号ｘの各周波数成分Ｘ(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)を単位期間毎に順次に算定する。本実施形態の指標算定部５４は、音響信号ｘの各周波数成分Ｘ(k,m)の強度（パワー）|Ｘ(k,m)|²の時系列を平滑化することで各周波数の指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)を算定する。具体的には、以下の数式(1A)および数式(1B)で表現される通り、音響信号ｘの強度|Ｘ(k,m)|²の指数移動平均が指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)として算定される。

記号α_i（ｉ＝１,２）は、指数移動平均の平滑化係数（すなわち、強度|Ｘ(k,m)|²の平滑化の時定数を規定する係数）であり、１未満の正数に設定される。平滑化係数α_iは、過去の指標値Ｒ_i(k,m-1)に対する最新（現在）の強度|Ｘ(k,m)|²の加重値に相当する。

図２の係数算定部５２は、平滑化係数α₁および平滑化係数α₂を設定する。平滑化係数α₂は、平滑化係数α₁を上回る数値に設定される（α₂＞α₁）。したがって、指標値Ｒ₁(k,m)における強度|Ｘ(k,m)|²の平滑化の時定数τ₁は、指標値Ｒ₂(k,m)における強度|Ｘ(k,m)|²の平滑化の時定数τ₂を上回る（τ₁＞τ₂）。係数算定部５２が各平滑化係数α_iを算定する具体的な動作については後述する。指標算定部５４は、係数算定部５２が設定した各平滑化係数α_iを適用した音響信号ｘの強度|Ｘ(k,m)|²の指数移動平均を指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)として算定する。

図３の部分(B)には、音響信号ｘの１個の周波数の周波数成分Ｘ(k,m)から算定される指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)の時間変化が例示されている。図３の部分(A)のように強度|Ｘ(k,m)|²が指数減衰する室内インパルス応答（ＲＩＲ：Room Impulse Response）を音響信号ｘとして音響処理装置１００に供給した場合の指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)が図３の部分(B)には図示されている。

図３の部分(B)から理解される通り、指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)は、音響信号ｘの強度|Ｘ(k,m)|²に追従して経時的に変化する。ただし、前述の通り、指標値Ｒ₁(k,m)の時定数τ1は指標値Ｒ₂(k,m)の時定数τ2を上回るから、指標値Ｒ₁(k,m)は、指標値Ｒ₂(k,m)と比較して低い追従性で音響信号ｘの強度|Ｘ(k,m)|²の変動に追従する。具体的には、図３の部分(B)から把握される通り、室内インパルス応答の開始の時点ｔ0の直後の区間では、指標値Ｒ₂(k,m)が指標値Ｒ₁(k,m)を上回る変化率で急峻に増加する。そして、指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)は、時間軸上の相異なる時点で極大値に到達し、指標値Ｒ₂(k,m)は指標値Ｒ₁(k,m)を上回る変化率で減少する。

以上のように指標値Ｒ₁(k,m)と指標値Ｒ₂(k,m)とは相異なる変化率で変化するから、指標値Ｒ₁(k,m)と指標値Ｒ₂(k,m)との大小は時間軸上の特定の時点ｔxで反転する。すなわち、時点ｔ0から時点ｔxまでの区間ＳAでは指標値Ｒ₂(k,m)が指標値Ｒ₁(k,m)を上回り、時点ｔx以降の区間ＳBでは指標値Ｒ₁(k,m)が指標値Ｒ₂(k,m)を上回る。区間ＳAは、室内インパルス応答のうち直接音および初期反射成分が存在する区間に相当し、区間ＳBは、室内インパルス応答のうち後期残響成分が存在する区間に相当する。

図２の調整値算定部５６は、指標算定部５４が算定した指標値Ｒ₁(k,m)と指標値Ｒ₂(k,m)とに応じた調整値Ｇ(k,m)を各周波数について単位期間毎に算定する。具体的には、調整値算定部５６は、以下の数式(2)で表現される通り、指標値Ｒ₁(k,m)に対する指標値Ｒ₂(k,m)の相対比を調整値Ｇ(k,m)として算定する。ただし、指標値Ｒ₁(k,m)に対する指標値Ｒ₂(k,m)の相対比が数値１を上回る場合には、調整値Ｇ(k,m)は数値１（調整値Ｇ(k,m)の上限値）に設定される。なお、調整値Ｇ(k,m)の上限値は任意であり、例えば数値１を下回る所定の正数（例えば０.９）にも設定され得る。

指標値Ｒ₁(k,m)と指標値Ｒ₂(k,m)とが図３の部分(B)のように変化する場合の調整値Ｇ(k,m)の時間変化が図３の部分(C)に例示されている。図３の部分(C)から理解される通り、指標値Ｒ₂(k,m)が指標値Ｒ₁(k,m)を上回る区間ＳA（直接音および初期反射成分が存在する区間）では調整値Ｇ(k,m)は最大値１に設定され、指標値Ｒ₁(k,m)が指標値Ｒ₂(k,m)を上回る区間ＳB（後期残響成分が存在する区間）では調整値Ｇ(k,m)は経時的に減少する。したがって、指標算定部５４が算定した調整値Ｇ(k,m)を残響調整部３４が音響信号ｘに作用させることで、音響信号ｘの残響成分を抑圧した音響信号ｙが生成される。以上に説明した通り、本実施形態では、音響信号ｘの強度|Ｘ(k,m)|²に相異なる時定数で追従する指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)に応じて調整値Ｇ(k,m)が算定されるから、特許文献１や非特許文献１の技術と比較して簡便に音響信号ｘの残響成分を抑圧できるという利点がある。

係数算定部５２による平滑化係数α_iの設定について以下に詳述する。以下の説明では、数式(3)で表現される通り、数式(1A)および数式(1B)の係数｛１−α_i｝を便宜的に係数（忘却係数）ζ_iに置換する。

数式(3)をＺ変換すると以下の数式(4)が導出される。

数式(4)から以下の数式(5)が導出される。

数式(5)で表現される指標値Ｒ₁[z]および指標値Ｒ₂[z]を前掲の数式(2)に適用することで、調整値Ｇ[z]を表現する以下の数式(6)が導出される。

数式(6)から理解される通り、調整値Ｇ[z]は、係数ζ₁に対応する数式(7A)の調整成分Ｇ₁[z]と、係数ζ₂に対応する数式(7B)の調整成分Ｇ₂[z]とに分解される（Ｇ[z]＝Ｇ₁[z]・Ｇ₂[z]）。

ところで、音響信号ｘに付与された残響効果の振幅-周波数特性は、変調角周波数ωと残響時間ＴRとを変数とする数式(8)の変調伝達関数Ψ(ω)で近似される。なお、数式(8)の変調伝達関数Ψ(ω)については、例えば、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にも詳述されている。

図４は、残響時間ＴRを相違させた複数の場合（ＴR＝０.１,０.３,０.５,１,２［sec］）について変調周波数ｆ（ω＝２πｆ）と変調伝達関数Ψ(ω)との関係を併記したグラフである。図４から理解される通り、変調周波数ｆ（変調角周波数ω）の高域側ほど変調伝達関数Ψ(ω)の数値は減少し、残響時間ＴRが長いほど高域側での変調伝達関数Ψ(ω)の減少が顕著である、という概略的な傾向がある。

以上に説明した変調伝達関数Ψ(ω)と調整値Ｇ[z]との対比を検討する観点から、数式(7A)および数式(7B)に変調角周波数ωを導入すると、以下の数式(9A)および数式(9B)が導出される。数式(9A)および数式(9B)の記号ｊは虚数単位を意味し、記号ＴSは、前述の通り、時間軸上で相前後する各単位期間の時間差（フレームシフト）に相当する。

本実施形態の係数算定部５２は、調整値Ｇ(k,m)の作用で音響信号ｘの残響成分が有効に抑圧されるように平滑化係数α₁（係数ζ₁）および平滑化係数α₂（係数ζ₂）を算定する。具体的には、調整値Ｇ(k,m)の作用が残響成分の振幅-周波数特性の逆特性（変調伝達関数Ψ(ω)の逆特性）に近似するように平滑化係数α₁および平滑化係数α₂が算定される。係数ζ₁および係数ζ₂の具体的な算定方法を以下に詳述する。

＜平滑化係数α₁（係数ζ₁）の算定＞
数式(9A)の調整成分Ｇ₁(ω)の振幅-周波数特性は、以下の数式(10)で表現される。

数式(10)の右辺の余弦項（cosωＴS）をTaylor展開することで、以下の数式(11)が導出される。

数式(8)の変調伝達関数Ψ(ω)で表現される残響効果を音響信号ｘ（周波数成分Ｘ(k,m)）に対する調整成分Ｇ₁(ω)の作用で抑圧するためには、調整成分Ｇ₁(ω)の振幅-周波数特性が変調伝達関数Ψ(ω)の逆特性に近似する必要がある（|Ｇ(ω)|・Ψ(ω)＝１）。数式(11)の総和項（Σ）を便宜的に無視したうえで、数式(8)の変調伝達関数Ψ(ω)と数式(11)の調整成分Ｇ₁(ω)との類似性に着目すると、以下の数式(12)が成立する場合に、調整成分Ｇ₁(ω)の振幅-周波数特性が変調伝達関数Ψ(ω)の逆特性に近似すると理解できる。

数式(12)は、係数ζ₁の２次方程式である。数式(12)に２次方程式の解の公式を適用すると、残響時間ＴRと時間差ＴSとに応じて係数ζ₁を算定する以下の数式(13)が導出される。なお、数式(13)の導出（２次方程式の求解）では、係数ζ₁の値域（ζ₁＜１）を考慮して解の符号を採択した。

本実施形態の係数算定部５２は、変数設定部４２が設定した時間差ＴSおよび残響時間ＴRを数式(13)に適用することで時間差ＴSおよび残響時間ＴRに応じた係数ζ₁を算定し、当該係数ζ₁に応じた平滑化係数α₁（α₁＝１−ζ₁）を算定する。

＜平滑化係数α₂（係数ζ₂）の算定＞
図５は、数式(13)で算定された係数ζ₁を数式(9A)に適用することで算定される調整成分Ｇ₁(ω)の振幅-周波数特性である。図５および後掲の図６から図８では、残響時間ＴRを相違させた複数の場合（ＴR＝０.３，０.８，１.２）について周波数特性が併記されている。図５から理解される通り、変調周波数ｆ（変調角周波数ω）の高域側ほど調整成分Ｇ₁(ω)の数値は増加し（高域通過特性）、残響時間ＴRが長いほど高域側での調整成分Ｇ₁(ω)の増加が顕著である、という傾向がある。

図６は、図５の調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝の振幅-周波数特性である。図６から理解される通り、変調周波数ｆ（変調角周波数ω）の高域側ほど調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝は増加するという傾向がある。したがって、調整成分Ｇ₁(ω)を調整値Ｇ(k,m)として音響信号ｘの周波数成分Ｘ(k,m)に作用させた場合、音響信号ｙでは、変調周波数ｆの高域側の音響成分が過剰に強調される可能性がある。以上の事情を考慮して、本実施形態では、変調周波数ｆの高域側ほど積｛|Ｇ₁(ω)|Ψ(ω)｝が増加するという以上の傾向が抑制されるように調整成分Ｇ₂(k,m)の係数ζ₂（平滑化係数α₂）を算定する。

具体的には、本実施形態の係数算定部５２は、調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝の逆特性に調整成分Ｇ₂(ω)の振幅-周波数特性が近似するように係数ζ₂（平滑化係数α₂）を算定する。係数ζ₂の算定には公知の技術が任意に採用され得るが、例えば以下に例示される通り、自己回帰（AR：Auto-Regressive）モデルを規定する自己回帰係数の算定に利用されるYule-Walker法が好適である。

具体的には、係数算定部５２は、数式(13)で算定した係数ζ₁を数式(9A)に適用することで算定される調整成分Ｇ₁(ω)と数式(8)の変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝の逆特性（１/｛|Ｇ₁(ω)|Ψ(ω)｝）のパワースペクトルを逆フーリエ変換することで自己相関関数を算定し（Wiener-Khinchinの定理）、自己相関関数を適用したYule-Walker方程式（正規方程式）から自己回帰係数を算定する。Yule-Walker方程式の解法としてはDurbinの再帰法（Levinson-Durbinアルゴリズム）が好適に利用される。

調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝の逆特性は、以下の数式(14)で表現されるIIR（Infinite Impulse Response）フィルタで近似され、以上に説明したYule-Walker法で算定される自己回帰係数は、数式(14)のIIRフィルタの係数ａ₁に相当する。

数式(7B)で表現される調整成分Ｇ₂[z]と数式(14)のIIRフィルタとの類似性に着目すると、数式(14)の係数ａ₁（Yule-Walker法で算定された自己回帰係数）が調整成分Ｇ₂[z]の係数ζ₂に対応する（ａ₁＝−ζ₂）と理解できる。なお、数式(7B)の分子の係数｛１−ζ₂｝は、形式的な係数（ゲイン項）であり振幅-周波数特性の本質的な傾向には関与しないから便宜的に無視することが可能である。また、IIRフィルタの係数ａ₁は負数であるから、数式(7B)の係数ζ₂の値域（ζ₂＞０）とも整合する。

以上の説明から理解される通り、係数算定部５２は、調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝の逆特性に近似する周波数特性のフィルタ（IIRフィルタ）の係数ａ₁に応じて調整成分Ｇ₂(k,m)の係数ζ₂を算定し、当該係数ζ₂に応じた平滑化係数α₂（α₂＝１−ζ₂）を算定する。

図７は、以上の手順で算定された係数ζ₂を数式(9B)に適用することで算定される調整成分Ｇ₂(ω)の振幅-周波数特性である。図７から理解される通り、調整成分Ｇ₂(ω)の振幅-周波数特性には、調整成分Ｇ₁(ω)の振幅-周波数特性とは逆に、変調周波数ｆの高域側ほど調整成分Ｇ₂(ω)の数値が減少し（低域通過特性）、残響時間ＴRが長いほど高域側での調整成分Ｇ₂(ω)の減少が顕著である、という傾向がある。

図８は、以上の手順で算定された係数ζ₁および係数ζ₂を数式(9A)および数式(9B)に適用することで算定される調整値Ｇ(ω)と数式(8)の変調伝達関数Ψ(ω)との積｛|Ｇ(ω)|Ψ(ω)｝の振幅-周波数特性である。図６と図８とを対比すると、係数算定部５２が算定した係数ζ₂に応じた調整成分Ｇ₂(ω)を加味することで、調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)｝が変調周波数ｆの高域側ほど増加する（高域側が過剰に強調される）という図６の傾向が抑制されることが理解できる。

図９は、記憶装置２４に記憶されたプログラムに従って演算処理装置２２が実行する残響調整処理のフローチャートである。入力装置（図示略）に対する利用者からの指示（残響調整処理の開始指示）を契機として残響調整処理が開始される。残響調整処理を開始すると、変数設定部４２は、各単位期間の時間差ＴSと残響時間ＴRとを設定する（Ｓ1）。解析処理部４４の係数算定部５２は、変数設定部４２が設定した時間差ＴSおよび残響時間ＴRに応じて平滑化係数α₁（係数ζ₁）と平滑化係数α₂（係数ζ₂）とを算定する（Ｓ2）。

以上の手順で各変数が設定されると、単位期間毎に以下のステップＳ3からステップＳ7の処理が実行される。まず、周波数分析部３２は、変数設定部４２が設定した時間差ＴSに応じて時間軸上に画定された単位期間について音響信号ｘの各周波数の周波数成分Ｘ(k,m)を算定する（Ｓ3）。そして、解析処理部４４の指標算定部５４は、係数算定部５２が算定した平滑化係数α₁および平滑化係数α₂と周波数分析部３２が算定した周波数成分Ｘ(k,m)とを適用した数式(1A)および数式(1B)の演算で指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)を算定する（Ｓ4）。調整値算定部５６は、指標算定部５４が算定した指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)を適用した数式(2)の演算で各周波数の調整値Ｇ(k,m)を算定する（Ｓ5）。

残響調整部３４は、指標算定部５４がステップＳ5で算定した調整値Ｇ(k,m)を音響信号ｘの各周波数成分Ｘ(k,m)に作用させることで周波数成分Ｙ(k,m)を算定する（Ｓ6）。波形生成部３６は、残響調整部３４が算定した各周波数成分Ｙ(k,m)から単位期間の音響信号ｙを生成する（Ｓ7）。音響信号ｘの全部の単位期間について以上の処理（Ｓ3〜Ｓ7）が反復される（Ｓ8：NO）ことで、音響信号ｘの残響成分を抑圧した音響信号ｙが生成される。

以上に説明した通り、本実施形態では、残響時間ＴRおよび時間差ＴSに応じて平滑化係数α₁（係数ζ₁）および平滑化係数α₂（係数ζ₂）が算定されるから、平滑化係数α₁および平滑化係数α₂を、残響時間ＴRや時間差ＴSに依存しない所定値に固定した構成と比較して、音響信号ｘの残響成分を高精度に調整可能な調整値Ｇ(k,m)を算定できるという利点がある。

また、調整成分Ｇ₁(k,m)の作用が残響成分の振幅-周波数特性（変調伝達関数Ψ(ω)）の逆特性に近似するように係数ζ₁が算定されるから、音響信号ｘの残響成分は有効に調整される。しかも、調整成分Ｇ₁(ω)と残響効果の振幅-周波数特性（変調伝達関数Ψ(ω)）との積｛|Ｇ₁(ω)|Ψ(ω)｝の逆特性に調整成分Ｇ₂(ω)の振幅-周波数特性が近似するように係数ζ₂が算定されるから、調整成分Ｇ₁(ω)と変調伝達関数Ψ(ω)との積｛|Ｇ₁(ω)|Ψ(ω)において変調周波数ｆの高域側が強調される傾向は抑制される。したがって、音響信号ｘの残響成分を調整値Ｇ(k,m)で高精度に調整できるという効果は格別に顕著である。

＜変形例＞
以上に例示した形態は多様に変形され得る。具体的な変形の態様を以下に例示する。以下の例示から任意に選択された２以上の態様は適宜に併合され得る。

（１）平滑化係数α₁（係数ζ₁）および平滑化係数α₂（係数ζ₂）の算定方法は前述の例示に限定されない。例えば、前述の形態では係数ζ₂の算定に自己回帰モデルを利用したが、以下に例示される通り、自己回帰移動平均（ARMA：Auto-Regressive Moving Average）モデルを利用して係数ζ₁と係数ζ₂とを算定することも可能である。前述の形態では、係数ζ₁と係数ζ₂とを段階的に算定したが、自己回帰移動平均モデルを利用した構成によれば、係数ζ₁と係数ζ₂とを一括的に算定することが可能である。

以下の数式(15)の左辺で表現される１次の自己回帰移動平均モデルの伝達関数は、数式(15)の右辺のように変形される。

数式(6)で表現される調整値Ｇ[z]と数式(15)の変形後の自己回帰移動平均モデルとの類似性に着目すると、自己回帰移動平均モデルの各係数（ａ₁，ｂ₀，ｂ₁）と係数ζ₁およびζ₂との間に以下の数式(16)の対応があると理解できる。

なお、数式(16)の記号σは、変数の自由度を補償するための便宜的な係数を意味する。数式(16)から以下の数式(17)が導出される。

係数算定部５２は、残響成分の振幅-周波数特性の逆特性（数式(8)の変調伝達関数Ψ(ω)の逆特性）に近似するように自己回帰移動平均モデルの各係数（ａ₁，ｂ₀，ｂ₁）を算定する。そして、係数算定部５２は、自己回帰移動平均モデルの各係数を数式(17)に適用することで係数ζ₁と係数ζ₂と係数σとを算定し、係数ζ₁に応じた平滑化係数α₁と係数ζ₂に応じた平滑化係数α₂とを算定する。自己回帰移動平均モデルの各係数（ａ₁，ｂ₀，ｂ₁）の算定には、修正型Yule-Walker（MYW：Modified Yule-Walker）法が好適に利用される。以上の説明から理解される通り、自己回帰移動平均モデルを利用することで、係数ζ₁（平滑化係数α₁）と係数ζ₂（平滑化係数α₂）とを一括的に算定することが可能である。

他方、指標算定部５４は、係数算定部５２が算定した平滑化係数α₁および平滑化係数α₂に応じた指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)を算定し、調整値算定部５６は、係数σを含む以下の数式(18)の演算で調整値Ｇ(k,m)を算定する。

なお、以上の説明では１次の自己回帰移動平均モデルを例示したが、同様の構成を、以下に例示する通り、高次の自己回帰移動平均モデルまで拡張ないし一般化することも可能である。以下の数式(19)の左辺で表現されるＮ次の自己回帰移動平均モデルの伝達関数は、数式(19)の右辺のように変形される。

また、調整値Ｇ[z]についても以下の数式(20)のように高次に拡張される。

数式(20)の調整値Ｇ[z]と数式(19)の変形後の自己回帰移動平均モデルとの類似性に着目すると、係数ζ_1nと係数ζ_2nと係数σとを算定するための以下の数式(21)が導出される。

係数算定部５２は、例えば前述の修正型Yule-Walker法を利用して、残響成分の振幅-周波数特性の逆特性に近似するように自己回帰移動平均モデルの各係数（ａ_n，ｂ₀，ｂ_n）を算定し、各係数を数式(21)に適用することで係数ζ_1nと係数ζ_2nと係数σとを算定する。指標算定部５４は、係数算定部５２が算定した係数ζ_1nを以下の数式(22A)に適用することで指標値Ｒ₁(k,m)を算定し、係数ζ_2nを以下の数式(22B)に適用することで指標値Ｒ₂(k,m)を算定する。調整値算定部５６は、係数算定部５２が算定した係数σと指標算定部５４が算定した指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)とを前掲の数式(18)に適用することで調整値Ｇ(k,m)を算定する。

以上の例示から理解される通り、調整値Ｇ(k,m)の作用が残響成分の振幅-周波数特性の逆特性（変調伝達関数Ψ(ω)の逆特性）に近似するように、音響信号ｘの強度|Ｘ(k,m)|²の平滑化に適用される係数（ζ₁，ζ₂，ζ_1n，ζ_2n）を算定する構成が好適であり、残響成分の振幅-周波数特性の逆特性を表現する時系列モデルの種類（FIR型／IIR型）や時系列モデルの各係数の算定方法は任意である。

（２）前述の形態における周波数分析部３２と残響調整部３４と波形生成部３６とは、音響信号ｘを時間領域で処理する図１０の要素（周波数分析部６２，残響調整部６４，波形生成部６６）に置換され得る。なお、解析処理部４４の構成および動作は前述の形態と同様である。変数設定部４２は、前述の形態と同様に残響時間ＴRを設定するほか、音響信号ｘのサンプリング周期を時間差ＴSとして指定する。

周波数分析部６２は、帯域分割部６２２と包絡抽出部６２４とを含んで構成される。帯域分割部６２２は、信号供給装置１２から供給される音響信号ｘを、相異なる周波数帯域に対応する複数の帯域成分ｘ(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)を合成（例えば加算）することで音響信号ｙを生成する。図１０の構成でも、前述の形態と同様の効果が実現される。

（３）前述の形態では、音響信号ｘの強度|Ｘ(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)は、指標値Ｒ₂(k,m)と比較して低い追従性で音響信号ｘの強度|Ｘ(k,m)|²の変動に追従する。以上の説明から理解される通り、指数移動平均に適用される平滑化係数（α₁，α₂）や係数（ζ₁，ζ₂）に加えて単純移動平均の平均個数（Ｍ₁，Ｍ₂）も、移動平均に適用される移動平均係数の概念に包含される。

（５）前述の形態では、音響信号ｘの残響成分を抑圧する調整値Ｇ(k,m)を例示したが、音響信号ｘの残響成分を強調（抽出）する場合にも本発明は適用される。例えば、数式(2)で算定される調整値Ｇ(k,m)を所定値（例えば１）から減算した調整値｛１−Ｇ(k,m)｝を音響信号ｘに作用させれば、残響成分を強調した音響信号ｙを生成することが可能である。以上の説明から理解される通り、調整値算定部５６は、音響信号ｘの残響成分を調整（抑圧または強調）するための調整値を算定する要素として包括される。

（６）指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)に応じて調整値Ｇ(k,m)を算定する方法は前述の例示に限定されない。例えば、指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)を変数とする所定の演算により調整値Ｇ(k,m)を算定する構成も採用される。以上の説明から理解される通り、調整値算定部５６は、音響信号ｘの残響成分を調整（抑圧または強調）するための調整値Ｇ(k,m)を指標値Ｒ₁(k,m)および指標値Ｒ₂(k,m)に応じて算定する要素として包括される。

（７）前述の形態では、変数設定部４２が各単位期間の時間差ＴSと残響時間ＴRとを可変に設定したが、残響時間ＴRおよび時間差ＴSの一方を所定値に固定することも可能である。したがって、係数算定部５２は、残響時間ＴRおよび時間差ＴSの少なくとも一方に応じて係数ζ₁および係数ζ₂を算定する要素として表現される。例えば既知の残響時間ＴRを適用する構成では残響時間ＴRの算定が省略されるから、演算能力が低い情報処理装置（例えば携帯機器）でも音響処理装置１００を実現することが可能である。なお、残響時間ＴRは、音響空間の気温等の要因にも依存するが、音響空間の音響特性に基本的には依存するから、ひとつの音響空間については、１回の演算で算定された残響時間ＴRを複数回にわたり継続的に適用することが可能である。なお、相異なる音響空間について事前に測定された複数の残響時間ＴRのうち音響信号ｘが収録された場所（例えば音響処理装置１００が使用される場所）に対応する残響時間ＴRを選択して平滑化係数α₁の算定に適用することも可能である。

（８）前述の各形態では、音響信号ｘの強度（パワー）|Ｘ(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)が音響処理装置１００から端末装置に提供される。

（１０）音響空間内での反射や散乱に起因した狭義の残響成分に加えて、例えば楽器の演奏音等の響き成分（共鳴成分）も残響成分に含意される。具体的には、ピアノ等の鍵盤楽器の響板による共鳴成分やバイオリン等の弦楽器の共鳴成分（胴鳴り，箱鳴り）の調整にも本発明を適用することが可能である。すなわち、本発明の残響成分は、経時的に減衰する成分（減衰成分）を意味する。

１００……音響処理装置、１２……信号供給装置、１４……放音装置、２２……演算処理装置、２４……記憶装置、３２，６２……周波数分析部、３４，６４……残響調整部、３６，６６……波形生成部、４２……変数設定部、４４……解析処理部、５２……係数算定部、５４……指標算定部、５６……調整値算定部。

Claims (6)

1. 残響時間に応じた第１移動平均係数を算定する係数算定手段と、
   前記係数算定手段が算定した第１移動平均係数を適用した音響信号の強度の移動平均で第１指標値を算定する指標算定手段と、
   前記音響信号の残響成分を調整するための調整値を、前記指標算定手段が算定した第１指標値と、当該第１指標値を上回る追従性で前記音響信号の強度に追従する第２指標値とに応じて算定する調整値算定手段と
   を具備する音響処理装置。
2. 前記係数算定手段は、前記残響時間と、前記音響信号のうち前記強度が算定される各単位期間の時間差とに応じて、前記第１移動平均係数を算定する
   請求項１の音響処理装置。
3. 前記指標算定手段は、前記第１移動平均係数を適用した前記音響信号の強度の移動平均で前記第１指標値を算定するとともに、前記第１移動平均係数とは相違する第２移動平均係数を適用した前記音響信号の強度の移動平均で前記第２指標値を算定し、
   前記係数算定手段は、前記第１移動平均係数に対応する第１調整成分と前記第２移動平均係数に対応する第２調整成分とを含む前記調整値のうちの前記第１調整成分が、前記音響信号に付与された残響効果の振幅-周波数特性の逆特性に近似するように、前記第１移動平均係数を算定する
   請求項１または請求項２の音響処理装置。
4. 前記係数算定手段は、前記第２調整成分が、前記第１調整成分と前記音響信号に付与された残響効果の振幅-周波数特性との積の逆特性に近似するように、前記第２移動平均係数を算定する
   請求項３の音響処理装置。
5. 前記調整値算定手段は、前記第１指標値を分母に包含するとともに前記第２指標値を分子に包含する演算により前記調整値を算定する
   請求項１から請求項４の何れかの音響処理装置。
6. コンピュータが、
   残響時間に応じた第１移動平均係数を算定し、
   前記算定した第１移動平均係数を適用した音響信号の強度の移動平均で第１指標値を算定し、
   前記音響信号の残響成分を調整するための調整値を、前記算定した第１指標値と、当該第１指標値を上回る追従性で前記音響信号の強度に追従する第２指標値とに応じて算定する
   音響処理方法。

Images