情報処理学会研究報告
IPSJ SIG Technical Report

無矛盾性規準に基づく連続ウェーブレット変換
スペクトログラムへの位相推定法と高速化
中村 友彦1,a)

亀岡 弘和1,2,b)

概要：人間の聴覚システムは対数的な周波数サブバンドを持つフィルタバンクとみなせると言われており，
連続ウェーブレット変換（Continuous wavelet transform, CWT）により同様のフィルタバンクを構成することができる．したがって，人間の聴覚システムの重要な機能を模倣した信号処理手法を開発するために，
CWT により得られる振幅スペクトログラム（振幅 CWT スペクトログラム）をモデル化，解析，処理するアプローチは有効であると考えられる．このアプローチを実現するためには，位相情報を含まない振幅
CWT スペクトログラムを時間領域の信号に変換する必要がある．これは，特に音響信号を出力するアプリ
ケーションにおいて重要である．そのため，本論文では音響信号を得るために，与えられた振幅 CWT スペ
クトログラムから位相を推定する高速アルゴリズムを提案する．定量評価実験により，従来のアルゴリズムに比べて提案アルゴリズムは，ほぼ同程度の高品質な音響信号を約 100 倍高速に得られることが示された．

1. はじめに連続ウェーブレット変換（Continuous wavelet transform,
CWT）は，対数周波数での時間周波数表現が可能な解析法である．人間の聴覚フィルタバンクは，特に高周波帯域に関し対数周波数上で一様な解像度を持つことが知られている [1]．このような人間の聴覚システムの重要な機能を模倣した信号処理を実現するには，CWT によって得られる
スペクトログラム（CWT スペクトログラム）を解析する手法が有望である．実際，近年の研究により振幅 CWT ス
ペクトログラム上での多重基本周波数推定などの多重音解
析手法が有効であることが示さている [2–5]．したがって，
振幅 CWT スペクトログラム上での音源分離や音響信号加工も有効であると考えられる．しかし，音源分離や音響信
号加工では，処理された振幅 CWT スペクトログラムを音響信号に変換する必要があり，そのためには与えられた振幅 CWT スペクトログラムに対し適切な位相を推定しなければならない．本論文では，この変換のために必要な位相を，与えられた振幅 CWT スペクトログラムから推定する手法を提案する．
振幅 CWT スペクトログラムからの位相推定アルゴリズ
ムは，Irino らによって提案された [6]．らのアルゴリズムでは，逆 CWT と CWT を行って得られた CWT スペクトログ
ラムの振幅部分（振幅 CWT スペクトログラム）を所望の振幅 CWT スペクトログラムと置換する処理を繰り返し行う．
1
2
a)
b)

東京大学 大学院情報理工学研究科,
東京都文京区本郷 7-3-1, 113-0033
NTT コミュニケーション科学基礎研究所,
神奈川県厚木市森の里若宮 3-1, 243-0198
nakamura@hil.t.u-tokyo.ac.jp
kameoka@hil.t.u-tokyo.ac.jp

c 1959 Information Processing Society of Japan⃝

ただし，CWT の計算時間は短時間 Fourier 変換（short-time
Fourier transform, STFT）に比べ非常に長いので，このアル
ゴリズムは長時間の計算を必要とする（5 節参照）
．そのため，実際のアプリケーションへの用いるには計算量の削減が欠かせない．
そこで，我々が以前開発した STFT での位相推定法 [7]から計算量削減のヒントを得る．我々は，STFT により得られたスペクトログラム（STFT スペクトログラム）の冗長性に着眼し，振幅 STFT スペクトログラムからの高速位相推定アルゴリズムを導いた．フレーム間にオーバーラッ
プする部分がある場合は，異なるフレーム間で音響信号の要素が共有されるため，得られた STFT スペクトログラム
は冗長性を持つ．この冗長性から，実際の音響信号から得られた STFT スペクトログラムが満たす条件を導き，その条件を基に STFT スペクトログラムらしさを測る規準（無矛盾性規準）を提案し，STFT スペクトログラムに対する位相推定問題を定式化した．この定式化に補助関数法と呼ばれる原理を適用して得られたアルゴリズムは，Griﬃn ら
による位相推定アルゴリズム [8] と一致することも示した．
この無矛盾性規準に基づく方法論は，アルゴリズムの高速化や収束性の証明にも有用であった．ここで，CWT スペク
トログラムも冗長性を持つことや，Griﬃn らによる STFT
での位相推定法 [8] と Irino らによる CWT での位相推定法
は統一的に記述できること [9] も鑑みると，この方法論はCWT にも有用だと考えられる．
我々はこの方法論に従い，CWT 振幅スペクトログラム
からの位相を推定する高速アルゴリズムを提案し，その収束性も示す．まず，実際の音響信号から得られた CWT ス
ペクトログラムが満たす条件（無矛盾条件）を導出し，その条件から CWT スペクトログラムらしさを測る規準とし

