2-10-12高速近似連続ウェーブレット変換による振幅スペクトログラムに対する実時間位相推定法 ∗◎中村友彦（東大院・情報理工），亀岡弘和（東大院・情報理工，NTT・CS 研）

1 はじめに定 Q フィルタバンクとしても知られる連続ウェー
ブレット変換（CWT）は，対数周波数スケールで一様な解像度をもつ時間周波数表現を与える．人間の聴覚システムも特に高周波数帯域で対数的な周波数解像度を持つことが知られており，聴覚システムに着想を得た音響信号処理手法の開発には CWT で得られたスペクトログラム上での処理が望ましいと考えられる．実際に，複素スペクトログラムを用いたマル
チチャネル音源分離 [1]，振幅スペクトログラムを用いた多重基本周波数推定 [2–4] や歌声分離 [5] などでCWT の有用性が確認されている．音源分離や音響信号加工などの時間領域の音響信号を生成することが目的のアプリケーションにおいては，分離や加工された振幅スペクトログラムを音響信号に変換する必要がある．そのため，与えられた振幅スペクトログラムに対し適切な位相を与える必要がある．
最初に提案された振幅スペクトログラムからの位相推定アルゴリズム [6] は，所望の振幅スペクトログラムに現在の位相の推定値を割り当てた複素スペクトログラムに逆 CWT と CWT を行うステップと，複素スペクトログラムの振幅部分のみを所望の振幅スペクトログラムと置換するステップからなり，これらのステップを繰り返し行うことで位相を推定する．しかし，CWT は計算量が高いためこのアルゴリズムは長時間の計算を要し [7]，実際のアプリケーションに用いるためには計算量の削減が重要となる．また，アルゴリズムの収束性も実用上重要である．CWT と逆CWT の様々な高速計算法は近年相次いで提案されて
おり [8–10]，CWT の代わりにこれらの高速計算法を利用することで計算法を削減できるはずである．しかし，そのようなアルゴリズムでは収束性が保証されるかどうかが不明であった．
そこで，我々は補助関数法 [11] と呼ばれる最適化原理に基づき導出したアルゴリズムが [6] で提案されたアルゴリズムと一致することを示し，冗長な線形変換であればアルゴリズムの収束性が保証されることを明らかにした [7]．この結果を元に，CWT の代わりに高速近似 CWT [8] を用いることで，収束性を保証しつつ計算量の削減されたアルゴリズム（高速位相推定アルゴリズム）を提案した．
高速近似 CWT を含め多くの CWT 高速化アルゴリ
ズムでは信号全体の高速 Fourier 変換（FFT）を用いることが前提となるため，信号長が長いほど高速位相推定アルゴリズムの空間計算量は増加する．その∗ため、実際的な信号長のスペクトログラムを計算するためには多くのメモリが必要となり，オーディオプレーヤーなどのメモリ制限のある装置では処理可能な音響信号が限られる．また，原理上実時間動作を要求するアプリケーションにそのまま適用することはできない．
そこで本稿では，高速位相推定アルゴリズムを拡張し，信号長によらず一定の空間計算量をもち音響信号を逐次処理可能な実時間位相推定アルゴリズムを提案する．短時間 Fourier 変換で得られた振幅スペクトログラムからの実時間位相推定アルゴリズム [12] と同様に，提案アルゴリズムでは音響信号をオーバーラップのある固定長のフレームに分割し，各フレームに高速近似 CWT を適用して得られたスペクトログラ
ムを対象とする．高速位相推定アルゴリズムは位相全体へ定数を加えたものも解として許容しうるため，隣り合うフレーム同士で独立に高速位相推定アルゴ
リズムを適用した場合には少なくとも異なる初期位相をとりうる．そのため，フレーム同士のオーバーラップ部分で不連続な信号成分となる位相が推定される可能性がある．そこで，提案アルゴリズムでは隣り合うフレーム間のオーバーラップ部分の信号成分が同一である（無矛盾）ことを考慮することにより，不連続な信号となることを抑制する．

