3-5-7

スペシャル・セッション〔音楽音響と音楽関連産業〕

調波時間因子分解に基づく音楽事前情報付き多重音解析 ∗☆四方 紘太郎 (東大・工)，高宗 典玄，中村 友彦 (東大院・情報理工)，亀岡 弘和 (東大院・情報理工/NTT CS 研)

1 はじめに複数の音源の信号が混合された観測信号から個々の音源に関する情報 (基本周波数，発音開始時刻，パワーなど) を抽出する処理である多重音解析は，音楽
情報処理における重要課題の一つであり，自動採譜や音楽音響信号加工などの基礎技術となりうる．
マイクロホンアレイ入力からのブラインド音源分
離では音源の空間的な手がかりを有効利用することができるが，モノラル音響信号を対象とした音源分離や多重音解析では空間的な手がかりに代わる何らかの仮定が必要である．聴覚情景分析における知見をヒントにしたアプローチである調波時間構造化クラスタ
リング (Harmonic-Temporal Clustering; HTC) [1,2] は，人間が音をひとまとまりの音 (音脈) として知覚する
要件 (調波性，連続性，同時性，同期性など) を時間周波数成分の局所的な制約として記述し，当該要件を満たすように観測信号の時間周波数成分を時間周波数
平面上でクラスタリングしようというアイディアに基づいている．HTC ではこのアイディアを，各音脈に対応するスペクトログラムを拘束つき混合正規分布
モデルとして記述し，それらを重畳したもので観測ス
ペクトログラムにフィッティングするアプローチとして実現している．一方，モノラル音響信号を対象とした多重音解析手法として有効なアプローチとして近年注目されている非負値行列因子分解 (Non-negative
Matrix Factorization; NMF) 法 [3] では，限られた種類の音高の楽音がそれぞれ異なるタイミングで繰り返し生起するという音楽特有の性質に着目し，限られた種類のスペクトルテンプレートの適当な重み付き和ですべての時刻の観測スペクトルを表せるはず，という仮定がベースとなっている．従って，観測スペク
トログラムを非負値行列と見なし，これを二つの非負値行列の積（各スペクトルテンプレートを表す基
底行列と，各荷重係数を要素にもつアクティベーショ
ン行列）に分解することにより観測スペクトログラム
から各スペクトルテンプレートと各時刻におけるそれらの荷重係数を同時推定することができ，観測ス
ペクトルを音高ごとのスペクトルに分解することが可能となるわけである．上記のとおり上述の二つのアプローチでは着目している手がかりが異なる．相対的には前者のアプローチでは局所的，後者のアプ
ローチでは大域的な楽音の性質に基づいていると言え，いずれの性質も多重音解析を解決する上で本質的かつ有用な手がかりとなる．本稿では，上述の二つの性質を同時に取り入れた新しいスペクトログラム
モデル，および当該モデルに基づく多重音解析手法「調波時間因子分解法」を提案する．
ところで，近年の音楽情報検索に関する国際会議
や音楽情報検索に関わる各種タスクの国際コンテス
ト [4] において，調推定，和音推定，拍推定，オンセッ
ト推定などの手法の研究が急速に進展している．調，和音，拍，オンセットなどの情報が高い精度で得られるのであれば，多重音解析において極めて有用な補助情報となりうる．そこで本稿ではさらに，調推定，
和音推定，拍推定，オンセット推定の各手法により得られる各種情報を前記提案モデルのパラメータ推論において補助情報として活用する枠組を提案する．
音源分離において，ユーザのアシストにより分離精∗度を向上できるようにすることを目的としたユーザ
ガイドつき音源分離 [5, 6] や，楽譜を補助情報とする
楽譜ガイドつき音源分離 [7, 8] などの研究が進められているが，本稿で提案する枠組も補助情報ガイドつき音源分離の一種として位置づけられる．
以下，正規分布と Dirichlet 分布，Poisson 分布の確率密度関数を N, Dir, Pois と表記する．