1

情報処理学会研究報告
IPSJ SIG Technical Report

W + s = argmin ∥s − W f̃∥22
逆CWT

無矛盾でない例

CWT
CWT
離散時間信号
の集合（
）

無矛盾な例

逆CWT

CWTスペクトログラムの集合（）

図 1 連続ウェーブレット変換によるスペクトログラムに対する無矛
盾条件．赤矢印と青矢印はそれぞれ無矛盾条件を満たす場合と満たさない場合を示している．

て無矛盾性規準を導き位相推定を最適化問題として定式化する．次に，補助関数法と呼ばれる原理を用いて最適化ア
ルゴリズムを導出し，従来法との関連を明らかにする．その後，CWT の高速近似アルゴリズム [10] を用いて高速な
アルゴリズムを与える．また，位相推定のためのスケールとマザーウェーブレットに関する要請も明らかにする．以下では，計算機への実装のため離散時間及び離散スケールドメインでの CWT と逆 CWT を考える．

2. 連続ウェーブレット変換スペクトログラムの無矛盾性
2.1 無矛盾条件
CWT のスケールはウェーブレットの局所的な周期に対応しており，本論文では対数周波数軸上で解析周波数の間隔が一様となるようなスケールの場合を考える．スケールと時間
シフトのインデックスをそれぞれ l ∈ [0, L−1] と t ∈ [0, T −1]とし，L 個の複素ベクトル sl := [sl,0 , sl,1 , · · · , sl,T −1 ]⊤ ∈ CT⊤⊤
LTを結合した複素ベクトル s := [s⊤を0 , s1 , · · · , sL−1 ] ∈ C
導入する．このとき，s が直流成分のない複素離散信号∑
f = [ f0 , f1 , · · · , fT −1 ]⊤ ∈ F , (F := {f ′ ∈ CT ; t ft′ = 0}) のCWT スペクトログラムとなるためには，CWT の基底波形を各行に並べた行列を W ∈ CLT ×T として
s = Wf W0  W  1 
W :=  .  , .. 


WL−1


 ψl,0ψ l,T −1
Wl :=  .
 ..

ψl,1

(1)
ψl,1
ψl,0

···
···
..
.

ψl,2

···


ψl,T −1 ψl,T −2 
..  .
. 

ψl,0

(2)

を満たす必要がある．ここで，Wl ∈ CT ×T は CWT スペク
トログラムの l 番目の成分に対応する CWT を表す行列であり，ψl,t := ψ(t∆/al )/al はスケール al でスケーリングされた時刻 t∆ のマザーウェーブレットを表す．∆ はサンプリン
グ周期，ψ(·) ∈ C は連続時間領域でのマザーウェーブレッ
トであり，アドミッシブル条件を満たす．一方で，逆 CWTは W の擬似逆行列 W + によって

f = W + s,

(4)

f̃∈F

W + := (W H W)−1 W H(3)と表現される．ここで，H は複素共役転置を表す．また，
W + は擬似逆行列であるため，

