WO2009110118A1

WO2009110118A1 - 音信号の分離方法

Info

Publication number: WO2009110118A1
Application number: PCT/JP2008/065287
Authority: WO
Inventors: 茂樹嵯峨山; 順貴小野; 弘和亀岡; 賢一宮本; ジョナトンルルー
Original assignee: 国立大学法人東京大学
Priority date: 2008-03-05
Filing date: 2008-08-27
Publication date: 2009-09-11
Also published as: JP5294300B2; JP2009210888A; US20110058685A1

Abstract

時間周波数領域におけるスペクトル成分の滑らかさの異方性に着目して、音信号から分離信号を取得する。音信号のスペクトログラムを、時間周波数領域でのスペクトル成分の滑らかさに時間周波数平面上で方向性を備えた複数のサブスペクトログラムの和であると仮定し、各サブスペクトログラムの滑らかさの時間周波数平面上での方向に基づいて、時間周波数領域における前記音信号のスペクトル成分を、少なくとも１つのサブスペクトログラムに分配する少なくとも１つの分配係数を取得し、前記分配係数を用いて前記音信号のスペクトル成分から少なくとも１つのサブスペクトログラムを分離する。

Description

音信号の分離方法

本発明は、音信号、典型的には多重音信号、の分離方法に関するものである。本明細書では、本発明が適用される典型的な例として、音楽音響信号から打楽器音を分離、抽出することを中心に説明するが、本発明は、音楽音響信号からの打楽器音の分離に限定されるものではなく、例えば、機械や装置から発生する工業的な音を分離するものでもよい。

音楽検索、自動採譜などの音楽情報処理の問題においては、音楽音響信号から音高、和声、リズムパターン、テンポなど様々な情報を抽出・認識する必要があるが、これらは未だ難しいタスクであり、近年活発な研究がなされている。音楽音響信号は大きく、メロディーや和声に関連する調波成分と、リズムやドラムパートに関連する打楽器成分の２つに大別されるが、これらは全く異なる性質をもつ信号であり、これらが混在していることが音楽音響信号解析の困難さの一つの要因となっている。しかしながら、モノラル録音された音楽音響信号から調波／打楽器各成分を分離することは容易ではなく、従来は楽譜や楽器の情報なしには行うことができない問題であった。これらの分離がうまくできれば、打楽器やノイズなどの非調波成分を含んだ多声音楽信号の楽音分析における前処理、打楽器パートの強調や打楽器パターン変更といった音楽加工など、多くの分野への応用が期待される。

多重音信号から打楽器（非調波成分）を分離、抽出する手法として幾つかの研究が知られている（非特許文献１乃至３）。

非特許文献１は、特定打楽器の周波数特性テンプレートを用いた音源同定、除去に関するものである。非特許文献１では、除去対象楽器のテンプレートが必要となり、フレームワイズな処理を行っている。また、反復推定による音源同定を行っている。

非特許文献２は、楽器情報を用いない、フレームワイズな調波・非調波成分の分離手法に関するものであり、反復推定による分離手法が行われている。

非特許文献３は、学習データを用いて、あらかじめ調波音、打楽器音の周波数特性を学習するものであり、学習した特徴とのマッチングによって、フレームワイズに分離する。
吉井和佳, 後藤真孝, 奥乃博, "実世界の音楽音響信号に対するドラムスの音源同定を利用したドラムイコライズシステムINTER:Dの開発," 第３回情報処理科学技術フォーラム FIT2004. 亀岡弘和, 後藤真孝, 嵯峨山茂樹, "スペクトル制御エンベロープによる混合音中の周期および非周期成分の選択的イコライザ," 情報処理学会研究報告. 2006-MUS-65, pp77-84. M. Helen, T. Virtanen, "Separation of Drums from Polyphonic Music Using Non-negative Matrix Factorization and Support Vector Machine," In proc, 13th EUSIPCO, 2005.

本発明は、従来のフレームワイズな分析手法に対して、時間周波数領域におけるスペクトル成分の滑らかさの異方性に着目して、音信号を分離することを目的とするものである。本発明のより具体的な一つの目的は、楽器や楽譜に関する情報を全く用いずに、音楽音響信号から調波的な楽器音成分と打楽器的な非調波音成分を分離することにある。

本発明が採用した技術手段は、音信号を分離する方法において、前記音信号のスペクトログラムを、時間周波数領域でのスペクトル成分の滑らかさに時間周波数平面上で方向性を備えた複数のサブスペクトログラムの和であると仮定し、各サブスペクトログラムの滑らかさの時間周波数平面上での方向に基づいて、時間周波数領域における前記音信号のスペクトル成分を、少なくとも１つのサブスペクトログラムに分配する少なくとも１つの分配係数を取得し、前記分配係数を用いて前記音信号のスペクトル成分から少なくとも１つのサブスペクトログラムを分離する、音信号の分離方法、である。

本発明は、時間周波数領域における音信号のスペクトログラムのスペクトル成分の滑らかな方向に着目したものである。すなわち、周波数特性だけでなく、時間的な変化の違いを用いる点に特徴を備えている。本発明では、音信号のスペクトログラムを、時間周波数領域でのスペクトル成分の滑らかさに時間周波数平面上で方向性を備えた複数のサブスペクトログラムの和であると仮定する。すなわち、本発明の対象となる音信号のスペクトログラムを時間周波数領域で異方向に滑らかなスペクトル成分を備えた複数のサブスペクトログラムの和であると仮定する。各サブスペクトログラム（同じサブスペクトログラム）に属するスペクトル成分は時間周波数平面上で概ね同じ方向に滑らかであり、）スペクトル成分の滑らかな方向は、各サブスペクトログラム間で異なる。例えば、ある多重音信号のスペクトログラムは、時間周波数平面上で第１の方向に延出するスペクトル成分群からなる第１サブスペクトログラムと、時間周波数平面上で第２の方向に延出するスペクトル成分群からなる第２サブスペクトログラムの和であると仮定する。あるいは、ある多重音信号のスペクトログラムを、時間周波数平面上で第１の方向に延出するスペクトル成分群からなる第１サブスペクトログラムと、時間周波数平面上で第２の方向に延出するスペクトル成分群からなる第２サブスペクトログラムと、時間周波数平面上で第３の方向に延出するスペクトル成分群からなる第３サブスペクトログラムと、の和であると仮定する。

ここで、本発明は、スペクトログラムにおけるスペクトル成分の滑らかな方向の違いに着目した点に特徴を有するものであるが、分離信号を得る処理ステップにおいて、スペクトログラムを実際に画面に表示することを要しない。本発明においては、分析対象となる音信号が時間周波数領域に変換され、スペクトル成分が得られていればよい。時間周波数領域への変換手段は、典型的な例では、短時間フーリエ変換であるが、ウェーブレット変換、定Ｑフィルタバンク分析、その他のフィルタバンク分析でもよい。また、１つの態様では、本発明は、観測信号である音信号を時間周波数領域に変換するステップと、分離された各サブスペクトログラムに対応するスペクトル成分を時間領域に変換するステップを含み得るが、本発明の全ての処理が時間周波数領域で行われる場合もある。実際のスペクトログラムの計算では、短時間周波数分析によって離散的な時間と周波数ごとに成分が得られる。したがって、スペクトログラムにおける各スペクトル成分（時間周波数成分）は、時間bin（フレーム)と周波数binにより特定される時間周波数binである。

本発明において、音信号のスペクトル成分を分配するための分配係数は、典型的には、いわゆる時間周波数マスク（各時間周波数で０～１の値を取る）として設定され、入力スペクトログラムと時間周波数マスクの乗算によって分離が行われる。音信号の各スペクトル成分を、各サブスペクトログラムに対応するスペクトル成分に分配する分配係数は、一つの態様では、音信号のスペクトログラムが２つのサブスペクトログラムからなる場合に、０あるいは１を取るバイナリマスクである。バイナリマスクを用いることで、良好な耳聴こえが得られる可能性がある。もっとも、分配係数は０あるいは１に限定されるものではなく、その他の比率で分配してもよい。分配係数、すなわち、時間周波数マスクは、入力された音信号のスペクトル成分に基づいて設計される。また、音信号のスペクトログラムが３つのサブスペクトログラムからなる場合には、３つの分配係数の合計が１となるように、各分配係数が設計される。