2 高速近似ウェーブレット変換を用いた位相推定アルゴリズム
2.1 高速近似 CWT
まず，高速近似 CWT について簡単に紹介する．
CWT は，スケーリングされたアナライジングウェー
ブレットが各フィルタのインパルス応答に対応するフィルタバンクと解釈できる．Morlet や対数正規分
布型ウェーブレット [3] などのアナライジングウェー
ブレットを用いた場合，各サブバンドフィルタの周波数特性の主要な部分が局所的に分布する．そのため，主要な値が存在する周波数領域のみを計算に用いることで時間信号の CWT を高速に計算できる．高速近似 CWT はこの考えに基づき提案された．
T 次元の時間信号を f ，T 次元の DFT 行列を FT とする．高速近似 CWT では，時間信号全体の FFT，FT fを求めた後に各サブバンドフィルタの主要な部分が存在する領域 k ∈ [B, B + D − 1] に帯域制限を行う．ここで，k = 0, · · · , T −1 は角周波数インデックスである．ある対数周波数インデックス ω(= 0, · · · , Ω − 1) に対応するサブバンドフィルタについて考える．帯域制限を表

Real-Time Phase Estimation from Magnitude Spectrogram of Fast Approximate Continuous Wavelet Transform by
NAKAMURA, Tomohiko (The University of Tokyo) and KAMEOKA, Hirokazu (The University of Tokyo / Nippon
Telegraph and Telephone Corporation)

日本音響学会講演論文集- 933 -

2016年3月

す行列は，全要素が 0 の D×B 次元の行列 0D×B と D 次
元の単位行列 ID を用いて，L := [0D×B , ID , 0D×(T −D−B) ]と書ける．帯域制限された時間信号の FFT に対して，帯域制限された周波数帯域の当該フィルタの周波数特性 ψω ∈ CD を乗算し，diag(ψω )LFT f を得る．ここで，diag(ψω ) は ψω を対角成分に並べた対角行列を表す．
帯域制限により正規化角周波数が [2πB/D, 2πB/D+2π]に分布するため，diag(ψω )LFT f を巡回的にシフトさせ，先頭の成分の位相が 0 となるようにする．この操作は以下の行列 C で書ける．

0(B−(p−1)D)×(pD−B)
IB−(p−1)DC :=  (1)I(pD−B)0(pD−B)×(B−(p−1)D)ここで， p は pD ≤ B + D < (p + 1)D となる最大の整数である．巡回シフトを行ったものに対し D 次元の逆 FFT を行うことで，当該サブバンドでのスペクトログラムが得られる．
これらをまとめると，対数周波数インデックス ωHに対応する高速近似 CWT Wω = F D
Cdiag(ψω )L を用⊤
いて，高速近似 CWT は W := [W0 , · · · , WΩ−1 ]⊤ と表せる．ここで，H は当該変数の Hermitian 共役を表す．
この逆変換（逆高速近似 CWT）は W の擬似逆行列
W + で表される．D, B はサブバンド毎に異なってもよいが，簡単のため以下では全てのサブバンドで D, Bは同一とする．

2.2 高速位相推定アルゴリズム
W は T 次元複素ベクトルから DΩ 次元複素ベクトルへの変換であるため，T < DΩ の場合，高速近似CWT で得られるスペクトログラムは時間信号の冗長な表現である．したがって，スペクトログラムの集合 W は CDΩ の部分空間である．ここで，逆高速近似 CWT の後高速近似 CWT を適用する操作 WW + を考える．WW + は W への直行射影であるため，DΩ次元の複素ベクトル y に対し y − WW + y は y から Wへの距離を表す．したがって，この距離が小さければ小さいほど y はより「スペクトログラムらしい」とみなせ， y = WW + y のとき y はスペクトログラムであり，対応する時間信号が存在する．そのため，与えられた振幅スペクトログラムに対して，この距離を最小化する位相を求める問題として位相推定問題は定式化できる．
振幅スペクトログラム a ∈ RΩT≥0 が与えられたとき，位相推定問題は位相の推定値 ϕ ∈ [0, 2π)ΩT を用いて以下のように定式化できる．
min sH (a, ϕ)(IΩT − WW + )s(a, ϕ)ϕ∈[0,2π)ΩT(2)