c 1959 Information Processing Society of Japan⃝をみたす．ここで，∥ · ∥2 は l2 ノルムである．
式 (1) から分かるように，CWT は T 次元複素ベクトル
から LT 次元の複素ベクトルへの変換であるため，CWT ス
ペクトログラムは冗長な表現である．したがって，図 1 に示すように CWT スペクトログラムの集合 W は LT 次元の複素ベクトルの集合 CLT の部分空間である．実際の信号の CWT スペクトログラムは逆 CWT と CWT を行った後にも元の CWT スペクトログラムと同一であるので，実際の信号の CWT スペクトログラムが満たす条件である無矛
盾条件は，0LT = s − WW + s(5)と定義できる．ここで，0LT は LT 次元の零ベクトルである．すなわち，W は (5) を満たす複素ベクトル s の集合
である．
また，W の各行を STFT の基底波形に置換すると，我々が以前提案した STFT スペクトログラムに対する無矛盾条
件と一致する [7].

2.2 無矛盾性条件と位相推定問題の関連
振幅 CWT スペクトログラムからの位相推定問題は，与えられた振幅 CWT スペクトログラムを持つような信号を再構成するための位相を推定する問題である．ここで，振幅 CWT スペクトログラムにランダムな位相を付与した複素ベクトルによって得られた信号は，図 1 の無矛盾でない
例のように，元の振幅 CWT スペクトログラムとは異なる可能性がある．したがって，再構成された信号の振幅 CWT
スペクトログラムが与えられた振幅 CWT スペクトログラ
ムと一致するためには，再構成された信号の CWT スペク
トログラムが無矛盾条件を満たさなくてはならない．そのため，位相推定に無矛盾条件が利用できる．
2.3 スケールとマザーウェーブレットへの要請
位相推定に無矛盾条件を用いる際に必要な条件の有無を調べるために，無矛盾条件について考察する．CWT は，
スケーリングされたマザーウェーブレットが各フィルタのインパルス応答に対応するフィルタバンクとして解釈できる．時間インデックスに対して (5) の両辺に離散 Fourier
変換（discrete Fourier transform, DFT）を適用すると，
0 = ŝ − Ŵ Ŵ + ŝ, Ŵ0  Ŵ1 
Ŵ =  .  , .. 


ŴL−1

Ŵl = FT Wl FT−1 ,(6)Ŵ + = (Ŵ H Ŵ)−1 Ŵ H

(7)が得られる．ここで，FT ∈ CT ×T は DFT 行列，ˆ は当該変数の DFT を表す．Wl は FT と FT−1 によって対角化され，対角化された行列 Ŵl の対角要素が l 番目のスケールに対応するフィルターの周波数特性を表す．この対角要素は，角周波数インデックスを k ∈ [0, T − 1] として2

情報処理学会研究報告
IPSJ SIG Technical Report

0 = ŝl,k −1 ∑ψ̂l,k ψ̂∗l′ ,k ŝl′ ,kC k l′と書ける．ここで，Ck は正規化係数，∗ は当該変数の複素共役を表す．
異なるフィルターのサブバンドが重複していれば(∃l , l′ , ψ̂l,k ψ̂l′ ,k , 0)，CWT スペクトログラムにとって(5) は非自明な制約である．しかし重複がなければ，振幅 CWT スペクトログラムに任意の位相を付与しても (5)が満たされるため，位相推定には無矛盾条件を用いることができない．したがって，位相推定のためにはスケールとマザーウェーブレットが適切に選択される必要がある．ス
ケールが適切に選択されていれば，Morlet や対数正規分
布型ウェーブレット [3]，auditory wavelet transform で用いられるウェーブレット [6] などが，適切なマザーウェーブ
レットとして使用できる．以下では，異なるフィルター同士の重複があるものと仮定する．
この CWT のサブバンドの重複に対する要請は，STFT
でのフレームの重複に対する要請に類似している．STFT
スペクトログラムに対する無矛盾条件は，重複したフレー
ムから得られたスペクトルが無矛盾でなくてはならないと言い換えることができる [7]．CWT の場合には，CWT ス
ペクトログラムの無矛盾条件が，重複したサブバンドを持つフィルターの出力が無矛盾でなくてはならないと解釈できる．

3. 無矛盾条件に基づく位相推定

逆CWT

CWT

離散時間信号の集合
図2

位相のみ更新CWTスペクトログラムの集合