１つの態様では、前記分配係数の取得は、前記音信号の各スペクトル成分について、各サブスペクトログラムの滑らかさの方向に基づいて、各サブスペクトログラムのスペクトル成分らしさのスコアを取得し、各スコアを指標として、分配係数を取得するステップと、からなる。

１つの態様では、前記スコアを取得するステップは、前記音信号のスペクトログラムを、時間周波数平面上で各スペクトル成分がエネルギーに対応する濃度値を備えた画像とみなし、前記音信号のスペクトログラムから各サブスペクトログラムに属するスペクトル成分の特徴を夫々抽出するフィルタを設定し、各スペクトル成分についての、各サブスペクトログラムに対応するフィルタ処理後のフィルタ出力をスコアとする、ものである。

１つの態様では、前記フィルタは、各サブスペクトログラムのスペクトル成分の滑らかな方向の濃度値を平滑化するローパスフィルタである。また、スペクトル成分の滑らかな方向の特徴を抽出するフィルタは、周波数領域のデジタルフィルタに限定されるものではなく、空間フィルタによっても設計し得ることは当業者に理解される。

１つの態様では、前記音信号のスペクトログラムを、２個のサブスペクトログラムの和であると仮定し、前記スコアを比較し、スコアが大きい方の分配係数を１、スコアが小さい方の分配係数を０とする。あるいは、フィルタ出力値の比に応じて、合計が１となるように分配係数を設定してもよい。

一つの態様では、前記複数のサブスペクトログラムは、周波数方向に滑らかなスペクトル成分からなる第１サブスペクトログラムと、時間方向に滑らかなスペクトル成分からなる第２サブスペクトログラムと、からなり、各サブスペクトログラムのスペクトル成分の滑らかな方向の特徴を抽出するフィルタは、実質的に時間方向の平滑化を行うフィルタと、実質的に周波数方向の平滑化を行うフィルタと、からなる。より具体的には、時間方向のみの１次元ローパスフィルタと、周波数方向のみの１次元ローパスフィルタ、あるいは、時間方向の遮断周波数ωt、周波数方向の遮断周波数ωfが大きく異なる２つの２次元ローパスフィルタ（一方はωt>>ωf、他方はωt<<ωf）などを含む。尚、スペクトル成分の方向が、周波数方向、時間方向でなくても、時間周波数領域においてある一定の方向に滑らかであれば、そのようなスペクトル成分を抽出するフィルタの設計が可能であることは当業者に理解され、そのフィルタ出力を指標として分配係数が設定できることも当業者に理解される。

一つの態様では、前記分配係数の取得は、分配係数をパラメータとして各サブスペクトログラムに分配された各スペクトル成分の滑らかさ指標の関数を含む目的関数を設定し、前記目的関数を最適化するパラメータを推定することで取得される。分配された各スペクトル成分の滑らかさ指標は、着目したスペクトル成分と、当該着目スペクトル成分の時間周波数平面上での近傍にある分配されたスペクトル成分とのエネルギーの差に基づいて決定される。着目スペクトル成分の近傍のスペクトル成分は、典型的には、時間周波数平面上で隣接するスペクトル成分であるが、近傍の範囲はこれに限定されるものではない。分配係数、すなわち、時間周波数マスクの設定は、滑らかさのコストをスペクトログラムの微分の関数で設計し、これを最小化する最適化問題として捉えることができる。

１つの態様では、前記滑らかさの指標の関数は、

である。
　ここで、Ｋ：サブスペクトログラムの数、
　　　　　ｉ：周波数方向のインデックス
　　　　　ｊ：時間方向のインデックス
　　　　　ｆ_ｋ（ｘ）：滑らかさをはかるコスト関数
　　　　　ａ_ｍ，ｎ：時間周波数領域でのある点近傍での重み係数
　　　　　ｍ：周波数方向の近傍を表すインデックス
　　　　　ｎ：時間方向の近傍を表すインデックスウ
　　　　　ｇ（ｘ）：滑らかさ指標における、スペクトログラムのレベル圧縮関数
　　　　　Ｑ^（Ｋ） _ｉ，ｊ：サブスペクトログラムのスペクトル成分
である。

１つの態様では、前記目的関数が、前記音信号のスペクトル成分とパラメータとしての分配係数によって分配された各スペクトル成分の和との距離指標の関数を含む。すなわち、滑らかさのコスト＋距離指標から目的関数が設定され、この目的関数を最小化するように分配係数を最適化する。一つの態様では、距離指標は、Iダイバージェンスである。Iダイバージェンスは、解析的な更新式を求めやすいという利点を有している。距離指標としては、パラメータの更新式が解析的に求められるような距離関数であれば、他の距離指標、例えば、ユークリッド距離(２乗誤差)やマハラノビス距離など、でもよい。分布間距離の要件は、２つの分布がどんな値であっても、関数の値が常に非負であり、かつ０のときにのみ両者の分布が完全に一致することである。

１つの態様では、前記音信号のスペクトログラムを、Ｋ個のサブスペクトログラムの和であると仮定した場合に、前記目的関数は、

である。
　ここで、Ｋ：サブスペクトログラムの数、
　　　　　ｉ：周波数方向のインデックス
　　　　　ｊ：時間方向のインデックス
　　　　　Ｄ（Ａ，Ｂ）：関数Ａと関数Ｂの距離指標
　　　　　φ（ｘ）：距離指標における、スペクトログラムのレベル圧縮関数
　　　　　Ｗ_ｉ，ｊ：観測スペクトル成分
　　　　　ｆ_ｋ（ｘ）：滑らかさをはかるコスト関数
　　　　　ａ_ｍ，ｎ：時間周波数領域でのある点近傍での重み係数
　　　　　ｍ：周波数方向の近傍を表すインデックス
　　　　　ｎ：時間方向の近傍を表すインデックスウ
　　　　　ｇ（ｘ）：滑らかさ指標における、スペクトログラムのレベル圧縮関数
　　　　　Ｑ^（Ｋ） _ｉ，ｊ：サブスペクトログラムのスペクトル成分
である。

１つの態様では、前記目的関数において、

である。

１つの態様では、前記目的関数において、

である。これは、後述する第２実施形態に対応する。

１つの態様では、前記目的関数において、

である。これは、後述する第３実施形態に対応する。

本発明の実施形態では、分離の耳聴こえを良くするための工夫を行っている。人間の聴覚では、音量（音響エネルギー）を対数的（0.3乗程度）に捉える。したがって、小さい音量の変化もある程度認識可能であり、少しでもエネルギーが残っていると、分離できてないように感じる。後述する第２の実施形態では、（１）I-Divergenceはエネルギーをやや対数的に扱う、（２）平方根の滑らかさコスト＝エネルギーをやや対数的に扱う、ことでこれに対応している。また、後述する第３の実施形態では、滑らかさはエネルギーを線形に扱っている。具体的には、レベル圧縮によって、あらかじめエネルギーを0.3乗程度に圧縮している。

第２実施形態、第３実施形態におけるリアルタイム調波音・打楽器音分離について説明する。本来は、入力すべての時間周波数成分を用いて分離するものであり、第２実施形態、第３実施形態の手法では反復処理により、時間を要しうる。しかしながら、隣のフレームのみとの滑らかさを定義することで、比較的高速に計算を可能とし、実時間処理を実現している。すなわち、隣接する時間周波数bin同士のエネルギーの滑らかさを最小化するようにして分配係数を算出している。

具体的には、分析領域のシフトをしながらＥＭ的な反復処理を行う。図９Ａに示すように、所定の分析スペクトログラム領域に対して、１フレームの入力があると、分析スペクトログラム領域で反復更新が行われて次の１フレームの出力のための分配係数が決定され、当該分配係数によって分配されたスペクトル成分が出力され、時間領域に変換される。まとめると、前記音信号を初期分析区間で時間周波数領域に変換してスペクトル成分を取得するステップと、前記音信号を１フレーム分だけ時間周波数領域に変換してスペクトル成分を取得し、当該スペクトル成分を前記分析区間に加えるステップと、前記分析区間のスペクトル成分を用いて、パラメータの推定を行うステップと、前記分析区間で最も古い１フレームのスペクトル成分を、推定されたパラメータを用いて分離し、分離されたスペクトル成分を時間周波数領域に変換する。