2 音楽スペクトログラムの確率モデル化
2.1 楽音信号のウェーブレット変換
本節ではまず，[1] に倣い，楽音信号のウェーブレッ
ト変換を導出する．
楽音の多くは，擬似周期信号 (局所的には周期的と見なせ，周期や調波成分のパワーが滑らかに時間変化する信号) と見なせる．そこで，音源 k の信号モデ
ルとして n 次調波成分の瞬時位相が nθk (u) + φk,n ，瞬時振幅が ak,n (u) の擬似周期信号の解析信号表現N∑
fk (u) =ak,n (u)e j(nθk (u)+φk,n )(1)n=1を考える．ただし，u は時刻，φk,n は初期位相である．
紙面の都合上導出を省略するが， fk (u) のウェーブレッ
ト変換 Wk (x, t)∑は，ak,n (t)Ψ∗ (ne−x µk (t))e j(nθk (t)+φk,n )Wk (x, t) =(2)nと表される．ただし， x は対数周波数，t は時刻，Ψ
は中心周波数が 1 のアナライジングウェーブレットのFourier 変換，µk (t) = θ̇k (t) は瞬時基本周波数をそれぞれ表している．今，Ψ が対数正規分布型の関数となる
ようなアナライジングウェーブレットを選び，また，Ωk (t) = ln µk (t) と置くと，Wk (x, t) は∑(x−Ωk (t)−ln n)2(3)
Wk (x, t) =ak,n (t)e− 4σ2 e j(nθk (t)+φk,n )nと表される．ただし，σ は対数正規分布のスケールパ
ラメータに対応する定数である．ここで，n, n′ (n , n′ )の指数項の重なりがほとんどない (調波成分が互いに重ならない) と仮定できるならば，|Wk (x, t)|2 は近似的に，∑(x−Ωk (t)−ln n)2|Wk (x, t)|2 ≃|ak,n (t)|2 e− 2σ2(4)nと表せる．ここまでは HTC で採用された調波時間構造モデルと同一であり，混合正規分布モデル (Gaussian
mixture model; GMM) と同形な関数となっていることが分かる．

2.2 音源スペクトログラムモデル
ここで，NMF におけるモデル化の考え方を上記モデルに取り入れる．NMF では，各楽音のスペクトルの形状は時不変で，スケールのみが時間変化すると仮定されるが，上記調波時間構造モデルにおいて，|ak,n (t)|2を時刻 t に依らない変数と調波成分インデックス n に依らない変数の積の形√(5)|ak,n (t)|2 = wk,n Uk (t)/ 2πσHarmonic Temporal Factor Decomposition Incorporating Musical Prior Information for Multi-pitch analysis by Kotaro
Shikata, Norihiro Takamune, Tomohiko Nakamura (The University of Tokyo), Hirokazu Kameoka (The University ofTokyo/NTT CS Lab.)

日本音響学会講演論文集- 1049 -2014年3月

Fig. 1 バイオリンのビブラート音のパワースペクト
ログラムに分解できるとすると，|Wk (x, t)|2 は|Wk (x, t)|2 = Hk (x, t)Uk (t)(6)∑ wk,n(x−Ωk (t)−ln n)2
Hk (x, t) :=e− 2σ2(7)√2πσnと表される．上式の Hk (x, t) は GMM 型の関数で表される音源 k の時刻 t におけるスペクトル形状を表しており，これにアクティベーション Uk (t) が乗じられた形になっている．以下では w と U のスケールの任意性を除くため，∑wk,n = 1(8)nを仮定しておく．このとき，wk,n は調波成分のパワー比を表すパラメータとなる．既存モデルとの関連は 2.4
節 で詳しく述べるが，以上のモデルは HTC と NMFで用いられている音源スペクトログラムモデルの特長を併せ持っている．

2.3 観測スペクトログラムの確率モデル
これまでは，スペクトログラムモデルを連続時間および連続周波数の関数として導いたが，実際に計算機で算出されるスペクトログラムは離散時刻および離散周波数ごとの値として得られる．そこで，以下では，等間隔に離散化された対数周波数を xl (l = 1, . . . , L)，
等間隔に離散化された時刻を ti (i = 1, . . . , I) とし，観測スペクトログラムを Y(xl , ti ) と表す．
HTC や NMF と同様にパワースペクトルの加法性を仮定すれば，音楽スペクトログラムは複数の音源
スペクトログラム（式 ∑(6)）を重畳したもの
X(xl , ti ) =Hk (xl , ti )Uk (ti )(9)kとして表せる．このスペクトログラムモデル X(xl , ti )には，実際の音源信号における式 (1) で定義した擬似
周期性の仮定からの逸脱による誤差，調波間干渉を無視したことによる誤差，パワーの加法性の仮定に起因する誤差，背景雑音の存在に起因する誤差など，さまざまな要因の誤差が混在する．提案法の枠組では，これらの複合的な誤差要因を一つ一つ詳細にモデル化することはせず，まとめて一挙に確率的な現象と見なし，Y(xl , ti ) はY(xl , ti ) ∼ Pois(Y(xl , ti ); X(xl , ti )) (∀l, ∀i)(10)により生成されたものと仮定する．なお，この仮定の下で，X(xl , ti ) を変数と見なして最尤推定する問題
は，スペクトル間の乖離度を I ダイバージェンスと呼ぶ歪み尺度を規準とした Y と X の最適フィッティン
グ問題と等価となる．