位相推定の反復アルゴリズム．赤矢印と青矢印はそれぞれ(14) と (15) の更新ステップに対応している．

項に ϕ が含まれるため難しい．そこで，補助関数法 [11] と呼ばれる原理を用いて I(ϕ) を最小化する．補助関数法は，コスト関数 I(ϕ) が与えられたときに，補助変数 s̃ と補助関
数 I+ (ϕ, s̃) ≥ I(ϕ) を導入する．次に，I+ (ϕ, s̃) を ϕ と s̃ について交互に最小化することによって，間接的に I(ϕ) を単調減少させる方法である．特に，I(ϕ) は凸関数であるため局所最適解が大域最適化と同一であり，補助関数法に従えば反復的に I(ϕ) を最小とする解に漸近する．
擬似逆行列は (4) を満たすので，補助変数 s̃ ∈ W を導入して補助関数 I+ (ϕ, s̃) はI(ϕ) = min ∥s(a, ϕ) − W f̃∥22

(11)

= min
∥s(a, ϕ) − s′ ∥22′(12)≤ ∥s(a, ϕ) − s̃∥22 =: I+ (ϕ, s̃)

(13)

f̃∈F

3.1 位相推定問題の定式化
振幅 CWT スペクトログラムからの位相推定問題では，
まず振幅 CWT スペクトログラム a ∈ [0, ∞)LT が与えられる．しかし，この与えられた振幅 CWT スペクトログラム
はいかなる位相を付与しても無矛盾条件を満たすことができない可能性があるため，無矛盾条件を厳密に満たす位相を推定することは難しい場合もある．そこで我々は，位相推定問題を振幅 CWT スペクトログラムと推定された位相
ϕ ∈ [−π, π)LT によって得られた CWT スペクトログラムができるだけ無矛盾となるような位相を推定する問題として定式化する．この定式化では I(ϕ) ≥ 0 を ϕ について最小化する問題として，
I(ϕ) := ∥s(a, ϕ) − WW + s(a, ϕ)∥22 ,(9)

 jϕ0  e  ejϕ1 s(a, ϕ) := a ⊙  .  .. ejϕLT −1(10)と表す．ここで，s(a, ϕ) は推定 CWT スペクトログラムであり，⊙ は要素積を表す．導入した目的関数 I(ϕ) は，どの程度 s(a, ϕ) が無矛盾条件を満たすかを表現しており，I(ϕ)
が 0 に近づくほど s(a, ϕ) は近似的に無矛盾条件を満たす．

3.2 補助関数法による反復アルゴリズムの導出
I(ϕ) を ϕ について直接最小化することは，(9) の第 2c 1959 Information Processing Society of Japan⃝

所望のCWTスペクトログラムをもつ
複素ベクトルの集合

(8)

s ∈W

と定義できる．(13) の等号成立条件は s̃ = WW + s(a, ϕ) である．この等号成立条件と I+ (ϕ, s̃)/∂ϕ = 0LT の解から，更新式はs̃ ← WW + s(a, ϕ),(14)

ϕ ← ∠s̃(15)と導出できる．ここで，∠ は複素ベクトルの各要素の偏角を [−π, π)LT 上のベクトルとして与える．

3.3 従来法との関連
図 2 に示すように，(14) と (15) は s(a, ϕ) に逆 CWTと CWT を行うステップと，s(a, ϕ) の位相を s̃ と置換するステップを表している．このアルゴリズムは，Irino らによるアルゴリズム [6] と同一である．さらに，もし W の各行
を STFT の基底波形と置換すれば我々が以前導出した振幅
STFT スペクトログラムへの位相推定アルゴリズム [7] が得られる．したがって，STFT と CWT での位相推定アルゴ
リズムとその導出は統一的に記述できる．これは，Lopes
らによる考察 [9] と一致している．

4. 高速位相推定アルゴリズム
4.1 高速近似連続ウェーブレット変換
3 節で無矛盾条件に基づく反復アルゴリズムを導出した

3

情報処理学会研究報告
IPSJ SIG Technical Report

図3