目的関数におけるパラメータである分配係数を推定するアルゴリズムとしては、一つの好ましい態様ではＥＭアルゴリズムであるが、最急降下法やニュートン法等の他の最適化アルゴリズムを用いてもよい。また、ＥＭアルゴリズムを解くにあたって、補助変数を導入してもよい。

本発明において、サブスペクトログラムの数は、２以上であれば限定されないが、一つの態様では、音信号のスペクトログラムは、２つのサブスペクトログラムからなり、さらに、典型的な一つの態様では、前記複数のサブスペクトログラムは、周波数方向に滑らかなスペクトル成分からなる第１サブスペクトログラムと、時間方向に滑らかなスペクトル成分からなる第２サブスペクトログラムと、からなる。この場合、一つの態様では、前記多重音信号は打楽器音を含む音楽信号であり、前記第１サブスペクトログラムには、打楽器音のスペクトル成分が含まれる。すなわち、第１サブスペクトログラムは非調波音的成分（典型的には、打楽器音）であり、第２サブスペクトログラムは調波音的成分である。また、本発明の対象となる多重音信号のサブスペクトログラムのスペクトル成分の滑らかな方向は、周波数方向や時間方向に限定されるものではなく、時間周波数領域上で一定の方向に滑らかであれば、スペクトル成分の滑らかな方向に基づいて、多重音信号のスペクトログラムを複数のサブサブスペクトログラムに分離することが可能である。

本発明のハードウエア構成としては、パーソナルコンピュータ等のコンピュータ（具体的には、入力装置、表示装置を含む出力装置、ＣＰＵ，記憶装置（ＲＯＭ，ＲＡＭ等）、これらを接続するバス等、を備えている。）から構成することができる。したがって、本発明は、音信号を分離させるために、コンピュータを、請求項１乃至２６いずれかに記載された方法を実行させる、コンピュータプログラムあるいはコンピュータプログラムを記憶したコンピュータ可読媒体としても提供される。

本発明は、音信号を分離する装置としても提供され、前記装置は、各サブスペクトログラムの滑らかさの時間周波数平面上での方向に基づいて、時間周波数領域における前記音信号のスペクトル成分を、少なくとも１つのサブスペクトログラムに分配する少なくとも１つの分配係数を取得する手段と、前記分配係数を用いて前記音信号のスペクトル成分から少なくとも１つのサブスペクトログラムを分離手段と、を備える。典型的には、前記装置は、さらに、多重音信号を時間周波数領域に変換する手段と、分離された各サブスペクトログラムに対応するスペクトル成分を時間領域に変換する手段と、を備えている。また、本発明は、分離された少なくとも１つのサブスペクトログラムのスペクトル成分を強調あるいは抑圧するステップないし手段を備えていてもよい。例えば、本発明を用いて図１１に示すようなイコライザを実現することができる。
　本明細書において、後述の実施形態では多数の数式が使用されているが、数式番号はセクション毎に独立して付してある。

本発明では、時間周波数領域におけるスペクトル成分の滑らかさの異方性を利用することで、学習データや事前情報を用いることなく、多重音信号から分離信号を取得する音響信号から打楽器音を分離することができる。

本発明では、学習データや打楽器テンプレートなどの楽器固有の情報を用いることなく、音響信号から打楽器音を分離することができる。

［Ａ］本実施形態の概要
本実施形態では調波音と打楽器音の混在した音楽信号を分析対象とし、入力信号の短時間周波数解析によって得られるスペクトログラムをＷ（ｘ，ｔ）とする（ｘ:周波数、ｔ：時刻)。本実施形態の問題は、このＷ（ｘ，ｔ）を打楽器的な音程を持たない非調波成分Ｐ（ｘ，ｔ）と音程を持つ楽器のような調波成分Ｈ（ｘ，ｔ）の２つのスペクトログラムに分解することである。このとき満たすべき要件は、任意の時間周波数（ｘ，ｔ）において、

が成り立つことである。

本実施形態では、調波成分・打楽器成分の異方性に着目する。より具体的には、図１で示すようなポピュラー音楽の音響信号のスペクトログラムが、時間周波数領域において、一般的に周波数方向に形成される山脈ないし畝のようなスペクトル成分と、時間方向に形成される山脈ないし畝のようなスペクトル成分とからなることが多い点に着目する。前者は、打楽器のように時間方向には急峻に変化するが周波数方向にはブロード（滑らか）である成分Ｐ（ｘ，ｔ）に、後者は逆に周波数方向には急峻な形状だが時間方向には滑らかな成分Ｈ（ｘ，ｔ）に対応するとみなすことができ、また２成分は時間周波数平面上においてスパース（同じ時間周波数ｂｉｎ上に存在することが少ない）に存在しているとみなせる。

本実施形態では、入力信号のスペクトログラムを、時間周波数マスクによって２つのスペクトログラムに分解する。すなわち、前述したＰ（ｘ，ｔ）とＨ（ｘ，ｔ）のスパース性から、任意の時間周波数において０～１の値をとる時間周波数マスクｍ_Ｐ（ｘ，ｔ），ｍ_Ｈ（ｘ，ｔ）を設計することで、

とＷ（ｘ，ｔ）を分解できると考えられる。これらの分離スペクトログラムは式(1),(2),(3)の性質を満たす。

時間周波数マスクは、サブスペクトログラムを形成するスペクトル成分の滑らかな方向を検出するように設計される。本実施形態では、打楽器成分のスペクトル成分が周波数方向に滑らかであるという特徴、及び、調波成分のスペクトル成分が時間方向に滑らかであるという特徴を用いて、入力信号のスペクトログラムをそれぞれのスペクトル成分に分離する時間周波数マスクが設計される。０～１の値を取る時間周波数マスクは、一つの態様では、０か１の値を取るバイナリマスクである。

上述のように、本実施形態では、音楽音響信号の時間周波数スペクトログラム上において、調波成分は時間方向に滑らか、打楽器成分は周波数方向に滑らかであるという時間周波数スペクトログラム上の性質の違いを積極的に用いることにより、高速にこれらを分離する。具体的には、与えられた時間周波数スペクトログラムを、時間方向に滑らかな成分と周波数方向に滑らかな成分に分解する相補的な時間周波数マスクを設計し、音楽音響信号の時間周波数スペクトログラムに対して時間周波数マスキングを行って、調波成分・打楽器成分を分離する。設計方法として、１）２次元フィルタを用いる手法、２）Divergenceと滑らかさコストをＥＭアルゴリズム的手法で最小化する手法、３）レベル圧縮したスペクトログラムに対し滑らかさコストをＥＭアルゴリズム的手法で最小化する手法、の３つの実施形態について説明する。各実施形態の説明において、数式番号は、各実施形態毎に独自に付与される。

［Ｂ］第１実施形態
第１実施形態では、観測信号の時間周波数平面のスペクトログラムを画像とみなし、調波的な音と打楽器的な音の持つ一般的な性質の違いを利用した２次元フィルタを用いることで、楽器固有の情報なしで音楽信号から打楽器音と調波音を分離する。

［Ｂ－１］２次元フィルタ出力を用いたマスク設計
時間周波数マスクｍ_Ｐ（ｘ，ｔ）、ｍ_Ｈ（ｘ，ｔ）の設計について述べる。Ｗ（ｘ，ｔ）を画像とみなすと、Ｐ（ｘ，ｔ）とＨ（ｘ，ｔ）の特徴、すなわち、周波数方向のエッジ(縦方向のエッジ)と時間方向のエッジ(横方向のエッジ)、を個別に抽出するような２次元フィルタをかけることで、そのフィルタ出力結果の大小から各時間周波数成分がＰ（ｘ，ｔ）に属するかＨ（ｘ，ｔ）に属するかを決定できる。

Ｗ（ｘ，ｔ）の２次元フーリエ変換成分をＷ（バー）（ａ，ｂ）（ａ:周波数方向のフーリエ成分，ｂ:時間方向のフーリエ成分）とすると、Ｐ（ｘ，ｔ）特徴抽出フィルタＦ（バー）_Ｐ（ａ，ｂ）、Ｈ（ｘ，ｔ）特徴抽出フィルタＦ（バー）_Ｈ（ａ，ｂ）を用いることで、