日本音響学会講演論文集

2.4 従来モデルとの関連
本節では，先に提案したスペクトログラムモデル
X(xl , ti ) と従来モデルとの関連について述べる．
まず，Hk (xl , ti ) に関し，式 (7) のようなパラメトリッ
クな関数を仮定せずインデックス k, l, i ごとの Hk (xl , ti )の値をパラメータと見なせば式 (9) は「可変基底 NMF」
[9] において仮定されるスペクトログラムモデルと同一となる．また，Hk (xl , ti ) に対し時不変となるよ
うな拘束を置けば通常の NMF [3] において仮定されるスペクトログラムモデルと同一となる．さらに各
Hk に調波構造をなすスペクトルを仮定すれば，「調波
NMF」[10, 11] において仮定されるスペクトログラム
モデルと同一となる．次に，式 (7) において，Ωk (t) に対し時不変となるような拘束を置けば，[12, 13] のス
ペクトログラムモデルと同一となる．最後に，Uk (t)を拘束つき GMM 型の関数で記述すれば，HTC [1, 2]
において仮定されるスペクトログラムモデルと同一となる．以上より提案モデルは NMF と HTC の双方と親戚関係にあることが分かる．

3 音楽事前情報の組み込み
3.1 事前分布としての音楽事前情報
音響信号から調や和音，拍・オンセット時刻などの情報を推定するための手法の研究は近年急速に進展している [4]．調や和音，拍・オンセット時刻などの情報は多重音解析においては有用な補助情報になりうるため，既存手法を用いてこれらの推定を前段で行い，その結果を補助情報として活用すれば高い精度で多重音解析を行える可能性がある．上記の前段
処理では推定誤りを含むこともありえるが，補助情報としての信頼度を確率と捉えれば，各パラメータの事前確率として推論に組み込むことが可能である．
以上の補助情報を事前確率として表せれば，複数の推定結果の同時活用も可能である．例えば，調と和音の情報が与えられたときには，その調で出現しやす
く，かつその和音で出現しやい音高が選ばれやすいはずである．この「かつ」に相当する演算は，両方の条件を表す事前分布の積で表される．このように，複数の条件を同時に成立させるには，条件を表す事前分布の積として表現すれば良く，Products of Experts [14]の考え方が採用できる．
3.2 音楽事前情報の事前分布による設計
3.1 節の方針に従って，いくつかの音楽事前情報について事前分布を設計する．
3.2.1 基本周波数の事前分布
弦楽器や管楽器の楽器音は，1 つの音符を演奏していても基本周波数が変化しうる．例えば，バイオリンなどの奏法であるビブラートは，基本周波数がある音高に対応する基本周波数の周りで小刻みに振動しつつ，連続的に変化する (Fig. 1)．3.1 節で議論した通り，k 番目の音源スペクトログラムモデルの対数基
本周波数 Ωk (ti ) をビブラートで演奏されている音符の
基本周波数とすると，この音高の基本周波数の周りでの変動を表す確率分布 qg と基本周波数の連続的な変化を表す確率分布 q s の積として事前分布を設計できる．
このような確率分布として，
qg (Ωk (ti )) = N(Ωk (ti ); mk , νk2 )(11)
qc (Ωk (ti )|Ωk (ti−1 )) = N(Ωk (ti ); Ωk (ti−1 ), τ2k ) (12)を用いることができる．ここで，mk , νk は k 番目の音
源スペクトログラムモデルの音高に対応する対数基
本周波数と対数周波数軸上での分散を表し，τ2k は対数基本周波数の時間変化量の分散を表す．定性的に説明すると，qc は対数基本周波数の軌跡が時間に関し- 1050 -2014年3月