スケーリングされたマザーウェーブレットの離散 Fourier 変換が，角周波数軸上で局所的にしか値を持たない例．マザーウェーブレットとして対数正規分布型ウェーブレット [3] を用
いた．

角周波数:

(i) ダウン
サンプリング
角周波数
インデックス(ii) 循環シフト

角周波数
インデックス
角周波数:
図4

高速近似連続ウェーブレット変換 [10] の概要．

番目のフィルターの出力の FFTŴl f̂ を k ∈ [B, B+D−1] の成分に帯域制限し，T 点から D 点にダウンサンプルする．このダ
ウンサンプルされた要素の角周波数は [2πB/D, 2πB/D + 2π]に分布しているので，(ii) ダウンサンプルされた要素を並び替えて角周波数を [0, 2π] に正規化する．信号の DFT は周期的であるため，角周波数の正規化を行うためにはダウ
ンサンプルされた要素を巡回的にシフトさせ，零位相に対応する部分が先頭となるようにすればよい．ここで，床関
数 ⌊·⌋ を用いて n := ⌊(B + D)/D⌋ となる整数 n を定義する．
(ii) の操作は，図 4 に示すように [2πB/D, 2πn] の要素をダ
ウンサンプルされた要素の後ろに連結すればよい．この処理の後，並び替えられた要素を逆 FFT することによってD 点にダウンサンプルされた CWT スペクトログラムの l
番目の成分を得られる．ここで，B, D は各フィルタ毎に異なっていてもよい．
ダウンサンプルによって削除された角周波数要素は，従来の CWT で 0 とみなして計算していることに相当する．
したがって，高速近似 CWT の逆変換（逆高速近似 CWT）
は，高速近似 CWT によって得られた CWT スペクトログラ
ムの各スケールに対応する成分に FFT を適用後，(ii) とは逆の操作の後ダウンサンプルによって削除された角周波数
要素を 0 詰めし，従来の逆 CWT と同様に計算すればよい．