のようにフィルタ出力結果が得られる。この結果から時間周波数マスクｍ_Ｐ（ｘ，ｔ）、ｍ_Ｈ（ｘ，ｔ）は、

と得られる。

［Ｂ－２］特徴抽出２次元フィルタの設計
前節で述べた２次元フィルタに関して満たすべき要件を検討する。出力結果が各時間周波数成分においてＰ（ｘ，ｔ）らしさ、Ｈ（ｘ，ｔ）らしさの指標となるためには、フィルタ出力が非負の実数になることが望ましいが、必ずしも非負でなくてもよい。また入力スペクトログラムとフィルタ出力の時間周波数が対応している必要がある。前者の実現のためには、フィルタが任意の２次元分布の畳み込みＡ（ａ，ｂ）＊Ａ（ａ，ｂ）で表現される形状であればよく、またその形状がａ，ｂ両軸に対して線対称な実数分布になっていれば後者の性質も満たす。

Ｐ（ｘ，ｔ）、Ｈ（ｘ，ｔ）の特徴をそれぞれ抽出する２次元フィルタＦ（バー）_Ｐ（ａ，ｂ）、Ｆ（バー）_Ｈ（ａ，ｂ）としては様々な形状が考えられる。次に述べる実験では、要件を満たす最も簡単なフィルタとして、Ｆ_Ｐ（ａ，ｂ）は周波数方向のみ、Ｆ_Ｈ（ａ，ｂ）は時間方向のみのローパスフィルタ、

として設計し、ｇ（ａ）やｈ（ｂ）の１次元ローパスフィルタの断面形状としては三角窓やgaussianが利用できる。
三角窓型ローパスフィルタは、

gaussian窓型フィルタは、

と書ける。フィルタを通過した成分の2次元逆フーリエ変換により、P₀(x,t)、H₀(x, t)が得られ、出力結果から時間周波数マスクm_P(x, t)、m_H(x, t) が設計できる。

2次元フィルタは要件を満たす最も簡単なフィルタ形状といえる。三角窓は2つの矩形窓の畳み込みで、Gaussian窓は2つのGaussianの畳み込みで表現できるため、上述したようにフィルタ出力が非負である性質を満たす。この2次元フィルタによって、時間方向、周波数方向に滑らかな成分だけが通過できる。よって各時間周波数binにおいて出力結果の2つの非負の値を比較することで、そのbin が調波音成分らしいか打楽器音の成分らしいかの判断ができる。このときフィルタの持つパラメータとして、ローパスフィルタのカットオフ周波数に対応するσ_P、σ_Hがある。この値が小さいほどより滑らかな成分しか通過しないことになる。提案した2次元フィルタのインパルス応答を考えると、時間周波数方向のうちローパスでない一方がデルタ関数になり、もう一方がsinc関数の2乗(三角窓型の場合)やガウシアン(ガウシアン窓型の場合)になる。このことからこのフィルタ処理は、対象の時間周波数binのスペクトログラムにおいて、周波数方向または時間方向の一方に関してのみの荷重平均を取ることに相当する。逆にスペクトログラムの各時間周波数binの周辺で任意の荷重平均を取る処理は、全て上述した正定値フィルタをかけることに相当する。また、フィルタのインパルス応答が時間周波数軸双方に線対称(つまり時間、周波数双方向に関して偶関数)であれば、時間・周波数成分に関する加重平均の偏りがないため、分離スペクトログラムと元のスペクトログラムで時間周波数binの対応のずれは全くないと考えられる。この性質により、フィルタの出力結果からマスク関数を設計することが妥当と考えられる。

［Ｂ－３］評価実験
［Ｂ－３－１］実際の楽曲への適用結果
ポピュラー音楽の楽曲を用いた分離実験を行った。入力信号として、RWC 研究用音楽データベースよりRWC-MDB-P-2001 No.7より抜粋して使用した(16kHz サンプリング)。入力信号のスペクトログラム、提案アルゴリズム(ローパスフィルタの形状はGaussian)による分離結果を、図５の左図に示す。

結果から、Ｐ（ｘ，ｔ）は周波数方向にブロードな成分、Ｈ（ｘ，ｔ）は周波数方向に急峻だが時間方向に滑らかな成分に分離されたことが分かる。分離音を聴くと、スネアドラムなどの打楽器音はＰ（ｘ，ｔ）に分離されたが、バスドラムやハイハットに関しては特にDuration部分がＨ（ｘ，ｔ）に分離されることが確認された。また歌声においてピッチが連続的に変化する部分はＰ（ｘ，ｔ）、Ｈ（ｘ，ｔ）どちらにも分離され得るが、ローパスフィルタのカットオフ周波数を調整することにより、Ｈ（ｘ，ｔ）の方に多く分離することが可能である。

［Ｂ－３－２］MIDIを用いた定量評価実験
次に、提案アルゴリズムの定量評価実験を行なった。RWC研究用音楽データベースよりRWC-MDBP-2001 No.18の前奏部を入力とし、MIDI形式データをパート別に分離し、各パートをWAV形式に変換してその信号の和を入力とした(16kHz サンプリング)。そして第１実施形態の手法によって得た分離結果の信号と各パート信号との相関を計算することで、Ｐ（ｘ，ｔ）とＨ（ｘ，ｔ）に含まれるエネルギー比率を算出した。その結果を表１に示す。表より、ギターやピアノなどのメロディーや伴奏はＨ（ｘ，ｔ）に、スネアドラムやハイハットはＰ（ｘ，ｔ）分離したが、バスドラムがＨ（ｘ，ｔ）に分離される結果を得た。

第１実施形態では、打楽器音や調波音の特徴としてスペクトログラムの周波数、時間方向の連続性を用いるものであり、スネアドラムなどの打楽器音や、音程を持つ楽器音の分離には適していると考えられる。バスドラムやハイハットのように周波数分布に偏りを持ち比較的音長の長い打楽器音や、ピアノの打鍵音やベースの打弦音、ピッチの変化しやすい歌声などの分離については、特徴抽出２次元フィルタの形状の設計によって解決可能であると考えられる。

［Ｃ］第２実施形態
第１実施形態では、楽器や楽譜に関する情報を全く用いずに、音楽信号からの分離手法として、スペクトログラム上で画像処理的な２次元フィルタを用いた高速な直接計算手法を示した。第２実施形態では、スペクトログラムの滑らかさの異方性に基づいたＥＭアルゴリズムによる反復解法を提案し、計算時間や性能の評価を行うまた、このアルゴリズムを応用して実時間で分離するシステムを提案する。

［Ｃ－１］滑らかさコストの導入
スペクトログラムにおける調波的な成分と打楽器的な成分の異方性を利用して、Ｗ（ｘ，ｔ）からＨ（ｘ，ｔ）とＰ（ｘ，ｔ）を推定する問題を議論する。実装上（ｘ，ｔ）は離散的な座標として取得できるため、以下の議論では離散的な時間周波数領域（ｘ_ｉ，ｔ_ｊ）と定義して議論を行なう(Ｉ:周波数bin数、Ｊ:分析フレーム数)。

本実施形態では、スペクトログラムの滑らかさの異方性を、最小化すべきコストとして、隣り合う時間周波数binとのエネルギーの平方根の二乗誤差

のように表現する。平方根を取ることにより、エネルギーを対数的に捉える人間の聴覚特性により近い滑らかさコストの定式化を実現した。

［Ｃ－２］目的関数最小化によるパラメータ反復推定
［Ｃ－２－１］概要
観測スペクトログラムを調波成分・打楽器成分に分配する時間周波数マスクｍ_Ｈ（ｘ_ｉ，ｔ_ｊ），ｍ_Ｐ（ｘ_ｉ，ｔ_ｊ）を導入する。時間周波数マスクｍ_Ｈ（ｘ_ｉ，ｔ_ｊ），ｍ_Ｐ（ｘ_ｉ，ｔ_ｊ）は数２の条件を満たす。