ここで， s(a, ϕ) は振幅 a，位相 ϕ の複素ベクトルを表す．この問題は s(a, ϕ) 自体をパラメータ s =[s0,0 , s0,1 , · · · , s0,D , s1,0 , · · · , sΩ−1,D ]⊤ とみなせば，各時間周波数成分 sω,t に関し |sω,t |2 = a2ω,t という 2 次制約を持つ 2 次計画問題に書き換えられる．しかし，WW +の行数および列数はスペクトログラムの要素数と同日本音響学会講演論文集

Fig. 1: オンライン高速近似 CWT の処理フロー．
一であるため計算コストが非常に高く，実際的にはこの 2 次計画問題を直接解くことは難しい．
そこで，我々は [7] で補助関数法と呼ばれる最適化原理に基づき高速位相推定アルゴリズムを提案した．
導出の詳細は省くが，下記の 2 つの更新式を交互に更新することにより，広義単調減少させることができる．

s̃ ←WW + s(a, ϕ)(3)

ϕ ←∠s̃

(4)

ここで，∠s̃ は複素ベクトル s̃ の各要素の偏角を[0, 2π)ΩT のベクトルとして返す演算子である．(3) 式
は，現在のスペクトログラム推定値に対し逆高速近似 CWT を行ったのちに高速近似 CWT を適用することを表す．(4) 式は，(3) 式の操作で得られたスペクトログラム s̃ の位相を新たな位相の推定値とすることに対応する．

3 位相推定アルゴリズムの実時間化
3.1 オンライン高速近似 CWT
前節までの処理では高速近似 CWT に時間信号全体の FFT が必要であるため，音響信号が逐次入力される場合には適用できない．このように場合には逐次処理が可能な CWT の高速計算法 [9] を用いることができ，本節では [9] と同一の方法で高速近似 CWT を拡張する．この手法をオンライン高速近似 CWT と呼ぶ．
オンライン高速近似 CWT は 3 ステップの処理からなる（図 1）．
（i）フレーム長 2N ，フレームシフ
ト N として，入力音響信号の当該時刻でのフレームを切り出す．
（ii）当該フレームの信号成分に分析窓
h = [h0 , · · · , h2N−1 ]⊤ を掛けたものに対して，高速近- 934 -

2016年3月

似 CWT を適用して当該フレームでのスペクトログラムを得る．
（iii）得られた各フレームでのスペクトログラム同士を対応する時刻で重畳することにより所望のスペクトログラムが得られる．
逆変換では，ステップ（ii）で得られた当該フレームのスペクトログラムと，当該フレームと直前のフレームのオーバーラップ部分の時間信号があればよい．当該フレームのスペクトログラムに逆高速近似 CWT を適用し得られた時間信号に合成窓 v = [v0 , · · · , v2N−1 ]⊤を掛け，適切にオーバーラップ加算することで入力された時間信号が得られる．ここで，合成窓は時刻イン
デックス n = 0, · · · , N −1 において hn vn +hn+N vn+N = 1を満たす．そのため，空間計算量の主な増大要因である各フレームのスペクトログラムの要素数は信号長に依存せず，信号長に比べ小さな N を用いれば空間計算量が削減できる．
ここで，hn として
0
 n − N−M

2 


0.5 − 0.5 cos πM−1

hn = 1
 n − 3N−M

2 

π
0.5+0.5cosM−1
0(0 ≤ n < N−M2 )
N+M
( N−M
2 ≤n< 2 )3N−M
( N+M2 ≤n<2 )( 3N−M≤ n < 3N+M22 )( 3N+M≤ n < 2N)2(5)で定義される Tukey 窓を用いると，n = 0, 1, · · · , N − 1において hn + hn+N = 1 が成り立つため各 n =0, · · · , 2N − 1 で vn = 1 となる．M (0 < M < N) はオー
バーラップ量を調節するパラメータであり， M が大きいほどフレーム同士のオーバーラップ量が多い．