4.2 反復アルゴリズムの高速化(i) と (ii) のステップは，]
0B ×(D−B0 )IB0  [ 0D×B ID 0D×(T −D−B)
K :=  0
ID−B00(D−B0 )×B0 |
{z
}
|
{z
}
(i) ダウンサンプリング(ii) 循環シフト

が，CWT の計算量は STFT の計算量よりも非常に大きい
ため，実際のアプリケーションで使用するためには，計算量を削減し計算速度を向上する必要がある．
計算量を削減するためには，我々が以前提案した高速
近似 CWT を利用できる [10]．2.3 節のフィルタバンク
としての CWT の解釈から，図 3 に示すように高速近似
CWT はフィルタの周波数特性が角周波数軸上で局所的にしか値を持たないことを利用する．このような局所化
は，Morlet や対数正規分布型ウェーブレット [3] をマザー
ウェーブレットとして用いたときに現れる．l 番目のフィルタに対し，周波数特性 {ψ̂l,k }k が角周波数インデックスk ∈ [B, B + D − 1] (0 ≤ B, 0 < D ≤ T ) に主要な値を持つとき，高速近似 CWT は k ∈ [B, B + D − 1] (0 ≤ B, 0 < D ≤ T )の各周波数領域の {ψ̂l,k }k のみを用いてダウンサンプルされた CWT スペクトログラムの l 番目の成分を計算する．ただし，D は 2 の冪乗でなくてはならない．また，高速近似
CWT で得られる CWT スペクトログラムはダウンサンプ
ルによって従来の CWT によって得られる CWT スペクト
ログラムに比べ時間解像度が低下しているが，適切に近似精度を調整すれば実用に十分な時間解像度で高速に計算できる．
高速近似 CWT は，音響信号 f とスケーリングされたウェー
ブレットの高速 Fourier 変換（fast Fourier transform, FFT）を計算後，各フィルターについて 2 ステップの処理を行う．(i) lc 1959 Information Processing Society of Japan⃝(16)
によって表現できる．ここで， B0 := B − (n − 1)D，ID は
D × D の単位行列，0D×B は D × B の零行列である．高速近似 CWT によって得られる CWT スペクトログラムの l 番
目の成分 šl ∈ CD は，行列 K を用いて
šl = F D−1 K Ŵl FT f

(17)と書ける．逆高速近似 CWT も高速近似 CWT の擬似逆行列によって表現できるため，3 節で得られた反復アルゴリ
ズムの CWT と逆 CWT の代わりに高速近似 CWT と逆高
速近似 CWT を用いれば，収束性を保証したまま反復アル
ゴリズムを高速化できる．

4.3 時間及び空間計算量の削減
信号長 T が十分大きければ，CWT と高速近似 CWT の時間計算量は，FFT の計算量に依存する．FFT は O(T log2 T ) の計算量をもつため，CWT の計算量は O(T log2 T + LT log2 T )で あ る ．一 方 ，高 速 近 似 CWT で l 番 目 の フ ィ ル タ で
Dl ≤ T にダウンサンプルした場合，計算量は O(T log2 T +∑L−1
l=0 Dl log2 Dl ) に削減される．フィルタの周波数特性が局所化していればいるほど，計算量の削減は大きい．
空間計算量についても，Irino らのアルゴリズム [6] に比べ提案法では削減されている．信号が十分長い場合には，
空間計算量の主な増大要因は CWT スペクトログラムの要4情報処理学会研究報告
IPSJ SIG Technical Reportと表される．ここで，ω は角周波数，σ は標準偏差を表す．
本実験では，σ を 0.02 とし，各フィルタの中心周波数が20 cent 毎に 27.5 から 7040 Hz となるようにスケールを設定した．提案法では，対数周波数軸で中心周波数から ±3σの領域を用いて計算を行った．Intel Xeon CPU E31245 (3.3
GHz)，32 GB RAM をもつ計算機を用いた．
音質の評価尺度として，音声の世界標準客観評価尺度である perceptual evaluation of speech quality (PESQ) [13] を用いた．PESQ は −0.5 から 4.5 までの値をとり，大きければ大きいほど音質が高いことを示している．計算時間の評価尺度としては，イテレーション毎の計算時間を用いた．
5.1.2 結果
PESQ の平均値と標準誤差は，Irino らのアルゴリズムは
4.20 ± 0.08，提案法は 4.1 ± 0.1 であり，得られた信号の音
声品質は実用上十分高いことが示された．また，両アルゴリズムの差異は実用上大差ないものであった*1 .
計算量は主に信号長に依存するため，図 5 に信号長に対するイテレーション毎の平均計算時間と標準誤差を示した．提案法は Irino らのアルゴリズムに比べ 100 倍の速度で計算できており，例えば Irino らのアルゴリズムでは 15s の信号に対し約 10 s/iteration で計算が行われるのに対し，提案法では 0.1 s/iteration であった．

5.2 実験 2: 近似精度と音質の関連
5.2.1 実験条件
提案法は近似を含むため，近似精度に対する再構成された信号の音質を評価した．RWC 音楽ジャンルデータベー
ス [14] の 102 個の 16 kHz サンプリングの音楽音響信号*1

サ ン プ ル 音 は http://hil.t.u-tokyo.ac.jp/˜nakamura/demo/fastCWT.html で試聴できる．

c 1959 Information Processing Society of Japan⃝

計算時間 [s/iteration]

5.1 実験 1: 計算時間と音質
5.1.1 実験条件
計算時間と位相推定により得られた信号の音質を評価するため，Irino らのアルゴリズム [6] と提案法の比較評価実験を行った．与えられる振幅 CWT スペクトログラムとして，ATR 日本語音声データベースの A セット [12] の女性
話者（faf）115 ファイルと男性話者（mht）の 113 ファイルの音響信号から得られた振幅 CWT スペクトログラムを用いた．FFT は 2 の冪乗の長さの信号に対して非常に高速に動作するので，使用した信号に対し信号長が最も近い 2 の冪乗となるまで零詰めを行った．位相はランダムに初期化
し，両アルゴリズムとも 1000 イテレーションまで計算した．マザーウェーブレットとして，対数正規分布型ウェー
ブレット [3] を用いた．このウェーブレットは Fourier 変換されたドメインで定義されており，)(
(log ω)2(ω > 0)
exp
−


4σ2
(18)
ψ̂(ω) := 0(ω ≥ 0)提案法[Irino1993]

100
10
1
0.1
0.01
0.001
0

2

4

6

8 10 12 14 16 18
信号長 [s]

図 5 信号長に対するイテレーション毎の平均計算時間と標準誤差．

[Irino1993]

Objective difference grade

5. 定量評価実験

1000

-0.5
-1
-1.5
-2
-2.5
-3
-3.5
-4
0

100

200
300
400
イテレーション回数500

図 6 様々な近似精度 ([−Pσ, Pσ] (P = 1, 2, 3, 5)) の提案法と Irino らのアルゴリズム [6] による，Perceptual evaluation of audio quality
で得られた objective diﬀerence grades のイテレーションに対す
る平均と標準誤差の時間発展．

[Irino1993]

Objective difference grade素数である．Irino らのアルゴリズムでは LT 個の要素が存∑在するが，提案法では l Dl 個にまで削減されている．

-0.5
-1
-1.5
-2
-2.5
-3
-3.5
-4
0

10

20
30
計算時間 [s]4050

図 7 様々な近似精度 ([−Pσ, Pσ] (P = 1, 2, 3, 5)) の提案法と Irino らのアルゴリズム [6] による，Perceptual evaluation of audio quality
で得られた objective diﬀerence grades の計算時間に対する平均と標準誤差の時間発展．

ファイルの 30 s から 35 s の部分の振幅 CWT スペクトロ
グラムを，データとして用いた．マザーウェーブレット
は σ = 0.02 の対数正規分布型ウェーブレットとした．提案法の近似精度はステップ (i) での帯域制限の範囲に依存する．この範囲を，中心周波数から対数周波数軸上で±Pσ (P = 1, 2, 3, 5) と変化させて計算した．提案法に関するイテレーション数は 500，Irino らのアルゴリズムでは5

情報処理学会研究報告
IPSJ SIG Technical Report

できることを示した．定量評価実験により，提案法が Irino
らのアルゴリズムに比べ約 100 倍高速であり，音楽と音声
音響信号に対して実用上十分な音質で信号を得られることを示した．
今後は，和音や調，音階の変換のための音楽音響信号の加工手法のポストプロセシングとして使用する予定である．また，実時間での処理のためにオンラインアルゴリズ
ムの構築も行う予定である．

Perceptual evaluationof speech quality

4.5
4
3.5
3
2.5
2
[Irino1993]

1.5
0
図8

20

40
60
計算時間 [s]80100様々な近似精度 ([−Pσ, Pσ] (P = 1, 2, 3, 5)) の提案法と Irino らのアルゴリズム [6] による，Perceptual evaluation of speech quality
スコアの計算時間に対する平均と標準誤差の時間発展．

謝辞
本研究の一部は，JSPS 科研費 26730100 の助成を受けた
ものです。
参考文献
[1]

100 とした．計算は，Intel Core i3-2120 CPU (3.30 GHz) と
8GB RAM を持つ計算機で行った．他の条件は，5.1.1 節と同一である．
音質評価尺度として，音質評価の世界的な標準尺度である perceptual evaluation of audio quality (PEAQ) [15] の objective diﬀerential grade (ODG) を用いた．この値は −4 から0 まで値をとり，ODG が大きくなればなるほど音質が高い
ことを示している．
5.2.2 結果
図 6 にイテレーションに対する平均 ODG と標準誤差を示した．P = 3, 5 での提案法の ODG は，100 イテレーショ
ン後には −2.0 より大きく，十分高音質である*2 ．音質について，これらの結果は Irino らのアルゴリズムによる結果と大きく異ならず，P ≥ 3 とすれば提案法は Irino らのアルゴリズムとほぼ同等の音質をもつ音響信号を得られると言える．
計算時間の削減のためには，P が小さければ小さいほど削減できる．図 7 は RWC-MDB-G-2001 No. 1 の音源に対する，計算時間と ODG の時間発展であり，提案法により得られた信号の ODG が従来法よりも急速に高くなったことが確認できる．同様の結果が，ATR 日本語音声データ
ベース [12] の A セットの fafsc110（7 s の音源）でも見られた（図 8）．我々は結果から，P = 3 程度の近似精度で提案法を用いれば，実用的な時間で高音質な信号が得られると結論づけた．

6. 結論
本論文では，収束性のある高速な位相推定アルゴリズムを提案した．位相推定問題を無矛盾条件をベースにして定式化し，補助関数法によって反復アルゴリズムを導出した．
また，このアルゴリズムが Irino らによって示されていたアルゴリズム [6] と同一であることも確認した．さらに，無矛盾条件から位相推定のためにスケールとマザーウェーブ
レットが満たすべき条件を導出した．高速近似 CWT [10]を用いて，反復アルゴリズムの時間及び空間計算量が削減
*2

[2]

[3]
[4]

[5]

[6]

[7]

[8]

[9]

[10]
[11]
[12]

[13]

[14]
[15]

入野 俊夫：[レクチャー講演] 聴覚フィルタの測定と定式化について，日本音響学会聴覚研究会資料，Vol. 39, No. 6,
pp. 413–418 (2009).
Schmidt, M. N. and Mørup, M.: Nonnegative matrix factor 2D deconvolution for blind single channel source separation,
Independent Component Analysis and Blind Signal Separation, Springer, pp. 700–707 (2006).
Kameoka, H.: Statistical Approach to Multipitch Analysis,
PhD Thesis, The University of Tokyo (2007).
Muller, M., Ellis, D. P. W., Klapuri, A. and Richard, G.: Signal processing for music analysis, IEEE J. Sel. Topics. Signal
Process., Vol. 5, No. 6, pp. 1088–1110 (2011).
de León, J. P., Beltrán, F. and Beltrán, J. R.: A complex
wavelet based fundamental frequency estimator in singlechannel polyphonic signals, Proc. Digital Audio Eﬀects(2013).
Irino, T. and Kawahara, H.: Signal reconstruction from modified auditory wavelet transform, IEEE Trans. Signal Process.,
Vol. 41, No. 12, pp. 3549–3554 (1993).
Le Roux, J., Kameoka, H., Ono, N. and Sagayama, S.: Fast
Signal Reconstruction from Magnitude STFT Spectrogram
Based on Spectrogram Consistency, Proc. Int. Conf. Digital
Audio Eﬀects, pp. 397–403 (2010).
Griﬃn, D. and Lim, J.: Signal estimation from modified
short-time Fourier transform, IEEE Trans. Acoust., Speech,
Signal Process., Vol. 32, No. 2, pp. 236–243 (1984).
Lopes, D. M. and White, P. R.: Signal reconstruction fromthe magnitude or phase of a generalised wavelet transform,
Proc. Eur. Signal Process. Conf., pp. 2029–2032 (2000).
亀岡 弘和，田原 鉄也，西本 卓也，嵯峨山 茂樹：信号処
理方法及び装置 (2008). 特開 2008-281898.
Ortega, J. M. and Rheinboldt, W. C.: Iterative solution of
nonlinear equations in several variables, No. 30 (2000).
Kurematsu, A., Takeda, K., Sagisaka, Y., Katagiri, S.,
Kuwabara, H. and Shikano, K.: ATR Japanese Speech
Database as a Tool of Speech Recognition and Synthesis,
Speech Commun., Vol. 9, No. 4, pp. 357–363 (1990).
ITU-T: Recommendation P.862, Perceptual Evaluation of
Speech Quality (PESQ): An Objective Method for End-ToEnd Speech Quality Assessment of Narrow-Band Telephone
Networks and Speech Codecs (2001).
Goto, M.: Development of the RWC Music Database, Proc.
Int. Congress Acoust., pp. l–553–556 (2004).
ITU-T: Recommendation BS.1387-1, Perceptual Evaluationof Audio Quality (PEAQ): Method for Objective measurements of perceived audio quality (2001).

c.f.) 使用した音楽音響信号を MPEG-3 に 160 kbps でエンコード
した場合は，平均 ODG と標準誤差は −3.68 ± 0.03 であった．

c 1959 Information Processing Society of Japan⃝

6