分配されたエネルギー分布ｍ_Ｐ（ｘ_ｉ，ｔ_ｊ）Ｗ（ｘ_ｉ，ｔ_ｊ）、ｍ_Ｈ（ｘ_ｉ，ｔ_ｊ）Ｗ（ｘ_ｉ，ｔ_ｊ）と、Ｐ（ｘ_ｉ，ｔ_ｊ）、Ｈ（ｘ_ｉ，ｔ_ｊ）との近さを表す分布間距離の指標としてI-Divergenceを採用すると、式(1)(2)の滑らかさコストとの和による目的関数

を最小化する問題として定式化できる。

この目的関数から、時間周波数マスクを固定して式(3)を最小化するＨ（ｘ_ｉ，ｔ_ｊ）とＰ（ｘ_ｉ，ｔ_ｊ）の更新と、Ｈ（ｘ，ｔ）, Ｐ（ｘ，ｔ）を固定して式(3)を最小化するようなｍ_Ｐ（ｘ_ｉ，ｔ_ｊ）とｍ_Ｈ（ｘ_ｉ，ｔ_ｊ）の更新を交互に行なうことにより、目的関数(3)の最小化における局所最適解が得られる。以下に、Ｉダイバージェンスを用いた反復解法について詳述する。

［Ｃ－２－２］エネルギー二乗誤差を滑らかさコストとした解法
Ｉダイバージェンスを用いた反復解法について詳述する。以下の説明において、反復解法における数式番号については、説明の都合上、独自に付与する。ここで解きたいのは、入力スペクトログラムＷ（ｘ、ｔ）（ｘ：周波数、ｔ：時間フレーム）から、打楽器成分Ｐ（ｘ、ｔ）と調波成分Ｈ（ｘ、ｔ）に分離する問題である。これを時間周波数マスクｍ_Ｐ（ｘ、ｔ）とｍ_Ｈ（ｘ、ｔ）を用いて、ＥＭ的なアルゴリズムにより反復的に推定する手法を述べる。Ｗ（ｘ、ｔ）とＰ（ｘ、ｔ）＋Ｈ（ｘ、ｔ）の近さを表す分布間距離として、Ｉダインバージェンスを採用する。この距離指標は、二乗誤差よりも対数的にエネルギーを捉えられ、かつ対数の二乗誤差に比べてエネルギーの非常に小さい部分の誤差を無視できるため、人間の聴覚特性との親和性が高い。また、値は非負であり、Ｗ（ｘ、ｔ）＝Ｐ（ｘ、ｔ）＋Ｈ（ｘ、ｔ）のときに０となる。この分布間距離とＨやＰの滑らかさを表す項を足したものを目的関数Ｊ_１とすると、この分離問題は、Ｊ_１を最小化する問題として定式化できる。Ｊｅｎｓｏｎの不等式より、

という関係を満たす。ここで導入したマスク関数ｍ_Ｐ（ｘ、ｔ）、ｍ_Ｈ（ｘ、ｔ）は、

という制約があり、上述の不等式の符号は、

のときのみ成立する。

Ω_ＰとΩ_Ｈは、滑らかさに関する制約であり、

のような、隣の時間周波数成分との二乗誤差として定義できる。

不等式（１）より、マスク関数ｍ_Ｐ（ｘ、ｔ）、ｍ_Ｈ（ｘ、ｔ）を固定して、Ｊ２を最小化するように、Ｐ（ｘ、ｔ）、Ｈ（ｘ、ｔ）を推定し、次にＰ（ｘ、ｔ）とＨ（ｘ、ｔ）を固定して、式（３）からマスクを更新するような、この二つのステップのパラメータ更新を交互に繰り返すことで、目的関数Ｊ_１が単調減少することが保証できる。また、Ｊ_１≧０は明らかなので、局所最適解に収束する。

実際のデータでは、ｘ、ｔは離散的なので、離散系のモデルで更新式を導出する。Ｊ_２をＰ（ｘ_ｉ，ｔ_ｊ）＝Ｐ_ｉ，ｊとＨ（ｘ_ｉ，ｔ_ｊ）＝Ｈ_ｉ，ｊで偏微分すると、

となる。これを０とおいて、それぞれ、Ｐ_ｉ，ｊ、Ｈ_ｉ，ｊについて解くと２次方程式の２解が得られるが、Ｐ_ｉ，ｊ、Ｈ_ｉ，ｊが正であることから、

と求まる。

反復推定アルゴリズムとしては、
１．初期Ｐ（ｘ_ｉ，ｔ_ｊ）、Ｈ（ｘ_ｉ，ｔ_ｊ）を定める。
２．（３）でｍ_Ｐ（ｘ_ｉ，ｔ_ｊ）、ｍ_Ｈ（ｘ_ｉ，ｔ_ｊ）を更新する。
３．（８）（１２）で、Ｐ（ｘ_ｉ，ｔ_ｊ）、Ｈ（ｘ_ｉ，ｔ_ｊ）を逐次更新する。
４．収束したら終了、しなければ２に戻る。

最終的な分離結果は、推定したマスクｍ_Ｐ（ｘ，ｔ）、ｍ_Ｈ（ｘ，ｔ）を用いて、

のようにして得られる。さらに、聴覚的にはバイナリマスクの方が分離の聴こええが良好であることから、推定したマスクから、

のようにバイナリ化を行って分離する。この時、ｑが大きいほどバイナリ化の効果が大きく、ｑ→∞でバイナリマスクをかけたことと一致する。

［Ｃ－２－３］聴覚特性を考慮した滑らかさ制約の導入
前節の滑らかさに関する制約は、エネルギーの小さい部分と大きい部分とを均等考えた滑らかさを定義している。しかし、人間の聴覚はエネルギーを対数的に捉えることが多いため、分離が良好に行われないおそれがある。そこで、制約を、

のように、エネルギーの平方根の二乗誤差として与える。これは音響エネルギーを、より対数的に捉えた上での滑らかさを考えたことになる。今、分析信号のエネルギーが定数倍になったとき、つまり、Ｗ（ｘ、ｔ）、Ｐ（ｘ、ｔ）、Ｈ（ｘ、ｔ）をＡ倍するとき、ＩダイバージェンスはＡ倍になるが、同時に上述のコスト関数もＡ倍になる。このため、音量の異なる曲でもパラメータσ_Ｐ，σ_Ｈを変える必要がない。さらに、聴覚的にも同じく大きさを対数的に考えた分布間距離指標であるＩダイバージェンスとの親和性も高い。

このコストを用いた場合の更新式を考える。目的関数をＰ（ｘ_ｉ，ｔ_ｊ）＝Ｐ_ｉ，ｊで偏微分すると、

となる。これを０とおいて、Ｐ_ｉ，ｊについて解くと、

となる。
同様に、Ｈ_ｉ，ｊについても解け、

となる。

［Ｃ－２－４］二乗誤差項の補助関数法の利用
前述のＩダイバージェンスを用いた解法では、各Ｐ（ｘ、ｔ）、Ｈ（ｘ、ｔ）の更新式のなかに、隣の時間周波数ｂｉｎの値が必要であった。そこで、滑らかさに関数二乗誤差項に補助関数法を適用し、この問題を解決する。

一般的に、

が成り立つ。等号は、

のときのみ成り立つ。これを利用すると、滑らか制約の項は、

と上限関数が作れる。

よって、目的関数は、

と上限関数が作れることになる。

Ｐ（ｘ、ｔ）、Ｈ（ｘ、ｔ）の更新式は、

となり、
同様に、Ｈについても解け、

となる。

ｍ_{Ｐ，ｉ，ｊ}、ｍ_{Ｈ，ｉ，ｊ}、Ｃ_ｉ，ｊ、Ｄ_ｉ，ｊの更新については、

となる。

反復推定アルゴリズムとしては、
１．初期Ｐ（ｘ_ｉ，ｔ_ｊ）、Ｈ（ｘ_ｉ，ｔ_ｊ）を定める。
２．ｍ_Ｐ（ｘ_ｉ，ｔ_ｊ）、ｍ_Ｈ（ｘ_ｉ，ｔ_ｊ）、Ｃ_ｉ，ｊ、Ｄ_ｉ，ｊを更新する。
３．補助関数を固定した上で、Ｐ（ｘ_ｉ，ｔ_ｊ）、Ｈ（ｘ_ｉ，ｔ_ｊ）を逐次更新する。
４．収束したら終了、しなければ２に戻る。