3.2 実時間位相推定アルゴリズム
高速位相推定アルゴリズムは (2) 式から分かる通り位相全体へ定数を加えたものも解として許容するため，隣り合うフレーム同士で独立に高速位相推定アルゴリズムを適用した場合には少なくとも異なる初期位相をとりうる．そのため，フレーム毎に高速位相推定アルゴリズムを独立に適用すると，フレーム同士のオーバーラップ部分で不連続な信号成分となる位相が推定される可能性が高い．このような信号成分の不連続性を抑制するために，隣接フレーム間のオーバーラップ部分の信号成分が無矛盾であることを考慮したアルゴリズムを提案する．
直前のフレームで得られた時間信号の当該フレームに対応する部分を，N 次元の複素ベクトル g で表す．
g のうち直前のフレーム内の時刻に含まれない当該フレームの部分には 0 を詰める． s(a, ϕ) を当該フレームのスペクトログラムの推定値として再定義すると，
g が当該フレームの時間信号と無矛盾であれば，

W + s(a, ϕ) = diag(h)diag(v)(g + W + s(a, ϕ))(6)が成り立つ．そこで，(6) 式の右辺を s(a, ϕ) に対する日本音響学会講演論文集

時間信号への逆変換とみなし，(3) 式の代わりにs̃ ← Wdiag(h)diag(v)(g + W + s(a, ϕ))(7)を用いることで，フレーム同士のオーバーラップ部分で不連続な信号成分となることを抑制できると考えられる．本稿では，このアルゴリズムを実時間位相推定アルゴリズムと呼ぶ．
このアルゴリズムでは，当該フレームのスペクトログラムと g のみを保存すればよく，空間計算量はO(N + ΩD) である．そのため，空間計算量が信号長に依存せず，信号長が大きい場合に高速位相推定アルゴリズムに比べて空間計算量を格段に削減できる．

4 信号再構成実験による性能評価
4.1 実験条件
フレーム間の無矛盾性の考慮による効果と計算速度を評価するため，振幅スペクトログラムからの時間信号の再構成実験を行った．比較手法として，各フレームに対してそれぞれ高速位相推定アルゴリズムを独立に適用する方法（ベースライン法）を用いた．
実験データとして，RWC 音楽ジャンルデータベース[13] の 10 曲を用いた（サンプリング周波数 48 kHz）．
各曲の冒頭 30 s にオンライン高速 CWT を適用し得られた振幅スペクトログラムを入力とし，位相はランダムに初期化した．フレーム長は 170, 340, 680 ms（N = 212 , 213 , 214 ）とし，(5) 式で定義される分析窓を用いた（ M = N/4）．アナライジングウェーブレット
として対数正規分布型のウェーブレット [3] を用いた．
このウェーブレットの Fourier 変換は対数周波数領域で正規分布と同形であり，正規分布の標準偏差に対応するパラメータを 0.02 とした．フィルタの中心周波数が 25 cent 毎に 27.5 から 23679.5 Hz となるようスケールを設計し，高速近似 CWT での帯域制限の範囲は中心周波数から対数周波数上で [−2σ, 2σ] とした．
計算速度の評価指標として，フレームシフトに対する処理時間の比で定義される real time factor（RTF）を用いた．RTF が 1 以下であれば実時間で実行可能であり，低ければ低いほど高速である．再構成信号の音質の評価指標として，perceptual evaluation of audio
quality （PEAQ）[14] による objective diﬀerence grade（ODG）を用いた．ODG は −4 から 0 までの値をとり，ODG が大きいほど音質が高い．

4.2 結果図 2 に，提案法とベースライン法による再構成信
号の ODG と計算時間の結果を示す．提案法はベースライン法に比べ高い ODG で時間信号を再構成できており，フレーム間の無矛盾性が再構成信号の音質に重要であることが確認できる．実際に筆者がベースライン法の再構成信号を聴取したところ，フレームのオーバーラップ部分で音量が下がったように聴こえた．これは，オーバーラップ部分で互いに打ち消し合うような信号成分が得られたことが原因であると考- 935 -