て滑らかであることを意味し，qg は時間に関係なく
対数基本周波数が与えられた音高 mk 周辺に存在することを意味する．qg , qc を用いて Ωk (t) の事前分布は，
p(Ωk (ti )|Ωk (ti−1 )) ∝ qg (Ωk (ti ))αg qc (Ωk (ti )|Ωk (ti−1 ))αc(13)と書ける．ここで，αg , αc は qg , qc の事前分布への寄与を調節するパラメータである．

3.2.2 調と和音の事前分布
調性のある楽曲においては，曲の部分ごとに調や和音が存在し，調と和音に従って出現する音高に偏りがある．これは，各音高に対応する音源スペクトログ
ラムモデルのアクティベーション Uk (t) の事前分布としてモデル化できる．また，従来よく用いられてきた事前分布も統合的に扱える．先行研究では，時間に関するスパース性を仮定することがあり，経験的に有効であることが知られている．
これら 2 つの条件を満たすような事前分布は，楽曲全体の音量
C と時間方向に正規化された音量∑
Bk (t∑i ) ( i Bk (ti ) = 1)，音高方向に正規化された音量
Ak ( k Ak = 1) を用いて，
Uk (ti ) = CAk Bk (ti ),(14)
A := [A1 , · · · , AK ]⊤ ∼ Dir( A; β)(15)
Bk := [Bk (t1 ), · · · , Bk (tI )]⊤ ∼ Dir(Bk ; γk )(16)とモデル化できる．このように，調や和音による各音高の生起しやすさは， A の事前分布のハイパーパラ
メーター β を適切に設定することにより反映することが出来る．また，時間に関するスパース性も Bk の事前分布としてハイパーパラメータ γk を適切に設定し，反映することができる．

は，補助変数を用いて J(Θ) の下界（補助関数）を作り，その補助関数を補助変数と Θ について交互に最大化することにより J(Θ) を単調増加させる．
(19) の問題となる部分は，対数関数が凹関数であるため Jensen の不等式を 2 回用いて，ln p(Θ|Y)∑∑∑wk,n ϕk,n (xl , ti )≥ Y(xl , ti )
λi,l,k Uk (ti ) ln nλi,l,ki,lk∑∑wk,n ϕk,n (xl , ti )≥ Y(xl , ti )
λi,l,k ζi,l,k,n Uk (ti ) ln
λi,l,k ζi,l,k,n
i,l
k,n

(21)
と 補 助 変 数 λi,l,k , ζi,l,k,n ∈ [0, 1] を 使 用∑し て 下 界を求める∑こ と が で き る ．こ こ で ，∀i, l, k λi,l,k =1, ∀i, l, k n ζi,l,k,n = 1 である．等号成立条件は，∑
Uk (ti )wk,n ϕk,n (xl , ti )
λi,l,k = ∑ n,(22)n,k U k (ti )wk,n ϕk,n (xl , ti )wk,n ϕk,n (xl , ti )
ζi,l,k,n = ∑(23)
n wk,n ϕk,n (xl , ti )である．したがって， J(Θ) の補助関数は，
J + (Θ, {λi,l,k , ζi,l,k,n }k,n )∑∑wk,n ϕk,n (xl , ti )=c Y(xl , ti )
λi,l,k ζi,l,k,n Uk (ti ) ln
λi,l,k ζi,l,k,n
i,l
k,n
∑(
)
+
αg ln qg (Ωk (ti )) + αc ln q s (Ωk (ti )|Ωk (ti−1 ))
k,i

+

∑

ln Dir( A; β) +i

4 事後確率最大化によるパラメータ推定ア
ルゴリズム
4.1 目的関数
2，3 章の議論に基づき，音楽パワースペクトログ
ラム Y := {Y(xl , ti )}l,i が与えらたときのパラメータ推定法を考える．事前分布を含む最尤推定は MAP 推定とよばれ，パラメータ Θ = {wk,n , Ωk (ti ), Ak , Bk (ti ), C}iを用いて，
argmax ln p(Θ|Y) = argmax (ln p(Y|Θ) + ln p(Θ)) (17)ΘΘとなる事後確率 ln p(Θ|Y) を最大化する Θ を推定する．
ここで，ln p(Y|Θ) は尤度， p(Θ) は事前分布を表し，ln p(Y|Θ)


∑

∑
(18)
=c
Y(xl , ti ) ln 
Hk (xl , ti )Uk (t)
k