［Ｃ－２－５］滑らかさ制約の分散の更新
前節で定義した滑らかさのコスト関数では、コストの強さを定めるσ_Ｐ，σ_Ｈを定数としていた。ここでは、このパラメータも変数として更新式を求める。

滑らかさに関するペナルティは、正規分布を用いたＰとＨに関する事前分布

の対数をとって、

とおける。これを用いた目的関数においてσ_Ｐとσ_Ｈの更新式を導出する。

同様に、

となる。
このモデルの場合、Ｐ_ｉ，ｊやＨ_ｉ，ｊの更新式は前節と同じになる。

また、この分散を周波数ごとに異なる値を持たせることも考えられる。この場合、ペナルティ項は、

となる。これを用いた目的関数において、σ_Ｐ，ｉとσ_Ｈ，ｉの更新式を導出する。

同様に、

となる。この分散の場合は、Ｐ_ｉ，ｊのみ更新式が変化する。

［Ｃ－２－６］スパース性の導入
前節の滑らかさの制約の他に、Ｐ_ｉ，ｊやＨ_ｉ，ｊの大きさに関してなるべく０を多くするというスパース制約が導入可能である。これは聴覚的に分離の聴こえ方を良くする働きになると考えられる。前述では、後処理としてバイナリ化を行う方法に言及したが、この項の導入により反復推定中にバイナリ化が行われる。

スパース制約としては、ラプラス分布を仮定し、目的関数に

を足すことで実現する。

また、

でも可能であると考えられる。

前者の場合、更新式は、

となり、
同様に、Ｈ_ｉ，ｊについても解け、

と求まる。これを用いて既述の反復更新を行えばよい。

提案した反復推定においては、時間周波数マスクm_P(x, t),m_H(x, t)は0から1の連続的な値をとる連続値マスクとして推定される。しかし、調波音と打楽器音の成分は時間周波数平面上でスパースに存在していること、またバイナリマスクの方が耳で聴いた分離の性能が良いことなどを考えると、推定された連続値マスクをバイナリマスクに近づけることが有効であると考えられる。推定されたm_P(x, t),m_H(x, t)の大小関係から、バイナリマスクを

と設計できる。しかし、完全なバイナリマスクで設計すると、時間周波数方向におけるスペクトルの連続性が悪いため、耳で聴いたときの分離後の音声があまり良くないことも考えられる。そこで、バイナリの強さを表すパラメータγを用いて、

と設計できる。このγが大きいほどバイナリマスクに近付き、理想的にγ→∞のときに完全なバイナリマスクに、逆にγ＝１の場合には、もとの連続値マスクと一致する。

［Ｃ－３］実時間分離システムの実現
上記解法は、入力信号全体の時間周波数領域における反復解法であるため、一般的には実時間分離は難しい。しかし、スペクトログラムの滑らかさを、隣接した時間周波数binのみを用いた微分的なコストとして表現することで、局所的な分析領域でもある程度妥当な解が得られると考えられる。そこで、局所的な分析時間区間を用い、分析区間の移動とパラメータの反復更新(1～数回)を交互に行なうことで、実時間での調波音・打楽器音分離システムを実現した（図９Ａ）。実時間での調波音・打楽器音分離のステップは以下の通りである。
1. 初期分析区間の入力スペクトログラムを計算する。
2. 新たに入力スペクトログラムを1フレーム計算し、分析区間に加える。
3. 分析領域のスペクトログラムを用いて、分離スペクトログラム、時間周波数マスクを1～数回反復更新する。
4. 分析時間区間で最も古いフレームに対して、推定された時間周波数マスクによる分離を行ない、逆フーリエ変換によって分離信号を出力する。
5. 曲が終われば終了。そうでなければStep.2に戻る。

［Ｃ－４］評価実験
［Ｃ－４－１］実際の楽曲への適用
本節ではポピュラー音楽の実演奏信号を用いた定性的実験を述べる。入力信号として、RWC 研究用音楽データベースからRWC-MDB-P-2001 No.7より抜粋して使用した(16kHz サンプリング)。入力信号のスペクトログラム、提案アルゴリズムの分離結果を、図８に示す。

結果から、Ｐ（ｘ，ｔ）、Ｈ（ｘ，ｔ）が着目した性質を満たすように分離されたことが分かる。結果の音声を聴くと、実施例１の手法に比べ良く分離でき、特に調波音は非常にスムーズに聴こえた。しかし、第１実施形態と同様、ハイハットやバスドラムのduration部分がＨ（ｘ，ｔ）に分離されること、歌声のビブラートや子音がＰ（ｘ，ｔ）に分離されやすいことを確認した。

［Ｃ－４－２］パート別の分離に関する定量評価実験
次にパート別信号を用いた定量的な評価実験を行なった。RWC研究用音楽データベースよりRWC-MDBP-2001 No.18の前奏部8.1秒を入力とし、MIDI形式データをパート別に分離し、各パートをWAV形式に変換してその信号の和を入力とした(16kHzサンプリング)。そして第１実施形態、第２実施形態の手法によって得た分離信号と各パート信号との相関を計算することで、Ｐ（ｘ，ｔ）とＨ（ｘ，ｔ）に含まれるエネルギー比率を算出し、計算時間とともに比較した(表２、CPU3.6GHz のマシンで計算)。表２より、第２実施形態の手法は、第１実施形態の手法に比べて計算コストは増大するが、分離性能を大きく改善できることが分かる。しかし、両手法ともバスドラムは調波音側に分離された。

結果より、スペクトログラムの滑らかさの異方性に基づく解法が、第１実施形態による解法と同様の性質をもった分離を、実時間に比べて十分高速にかつより高い性能で実現したと言える。楽器の知識を用いずに簡便な特徴に基づいた解法のため、比較的音長の長いバスドラムやハイハットの打楽器音、ピアノの打鍵音、ピッチの変化しやすい歌声などは着目した特徴を満たしにくく、楽器分類の通念とは必ずしも対応しない可能性があるが、実時間演算で分離できるメリットは非常に大きいと考えられる。

［Ｄ］第３実施形態
第２実施形態では、Ｗ（ｘ，ｔ）からＨ（ｘ，ｔ）とＰ（ｘ，ｔ）を推定する問題を議論するものであったが、第３実施形態では、Ｈ（ｘ，ｔ）とＰ（ｘ，ｔ）を利用せずに、分配されたスペクトログラムの滑らかさコストを最小化する問題として議論する。

［Ｄ－１］調波音成分／打楽器音成分の事前モデル
Ｆ_ｈ，ｉモノラル音響信号ｆ（ｔ）の短時間フーリエ変換（ＳＴＦＴ）とすると、
Ｆ_ｈ，ｉ＝φ（｜Ｆ_ｈ，ｉ｜^２)となり、ここで、ｈ、ｉは、周波数ｂｉｎ、時間ｂｉｎのインデックスである。Ｆ_ｈ，ｉは、φ（Ａ）＝Ａの時には通常のスペクトログラムを表し、φ（Ａ）＝Ａ^γ（γ＜１）のような凸関数φ（Ａ）を設定することで、レンジ圧縮されたスペクトログラムが生成される。

スペクトログラムの調波成分は通常一定のピッチを備え、滑らかな時間エンベロプを備えた並行状の畝を形成する。これに対して、打音のエネルギーは短時間に集中しており、広域スペクトルエンベロープを備えた縦方向のリッジを形成する。したがって、音楽信号のスペクトログラムでは、典型的に、垂直構造と水平構造が現れる（図１）。また、水平状の複数の畝と垂直状の複数の畝との交差は極めて少ないと考えられる。したがって、ここでの目的は、以下のような適切な時間周波数バイナリマスｍ_ｈ，ｉを見つけることである。

ここで、Ｈ_ｈ，ｉ、Ｐ_ｈ，ｉはそれぞれ、スペクトログラムの調波成分、非調和（打楽器）成分を表す。マスクｍ_ｈ，ｉを設計する一つの手法は、ある事前分布に基づく最大事後推定（ＭＡＰ）推定を適用することである。水平方向、垂直方向にそれぞれ滑らかなＨ_ｈ，ｉ、Ｐ_ｈ，ｉのエンベロープに着目して、各成分について次の事前確率を仮定する。