2016年3月

えられる．短いフレーム長ほどベースライン法で得られた再構成信号の ODG が低い傾向があるが，これはフレーム長が短くなるほどフレーム数が増えオー
バーラップする信号成分が増えるためにフレーム間の無矛盾性が増しやすいからである．また，フレーム長 340, 680 ms とした提案法では実時間で ODG が−2 以上の再構成信号が得られており，実時間制約下でも実用上十分な音質の音響信号が再構成可能であった．提案アルゴリズムの収束性は理論的には保証されていないが，実験において各フレームでの目的関数は単調に減少したことを確認しており，収束性が保証される可能性が示唆された．

提案法
Objective difference grade-1

5 まとめ
本稿では振幅スペクトログラムからの実時間位相推定アルゴリズムを提案した．提案アルゴリズムでは固定長のフレーム毎に位相推定を行うことにより，空間計算量が信号長に依存せず長い信号に対しても位相を推定できる．隣接フレーム間でのオーバーラップ部分の信号成分に対する無矛盾性を考慮した更新式を用いることにより，各フレーム対して独立に高速位相推定を適用した場合よりも高音質な信号を振幅スペクトログラムから再構成できることを実験により確認した．また，音楽音響信号に対して実時間制約下でも実用上十分な音質の音響信号が提案アルゴリズムにより再構成できることも確認した．今後は，音声に関する性能評価実験や提案アルゴリズムの収束性に関する理論的な保証が課題である．

-1.5-2-2.5-3-3.5-400.2

0.4

0.6 0.8
1
Real time factor

1.2

1.4

1.2

1.4

1.2

1.4(a) フレーム長 170 ms
提案法

Baseline

Objective difference grade-1-1.5-2-2.5-3-3.5-40

謝辞 本研究は JSPS 科研費 26730100，15J0992 の助成を受けたものです．

0.2

0.4

0.6 0.8
1
Real time factor(b) フレーム長 340 ms
提案法

参考文献

Baseline

Objective difference grade-1[1] J. J. Burred et al., Proc. AES Convension, 2006.
[2] M. N. Schmidt et al., Proc. ICA-BSS, pp. 700–707, 2006.
[3] H. Kameoka, Ph.D. dissertation, The University ofTokyo, Mar. 2007.
[4] J. P. de León et al., Proc. DAFx, pp. 47–54, 2013.
[5] Y. Ikemiya et al., Proc. ICASSP, pp. 574–578, 2015.
[6] T. Irino et al., IEEE Trans. SP, 41 (12), pp. 3549–3554,
1993.
[7] T. Nakamura et al., Proc. DAFx, pp. 129–135, 2014.
[8] 亀岡他, “信号処理方法及び装置,” Nov. 2008, 特特開2008-281898.
[9] N. Holighaus et al., IEEE Trans. ASLP, 21 (4), pp. 775–785, Apr. 2013.
[10] C. Schörkhuber et al., Proc. AES Int. Conf. Semantic Audio, Jan. 2014.
[11] J. M. Ortega et al., Iterative solution of nonlinear equations in several variables, no. 30, 1970.
[12] J. Le Roux et al., Proc. ISCA Tutorial and ResearchWorkshops on Statistical and Perceptual Audition, pp.
23–28, 2008.
[13] M. Goto, Proc. Int. Congress Acoust., 1, pp. 553–556,
2004.
[14] “ITU-T recommendation BS.1387-1, Perceptual evaluation of audio quality (PEAQ): Method for objective measurements of perceived audio quality,” Sep. 2001.

日本音響学会講演論文集

Baseline-1.5-2-2.5-3-3.5-400.2

0.4

0.6 0.8
1
Real time factor(c) フレーム長 680 ms
Fig. 2: 提案法とベースライン法（Baseline）による様々なフレーム長での RTF に対する ODG の平均値と標準誤差．各点は，RTF が小さい方からそれぞれ各フレームで 0, 10, · · · , 200 反復だけアルゴリズムを実行した場合の結果である．

- 936 -

2016年3月