i,l

∑

=

i,l




∑

Y(xl , ti ) ln  wk,n ϕk,n (xl , ti )Uk (t)

(19)

k,n

ln p(Θ)∑()αg ln qg (Ωk (ti )) + αc ln q s (Ωk (ti )|Ωk (ti−1 ))
=c
k,i

+

∑

ln Dir( A; β) +i∑

ln Dir(Bk ; γk )

(20)

k

と書ける．ここで，=c は定数を除いて一致することを表す．(17) の右辺を arg maxΘ J(Θ) とすると， J(Θ)
を最大化するような Θ の更新則を求めれば良い．しかし，(19) には対数関数の中に k, n に関する和が存在するため，閉形式の更新則を求めることは容易ではない．
この場合には，補助関数法 [15] を用いることにより閉形式の更新則を導くことができる．補助関数法

日本音響学会講演論文集∑

ln Dir(Bk ; γk )

(24)

k

と書ける．

4.2 パラメータの更新則
補助関数法では，パラメータの更新と補助変数の更新を繰り返して最適化を行う．補助変数の更新は等号
成立条件 (22),(23) にしたがって更新すれば良い．パ
ラメータの更新則は，補助関数 J + (Θ, {λi,l,k , ζi,l,k,n }k,n )を各パラメータで偏微分した値が 0 となるように求めればよい．C は音楽パワースペクトログラム Y から一意に決まり，∑C = (x2 − x1 )Y(xl , ti )
(25)
i,l

である．パラメータの更新則はそれぞれ∑
l,i Y(xl , ti )λi,l,k ζi,l,k,n
wk,n ← ∑(26)
l,i,n Y(xl , ti )λi,l,k ζi,l,k,n
(α
)−1
αg
c
Ωk (ti ) ← 2 D⊤ D + 2 E I + (diag pl,k,n )
ν

τ
∑

 mk αc 2 1I +(27)