ベクトルＨ、Ｐは、それぞれ、Ｈ_ｈ，ｉ、Ｐ_ｈ，ｉの集合を表し、σ^２ _Ｈ、σ^２ _Ｐは、スペクトログラムの勾配の分散を表し、これらは、ＳＴＦＴのフレーム長やフレームシフトに依存するであろう。スペクトログラムの勾配の実際の分布はガウス分布とは異なるが、ガウス分布を仮定することで問題の定式化及び解法を容易としている。後述するように、φ（Ａ）を用いてスペクトログラムのダイナミックレンジを圧縮することで、実際の状態と仮定とのギャップをある程度埋めることができる。

したがって、ＭＡＰ推定の目的関数は、以下のように書ける。

ここで、ベクトルｍはｍ_ｈ，ｉの集合であり、定数項は簡略化のため省略してある。

［Ｄ－２］補助関数を用いた更新ルールの導出
式（５）はｍ_ｈ，ｉの定積分形式であり、最適なｍは、ｍを連続値の変数であるとすると、∂Ｊ／∂ｍ_ｈ，ｉ＝０で求められる。ここで、∂Ｊ／∂ｍ_ｈ，ｉ＝０をより簡単に解くために、補助関数手法を用いる。補助関数は例えば、ＮＭＦ（Non-negative matrix factorization）やＨＴＣ（Harmonic-Temporal Clustering）において用いられており、当業者において公知の手法である。

補助関数を設計するにあたり、すべてのＡ，Ｂ，Ｘについて、

が成り立つことが、

が非負であり、X =(A + B)=2の場合に０となることから示される。

式（５）について不等式を適用することで、以下の補助関数

を導入し、この関数は、全てのｍ、補助パラメータＵ，Ｖについて、

を満たす。
したがって、以下の更新は、

Ｊを単調増加させる（ｋ：更新回数）。

［Ｄ－３］更新規則
∂Ｑ（ｍ，Ｕ^{（ｋ＋１）}，Ｖ^{（ｋ＋１）}）／∂ｍ_ｈ，ｉ＝０を

に変形することで、ｍ_ｈ，ｉのみの単純な更新式を生成できる。
式（７）から、Ｑ（ｍ^（ｋ），Ｕ，Ｖ）を最大化するＵ_ｈ，ｉ、Ｖ_ｈ，ｉは、

によって与えられる。式（１４）（１５）を式（１３）に代入し、ｍ_ｈ，ｉに代えてＨ_ｈ，ｉ、Ｐ_ｈ，ｉを更新変数とすることで、分離アルゴリズムは次のようになる。

（１）入力信号ｆ（ｔ）のSTFTであるＦ_ｈ，ｉを計算する。
（２）以下の式を用いて、レベル圧縮したrange-compressed パワースペクトログラムを計算する。

（３）全てのｈ、i（ｋ＝０）について、以下の式のように初期値を与える。

（４）更新変数△(k)を以下のように計算し、

ここで、

である。
そして、以下の場合に従って、Ｈ_ｈ，ｉ、Ｐ_ｈ，ｉを更新する。

（５）ｋをIncrementする。ｋ＜ｋ_ｍａｘ（ｋ_ｍａｘ：最大更新回数)の場合、ステップ４に戻り、それ以外の場合、ステップ６に進む。
（６）時間周波数マスクをバイナリ化する。これは以下と等価である。

（７）Ｈ^{（ｋｍａｘ）} _ｈ，ｉ，Ｐ^{（ｋｍａｘ）} _ｈ，ｉを、以下の式にしたがって、時間領域の波形に変換する。

ここで、ISTFTは逆STFTである。

［Ｄ－４］評価実験
実施例３の手法を用いて幾つかの実験を行った。
入力信号として、RWC 研究用音楽データベースからRWC-MDBより抜粋して使用した(16kHz サンプリング)。実験パラメータを表３に示す。バランスパラメータα、圧縮パラメータγは経験的に決定した。

RWC 研究用音楽データベースからRWC-MDB-P-2001 No.7より抜粋した６．２５［ｓ］区間に対して得られた調波成分Ｈ^{（ｋｍａｘ）} _ｈ，ｉ，Ｐ^{（ｋｍａｘ）} _ｈ，ｉのスペクトログラムを、図２に示す。γ＝０．３とした。更新を繰り返すにしたがって、スペクトログラムのエネルギーが水平方向の畝と垂直方向の畝とに分離されていくのがわかる。6.25[s]長信号（繰り返し数：５０）の計算時間は、laptop-PC with 1.20GHz Pentiumで2.3[s]であった。この時間は、実時間処理の３倍の速度である。

提案アルゴリズムの定量評価実験を行なった。RWC 研究用音楽データベースよりRWC-MDB-P-2001 No.18 and RWC-MDBJ-2001 No.16の前奏部を入力とし、MIDI形式データをパート別に分離し、各パートをWAV形式に変換してその信号の和を入力とした(16kHz サンプリング)。そして得られた分離結果の信号と各パート信号との相関を計算することで、Ｐ（ｘ，ｔ）とＨ（ｘ，ｔ）に含まれるエネルギー比率を算出した。エネルギー比率は、以下の式で算出した。

ここで、

であり、<>は相関計算を示す。

結果を表４，５に示す。

楽器や楽譜の情報なしで、音楽音響信号を調波成分と打楽器成分に分離する技術は、自動採譜や音楽検索など、音楽信号解析の様々なタスクを容易にする基礎技術として有用であり、またメロディー／リズムパートの強調や抑圧といったイコライジングなどの音楽信号の加工をも可能にする。これに関連して、リアルタイム調波音・打楽器音分離システムのGUI画面を図１１に示す。本システムでは、リアルタイムに分離した調波音・打楽器音のパワースペクトルを表示し、両者の音量バランスを調整しながら再生するという加工機能を実現した。本発明の対象となる音信号は音楽音響信号に限定されるものではなく、例えば、機械や装置から発生する工業的な音において異常音を取得することにも応用できる。

ポピュラー音楽のスペクトログラムを例示する図である。スペクトル成分の大きさは濃淡で表示されている。スペクトログラムを３次元表示した図である。時間方向（図において右側の軸）、周波数方向（図において左側の軸）に沿って、それぞれ滑らかなスペクトル成分があることが観察できる。尚、図１Ａは図１とは独立した図であり、別個のスペクトル成分に係るものである。時間周波数スペクトログラムの観測モデルを示す図である。左図は、調波音のスペクトログラムであり、時間方向に滑らか・周波数方向に急峻なスペクトル成分からなる。右図は、打楽器音のスペクトログラムであり、時間方向に急峻・周波数方向に滑らかなスペクトル成分からなる。左図のスペクトル成分と右図のスペクトル成分は、時間周波数平面上でスパースに存在している。入力スペクトルグラムと時間周波数マスクの乗算による、当該入力スペクトログラムの分離を示す図である。第１実施形態を示すブロック図である。第１実施形態で用いられるフィルタを示し、左図は、Ｈ（ｘ，ｙ）の特徴抽出フィルタ、右図は、Ｐ（ｘ，ｙ）の特徴抽出フィルタを示す。図６に示すフィルタの断面形状を示す。横軸は、Ｈ（ｘ，ｙ）の特徴抽出フィルタでは、時間の2次元フーリエ変換成分、Ｐ（ｘ，ｙ）の特徴抽出フィルタでは、P用フィルタでは周波数の2次元フーリエ変換成分である。縦軸は実際にフィルタの大きさを示し、大きいほど、その成分を通過しやすい。真ん中が0なので、三角窓、ガウス窓共にローパスフィルタになっている。第２実施形態を示すブロック図である。リアルタイム調波音・打楽器音分離を説明する図である。リアルタイム調波音・打楽器音分離における分離過程のある段階を示す図である。表示されたスペクトログラムにおいて、古い時間フレーム側において、周波数方向に滑らかなスペクトル成分が分離されていることが観察できる。リアルタイム調波音・打楽器音分離における分離過程のある段階を示す図である。表示されたスペクトログラムにおいて、古い時間フレーム側において、時間方向に滑らかなスペクトル成分が分離されていることが観察できる。第３実施形態の手法に基づいて繰り返し更新された調波成分（左）、打音成分（右）のスペクトログラムを示し、上から、ｋ＝０、ｋ＝３、Ｋ＝１０、Ｋ＝５０、２値化後の画像である。アルタイム調波音・打楽器音分離システムのGUI画面を示す。Method1は第１実施形態の手法、Method2は第２実施形態の手法を示す。