(xl − ln n) pl,k,n νl,n∑l,i,n Y(xl , ti )λi,l,k ζi,l,k,n + βk − 1
Ak ← ∑ ∑(28)
Y(xl , ti )λi,l,k ζi,l,k,n + βk − 1)(k
∑ l,i,nl,n Y(xl , ti )λi,l,k ζi,l,k,n + γi,k − 1
Bk (ti ) ← ∑ ∑(29)(l,n Y(xl , ti )λi,l,k ζi,l,k,n + γi,k − 1)tと導出できる．ここで，diag はベクトルを行列の対角要素に順に並べる演算を表し，1I は I 次元で要素が全て 1 のベクトル，E I は I 次元の単位行列である．
pl,k,n はY(xl , ti )λi,l,k ζi,l,k,n
pl,k,n,i =(30)
σ2- 1051 -2014年3月

Fig. 2 KL-ダイバージェンス基準 NMF による基底
スペクトログラムの推定結果 (基底数 3)

Fig. 3 D♭,F,A♭ それぞれの音源スペクトログラムモ
デルと基本周波数の時間軌跡推定結果 (調に関する音
楽事前情報なし)．赤線は各時刻の Ωk (ti ) の推定値を表す．
を用いて，pl,k,n := [pl,k,n,1 , · · · , pl,k,n,I ]⊤ と書け，I × I
行列 D は 0 0 · · · · · · 0−1 1


−1
1


(31)
D = 
. . . . . . 
0−1 1である．これを調波時間因子分解法と名付ける．

5 多重音解析動作実験
本研究において提案した手法を MATLAB で実装
し，多重音解析動作実験を行った．本章ではそれに対し定性的な考察を試みた．

5.1 ビブラートの基本周波数推定実験
提案モデルは基本周波数の時間変動を考慮しており，ビブラート音の基本周波数の推移を追うことが可能であると期待される．そこで，RWC 楽器音データ
ベース [16] を用いて，D♭，F，A♭ のバイオリン音源のビブラート音を合成した音源（サンプリング周波数
16 kHz）を作成し，基本周波数の推定に関し，KL-ダ
イバージェンス基準の NMF(基底数 3, パラメータ更
新回数 100 回) との比較実験を行った．
スペクトログラムを音響信号から得る際に，高速近似連続ウェーブレット変換 [17] を用いた．この変換を用いると，スペクトログラムの時間シフトはある程度自由に設定でき，本実験では ti − ti−1 = 16 msとした．アナライジングウェーブレットは，対数正規分布型のウェーブレット [1] を用いて x1 = ln(55), xl −xl−1 = ln(2)/120 となるようにスケール値を設定した．
提案法のパラメータは (N, K, τ, ν, σ, αg , αc ) = (8, 73, 2×
104 , 1.5, ×104 , 0.02, 1, 1)，γk = (1 − 3.96 × 10−6 )1I ，β =(1 − 2.4 × 10−3 )1K とし，パラメータ更新回数は 2 回とした．
観測 (Fig. 1) で見られたビブラートは，Fig. 2 のよ
うに NMF のモデルスペクトログラムで平坦な形になり，ビブラートの基本周波数の時間的な変化を捉えられていないことがわかる．それに比べ，提案手法ではFig. 3 のように基本周波数パラメータが時間領域でビブラート時の基本周波数の推移と同様の軌跡を描いていることが確認できる．

Fig. 4 音楽事前情報なし (左), あり (右) のときのア
クティベーション推定結果
5.2 音楽事前情報の効果確認の実験
次に，調の情報を事前分布として取り入れた場合，調に対応する音階外の音が抑制されるかを確認する実験を行った．ここで調は D♭ であると仮定し，音階外である D の音に対するアクティベーションが抑制されているかを確認した．音階内の音について β = 1 − 2.4 ×
10−3 とし，音階外の音については β = 1 − 3.0 × 10−3とした．その他のパラメータは 5.1 節と同じである．
Fig. 4 右が示す通り，調の音楽事前情報を反映することによって，調の音楽事前情報を用いていない
Fig. 4 左に比べて音階外である D の音高のアクティ
ベーションが抑制された．したがって，音楽事前情報に音高の誤推定を抑制する効果があることを確認で
きた．

6 結論
本研究では，NMF と HTC の性質を同時に取り入れた新たな音源スペクトログラムモデルを提案し，そのモデルに基づく多重音解析手法である調波時間因
子分解法を提案した．また，音楽事前情報を Productof Experts の形で事前分布としてモデルに導入できることを示し，事後確率を最大化するような推定アルゴ
リズムを導出した．基本周波数の時間変化に関して
は，ビブラート音に対する基本周波数推定実験から，基本周波数パラメータが時間領域で基本周波数の推移に追従できることを確認した．また，音楽事前情報による多重音解析性能への効果を確認するために，調を音楽事前情報として用いた実験を行い，音階外の音高の誤推定が抑制されることを確認した．
今後は，提案手法の定量的な評価を行い，音楽事前情報を取り入れた高精度な多重音解析ソフトウェアを開発することが課題である．

参考文献[1] H. Kameoka, “Statistical approach to multipitch analysis,” Ph.D.
thesis, University of Tokyo, 2007.
[2] H. Kameoka et al., In Proc. of IPSJ SIG Technical Report,
pp. 27–32, 2003.
[3] P. Smaragdis et al., In Proc. of IEEE, 2003.
[4] http://www.music-ir.org/mirex/wiki/MIREX_HOME.
[5] P. Smaragdis et al., In Proc. of WASPAA, pp. 69–72, 2009.
[6] A. Ozerov et al., In Proc. of ICASSP, pp. 257–260, 2011.
[7] R. Hennequin et al., In Proc. of ICASSP, pp. 45–48, 2011.
[8] U. Simsekli et al., In Proc. of EUSIPCO, pp. 2639–2643, 2012.
[9] M. Nakano et al., In Proc. of LVA/ICA, pp. 149–156, 2010.
[10] S.A. Raczyński et al., In Proc. of ISMIR, pp. 381–386, 2007.
[11] E. Vincent et al., In Proc. of ICASSP, pp. 109–112, 2008.
[12] K. Yoshii et al., In Proc. of ISMIR pp. 309–314, 2010.
[13] 阪上 他, 情処全大, 4T-8, 2013.
[14] G.E. Hinton, Neural Comput., 14(8), pp. 1771–1800, 2002.
[15] 亀岡 他, 情処研報, 49, pp. 139–144, 8, 2006.
[16] M. Goto, In ICA, pp. I-553–556, 2004.
[17] 亀岡 他，“信号処理方法及び装置,” 特開 2008-281898, 11, 2008.

日本音響学会講演論文集- 1052 -2014年3月