Claims

　音信号を分離する方法において、
　前記音信号のスペクトログラムを、時間周波数領域でのスペクトル成分の滑らかさに時間周波数平面上で方向性を備えた複数のサブスペクトログラムの和であると仮定し、
　各サブスペクトログラムの滑らかさの時間周波数平面上での方向に基づいて、時間周波数領域における前記音信号のスペクトル成分を、少なくとも１つのサブスペクトログラムに分配する少なくとも１つの分配係数を取得し、前記分配係数を用いて前記音信号のスペクトル成分から少なくとも１つのサブスペクトログラムを分離する、
　音信号の分離方法。
　前記分配係数は、時間周波数マスクである、請求項１に記載の音信号の分離方法。
　前記分配係数の取得は、
　前記音信号の各スペクトル成分について、各サブスペクトログラムの滑らかさの方向に基づいて、各サブスペクトログラムのスペクトル成分らしさのスコアを取得し、
　各スコアを指標として、分配係数を取得するステップと、
　からなる、請求項１，２いずれかに記載の音信号の分離方法。
　前記スコアを取得するステップは、
　前記音信号のスペクトログラムを、時間周波数平面上で各スペクトル成分がエネルギーに対応する濃度値を備えた画像とみなし、前記音信号のスペクトログラムから各サブスペクトログラムに属するスペクトル成分の特徴を夫々抽出するフィルタを設定し、
　各スペクトル成分についての、各サブスペクトログラムに対応するフィルタ処理後のフィルタ出力をスコアとする、
　請求項３に記載の音信号の分離方法。
　前記フィルタは、各サブスペクトログラムのスペクトル成分の滑らかな方向の濃度値を平滑化するローパスフィルタである、
　請求項４に記載の音信号の分離方法。
　前記音信号のスペクトログラムを、２個のサブスペクトログラムの和であると仮定し、
　前記スコアを比較し、スコアが大きい方の分配係数を１、スコアが小さい方の分配係数を０とする、
　請求項３乃至５いずれかに記載の音信号の分離方法。
　前記分配係数の取得は、
　分配係数をパラメータとして各サブスペクトログラムに分配された各スペクトル成分の滑らかさ指標の関数を含む目的関数を設定し、
　前記目的関数を最適化するパラメータを推定する、
　請求項１，２いずれかに記載の音信号の分離方法。
　前記分配された各スペクトル成分の滑らかさ指標は、
　着目したスペクトル成分と、当該着目スペクトル成分の時間周波数平面上での近傍にある分配されたスペクトル成分とのエネルギーの差に基づいて決定される、
　請求項７に記載の音信号の分離方法。
　前記滑らかさの指標の関数は、

である、請求項７，８いずれかに記載の音信号の分離方法。
　ここで、Ｋ：サブスペクトログラムの数、
　　　　　ｉ：周波数方向のインデックス
　　　　　ｊ：時間方向のインデックス
　　　　　ｆ_ｋ（ｘ）：滑らかさをはかるコスト関数
　　　　　ａ_ｍ，ｎ：時間周波数領域でのある点近傍での重み係数
　　　　　ｍ：周波数方向の近傍を表すインデックス
　　　　　ｎ：時間方向の近傍を表すインデックス
　　　　　ｇ（ｘ）：滑らかさ指標における、スペクトログラムのレベル圧縮関数
　　　　　Ｑ^（Ｋ） _ｉ，ｊ：サブスペクトログラムのスペクトル成分
である。
　前記目的関数が、前記音信号のスペクトル成分とパラメータとしての分配係数によって分配された各スペクトル成分の和との距離指標の関数を含む、
　請求項７乃至９いずれかに記載の音信号の分離方法。
　前記音信号のスペクトログラムを、Ｋ個のサブスペクトログラムの和であると仮定し、前記目的関数は、

である、請求項７乃至１０いずれかに記載の音信号の分離手法。
　ここで、Ｋ：サブスペクトログラムの数、
　　　　　ｉ：周波数方向のインデックス
　　　　　ｊ：時間方向のインデックス
　　　　　Ｄ（Ａ，Ｂ）：関数Ａと関数Ｂの距離指標
　　　　　φ（ｘ）：距離指標における、スペクトログラムのレベル圧縮関数
　　　　　Ｗ_ｉ，ｊ：観測スペクトル成分
　　　　　ｆ_ｋ（ｘ）：滑らかさをはかるコスト関数
　　　　　ａ_ｍ，ｎ：時間周波数領域でのある点近傍での重み係数
　　　　　ｍ：周波数方向の近傍を表すインデックス
　　　　　ｎ：時間方向の近傍を表すインデックス
　　　　　ｇ（ｘ）：滑らかさ指標における、スペクトログラムのレベル圧縮関数
　　　　　Ｑ^（Ｋ） _ｉ，ｊ：サブスペクトログラムのスペクトル成分
である。
前記目的関数において、

である、請求項１１に記載の音信号の分離方法。
　前記目的関数において、

である、請求項１１，１２いずれかに記載の音信号の分離方法。
　前記目的関数において、

である、請求項１１，１２いずれかに記載の音信号の分離方法。
　前記パラメータを推定するステップは、
　パラメータの更新と、当該パラメータによって分配された各サブスペクトログラムに対応するスペクトル成分の更新とを交互に繰り返すものである、
　請求項７乃至１４いずれかに記載の分離方法。
　前記音信号のスペクトログラムを、２つのサブスペクトログラムの和であると仮定し、
　前記時間周波数領域で隣接するパラメータによって分配されたスペクトル成分間のエネルギーの差の関数は、

　である、請求項７乃至１３いずれかに記載の音信号の分離方法。
　前記音信号のスペクトログラムを、２つのサブスペクトログラムの和であると仮定し、
　前記目的関数は、

である、請求項７乃至１３，１６いずれかに記載の音信号の分離方法。
　前記音信号のスペクトログラムを、２つのサブスペクトログラムの和であると仮定し、
　前記目的関数は、

　である、請求項７乃至９，１４いずれかに記載の音信号の分離方法。
　前記音信号を初期分析区間で時間周波数領域に変換してスペクトル成分を取得するステップと、
　前記音信号を１フレーム分だけ時間周波数領域に変換してスペクトル成分を取得し、当該スペクトル成分を前記分析区間に加えるステップと、
　前記分析区間のスペクトル成分を用いて、パラメータの推定を行うステップと、
　前記分析区間で最も古い１フレームのスペクトル成分を、推定されたパラメータを用いて分離し、
　分離されたスペクトル成分を時間周波数領域に変換する、
　請求項７乃至１８いずれかに記載の音信号の分離方法。
　推定された分配係数を２値化するステップを含む、
　請求項７乃至１９いずれかに記載の音信号の分離方法。
　２値化の強度が可変である、
　請求項２０に記載の音信号の分離方法。
　前記複数のサブスペクトログラムの少なくとも１つは、周波数方向に滑らかなサブスペクトログラム、あるいは、時間方向に滑らかなサブスペクトログラムである、請求項１乃至２１いずれかに記載の音信号の分離方法。
　前記複数のサブスペクトログラムは、周波数方向に滑らかな第１サブスペクトログラムと、時間方向に滑らかな第２サブスペクトログラムと、を含む、請求項２２に記載の音信号の分離方法。
　前記周波数方向に滑らかなサブスペクトログラムは、非調波的な成分であり、時間方向に滑らかなサブスペクトログラムは、調波的な成分である、請求項２２、２３いずれかに記載の分離方法。
　前記音信号は音楽信号であり、前記非調波的な成分は、打楽器音である、請求項２４に記載の音信号の分離方法。
　前記分離方法は、分離された少なくとも１つのサブスペクトログラムのスペクトル成分を強調あるいは抑圧するステップを備えている、
　請求項１乃至２５いずれかに記載の音信号の分離方法。