『Online Phase Reconstruction via DNN-based Phase Differences Estimation』（IEEE/ACM TASLP 2023）に基づく位相復元手法をPythonで実装した

はじめに

Masuyama氏らによる位相復元の論文が出版されている．

　Y. Masuyama, K. Yatabe, K. Nagatomo and Y. Oikawa, "Online Phase Reconstruction via DNN-Based Phase Differences Estimation," in IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 31, pp. 163-176, 2023.

論文は以下のリンク先から読める．

ただしソースコードは公開されていない．そこで本論文で導入されているニューラルネットワークおよび位相復元アルゴリズムをPythonで実装したので，それらを紹介するのが本記事の主旨である．JSUTコーパスのうちbasic5000およびonopatopee300を利用して位相復元実験を行い，著者らの提案手法について簡単な検証を行ったので，それらの結果について報告する．

本論文で提案されている手法については，以下のスライドの19ページから説明が始まるので，読者はまずこちらに目を通すのをオススメする． speakerdeck.com

手法

図1に振幅スペクトルと位相の差分特徴量を可視化したものを示す（上記スライド21ページより引用）．位相スペクトルを直接推定するのは難しいが，図に示されるようにその差分特徴量には振幅スペクトルと類似した明確なパターン（調波構造）が多く見つかるので，差分特徴量自体の推定は比較的易しいだろうという観察・考察に基づく．論文では位相の時間方向の差分特徴量をTPD，周波数方向の差分特徴量をFPDと呼ぶ．それぞれ瞬時周波数と（負の）群遅延に対応している．TPDはより推定が容易なBPDへと修正される．

図2に2段階のオンライン位相推定の図解を示す（論文のFig 1 より引用）．位相スペクトルの差分特徴量をDNNに基づいて推定し，それらから元の位相スペクトルを再構築する（図2 (i)）．BPDとFPDがそれぞれ対数振幅スペクトルを入力とする別々のDNNで推定される．BPDとFPDは続く第2段階で具体的な位相復元処理に利用される．図2の(ii)にはその詳細が図解されている．つまるところ，差分特徴量と辻褄の合う複素STFT係数を最小二乗誤差の観点で推定している．

ところで論文の表題にある"Online"とは，当該フレームの位相復元の際に，時間的に当該フレーム以前のいくつかの音声フレームたち，そして”ほんの少し”未来のフレームたちの参照を許容することを意味する（厳密にはそれら音声フレームから抽出される振幅スペクトルたち）．未来のフレームを参照しようとすれば，入力がやって来るまでの待機時間が発生するために処理全体が遅延する．未来方向の参照フレーム数が大きいほど，高い復元精度が期待できる一方で遅延もますます大きくなる．遅延自体はなるべく抑えたいので，復元精度を保ちつつ参照フレーム数をなるべく小さな値で留めたいわけである．この未来方向に参照するフレーム数が0の場合，特に「因果的」と称する．

"Online"とは対照的に，"Offline"ではオンライン性やリアルタイム性は一旦無視して，位相推定に際して一度に全フレームの振幅スペクトルが手に入る前提で処理を行う．参照できる音声フレームつまり振幅スペクトルが豊富なので，一般に"Online"と比較して高品質な復元音声が手に入るが，しばしば多くの計算時間が必要である．

位相の差分を推定するニューラルネットワーク

差分特徴量推定のためのDNNアーキテクチャを図3に示す（論文のFig 6 より引用）．何より特徴的なのは"1-D Frequency Convolution"（FreqConv）である．データのモダリティが音声のみの場合，1-D convと言えば時間方向にフィルタが動いていく畳込みをしばしば思い浮かべる．しかしながら，未来方向を参照する畳込みは今回使えない．FreqConvは時刻を固定したうえで周波数方向にフィルタが動く1次元畳込みであり，未来方向のフレームは常に参照されないので今回の定式化（＝オンライン）にとって都合が良い．

アーキテクチャ図に戻って，入力に一番近いFreqConvの最初のチャネル数がどのように決まるのだろうか．その点，論文にもきちんと書いてある．過去方向に参照する時間フレーム数は本論文では3であり，これに現在時刻のフレームを加えた計4フレームが時間的に隣接している．これら隣接4フレームをFreqConvのチャネル次元とみなす．論文のFig. 2を以下に拡大して示す．

濃い緑色が当該フレームを示しており，薄くなっているのが過去3フレームを示している．FreqConvは図の上下方向に進むのであり，「4次元の行ベクトルが上下方向に積み重なっている」とみなせば，FreqConvのチャネル数が4ということも容易に理解できる．このように最下層は4chであるが，上層のチャネル数はこれよりも大きい．

周波数方向のサイズはスペクトルの周波数ビンの総数に等しい．たとえばFFT長が512であれば総ビン数は257（= 512 / 2 + 1）である．今回は周波数方向にダウンサンプリングするプーリング層が用意されていないので，下層から上層まで常に固定サイズの1次元畳込みが必要である．そしてこれを毎フレーム行うので，計算コストは決して小さくない．

訓練時にFreqConvをフレーム順に直列で行うと計算時間がとても延びてしまうので，全フレームを並列処理できるようにうまく実装を工夫しなければならない．オンライン性が特に求められるのはあくまで訓練済みモデルを用いた推論時であって，訓練時は積極的に並列処理を導入する．詳細は「実装」のセクションで述べる．

FreqConv以外の層であるFreqGatedConvおよびそれらの残差接続は図に示してある通りで，分かりやすい構成である．FreqConvさえ理解できれば実装の困難さはないだろう．差分特徴量もまた$2\pi$の周期性を有しているため，その損失関数にはvon Mises DNNやRPUと同じくcosine lossを利用している．先行研究で実績のある損失関数の採用は手堅いといえる．

位相復元アルゴリズム

解説スライドから復元アルゴリズムの核を説明している箇所（25ページ）を引用し，図5に示す．

推定されたBPD（から変換されたTPD）とFPDに基づき，フレーム毎に位相を復元する．位相差分は複素指数関数の観点では「比」になっているので，推定するべき複素STFT係数もまた隣接するSTFT係数との「比」を最小二乗誤差の意味で最適化するものとなっている．複素平面上でSTFT係数の最適化問題を解いているため，位相差分に存在する2$\pi$の不定性を回避している．

最終的には複素係数の線形方程式系を解くことに帰着される．この方程式系は係数行列がいわゆる帯行列（疎行列; sparse）であるため，密行列（dense）を解く場合と比べて大幅に計算コストが削減される．

実装

以下のリポジトリに置いた．Enjoy!

github.com

ファイル名	機能
config.yaml	各種の設定ファイル
preprocess.py	各種の前処理を実施するスクリプト
dataset.py	データローダの構築を実施するスクリプト
model.py	ネットワークを定義
factory.py	オプティマイザ，スケジューラ，loss関数を定義
training.py	モデルの訓練を実施するスクリプト
evaluate_scores.py	音質の客観評価指標を計算するスクリプト

evaluate_scores.py 内の compute_1st_stage および compute_2nd_stage 関数に位相復元手法が実装されている．

compute_1st_stage 関数ではDNNを駆動して対数振幅スペクトルからBPD（＝TPD），FPDを推定する．
compute_2nd_stage 関数ではBPDおよびFPDに基づいて構築される線形方程式系を frame-by-frame で繰り返し解いて，位相スペクトルを推定する．

既存のオンライン手法との比較のため、以下のスクリプトを用意した：

ファイル名	機能
evaluate_scores_spsi.py	Single Pass Spectrogram Inversion (SPSI) に基づく位相復元と客観評価指標の計算
evaluate_scores_rtpghi.py	Real-Time Phase Gradient Heap Integration (RTPGHI) およびPGHIに基づく位相復元と客観評価指標の計算
evaluate_scores_rtsisla.py	Real-Time Iterative Spectrogram Inversion with Look Ahead (RTISILA) に基づく位相復元と客観評価指標の計算

config.yamlで先読みフレーム数を0にすることで，オンライン設定が可能であり，さもなくばオフライン設定となる．RTPGHIはオンライン，PGHIはオフライン手法である．同様にRTISILAは先読みフレーム数の設定でオンライン・オフラインを切り替える．これら既存手法の実装は前回記事で紹介したLTFATおよびPHASERETによった．Oct2PyによりPythonから呼び出している．

oct2py経由でPHASERETを使って位相復元を行うためのPythonスクリプトを書いた話 - 備忘録

訓練時におけるテンソルの工夫

すでに述べた通り，訓練時は全ての音声フレームが与えられたオフラインな状況なので，最初の入力テンソルのshapeは

$\begin{align*} (B, T, K) \end{align*}$

である．ここで $B$ はミニバッチサイズ， $T$ は音声フレーム数， $K$ は周波数ビンのサイズである．本記事の実験では音声データを0.9秒ずつ複数のセグメントに分割しているので， $T$は約110である．また $K$ は 257である．このテンソルはFreqConvの入力に適合するshapeではなく，訓練の並列処理もできないのでうまく変形する必要が出てくる．そこでPyTorchのunfold関数をミニバッチ構成直前のcollate関数内で適用することで，shapeが

$\begin{align*} (B, T-3, K, 4) \end{align*}$

へと変わる．この4は過去フレーム3+現在フレームの1から来ている．4はFreqConvのチャネル数，$K$ がその畳込みに要するサイズなので，これをtransposeして

$\begin{align*} (B, T-3, 4, K) \end{align*}$

となる．さらにaxis 0 と axis 1 をflattenすれば

$\begin{align*} (B(T-3), 4, K) \end{align*}$

となる．このように大きなバッチ次元として改めてまとめてしまえば，時間方向の並列訓練が可能となる．このテンソルをFreqConvのConv1dモジュールに入力すれば良い（今回の実装はPyTorchなので）．テンソルのshapeは上層へ向かって，

$\begin{align*} (B(T-3), 4, K) \rightarrow (B(T-3), 64, K) \rightarrow \ldots \rightarrow (B(T-3), 1, K) \end{align*}$

と推移する．そして最上層でsqueezeしたのちにreshapeすれば，最終的なshapeは

$\begin{align*} (B, T-3, K) \end{align*}$

へと変わる．これが訓練時におけるBPDないしFPDのshapeを与える．unfold関数の副作用として時間フレーム数が小さくなるので，ロス計算時に正解として与えるBPDやFPDのフレーム対応には注意する．推論時には上記の込み入ったreshapeは必要ない．

補足：unfold関数の挙動

unfold関数に不慣れな読者のために動作例を示しておく．簡単のため，1次元のテンソルで考える．

>>> import torch
>>> x = torch.arange(0, 10)
>>> x
tensor([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> x.unfold(0, 4, 1) # axis=0 に対して size=4 のスライスを作る スライスのstepは1
tensor([[0, 1, 2, 3],
        [1, 2, 3, 4],
        [2, 3, 4, 5],
        [3, 4, 5, 6],
        [4, 5, 6, 7],
        [5, 6, 7, 8],
        [6, 7, 8, 9]])

例えば最上段の [0, 1, 2, 3] は4番目のフレーム '3' に対して過去3フレームが付随した形とみなせる．実装ではunfoldを適用するaxisを音声フレームに対応するところ，すなわちaxis=1に指定している．

位相復元の線形方程式を解くための工夫

先述の通り，係数行列は疎行列であり，具体的にはバンド幅1の帯行列となる．scipyには帯行列の線形方程式系を効率的に解くための solve_banded 関数（Link）があり，これを利用できるが，今回の実装では使わなかった．

せっかく行列が疎なのだから，同じくscipyのsparseパッケージを使って効率よく解いてみたいと思ったわけである．係数行列自体の計算もsparse arrayを用いて効率よく行える．また方程式を解くためにspsolve関数を用いた．

性能評価スクリプトの所要時間を短縮するための工夫

性能評価スクリプトの中では評価用データとなる300個の音声サンプル（後述）に対して位相復元を適用し，客観評価スコアを計算している．位相復元自体は音声サンプルごと独立に行われ，互いに干渉しないのは自明なので，積極的に並列化すべきと考える．今回はDNNによる推論部分を並列化した．

単一GPU上にDNNモデルを置いて推論を並列で行わせるためには，torch.multiprocessingの set_start_method を事前に呼び出しておく必要がある．

from torch.multiprocessing import set_start_method
...
set_start_method("spawn")

その後，Pythonの通常の並列処理の枠組みで推論を実装すれば良いので，難しくはない．今回はconcurrent.futuresによった．

評価実験

JSUTコーパスのbasic5000を訓練データに，onomatopee300を評価データに用いて位相復元実験を行った．客観評価指標にESTOI，(wideband) PESQ，LSC (log-spectral convergence)を採用した．

提案手法と比較するのは以下のオンライン手法たちである：

Single Pass Spectrogram Inversion (SPSI)
Real-Time Phase Gradient Heap Integration (RTPGHI)
Real-Time Iterative Spectrogram Inversion (RTISI)

論文でも言及されているが，本実験は完全にフェアな比較ではない．提案手法のDNNは訓練データの事前情報をパラメタの中に取り込んでいるが，比較手法はそれら事前情報を明示的に利用していないからである．

実験条件

計算機環境を示す．

計算機環境	バージョンなど
OS	Ubuntu 22.04
CPU	Intel i9-9900K
GPU	RTX2070
Python	3.10.12
PyTorch	2.2

ネットワークアーキテクチャの条件は以下の通りであり，論文に沿っている．FreqConv層は入力側と出力側で1つずつなので層数は2である．FreqConv層のカーネルサイズは1であるが，いわゆる1x1 convによるチャネル数変換に相当するとみなせばよい．パラメータ数は約206kであり，これまで実装してきたDNN系位相復元のネットワークよりも圧倒的に小さい．

項目	設定
FreqConv層の数	1 + 1
FreqGatedConv層の層数	5
FreqGatedConv層のチャネル数	64
FreqGatedConv層のカーネルサイズ	5
FreqConv層のカーネルサイズ	1
過去フレームの参照数	3
パラメータ総数	205,825

訓練の設定は以下の通りである．オプティマイザはRAdamだが，weight_decay はデフォルトでAdamと同じく勾配にdecayがかかり、これは実質的にL2正則化となるのがPyTorchの公式実装である．decoupled_weight_decay オプションを有効化（True）すると AdamW的な重み減衰として適用される．本実験ではdecoupled_weight_decay を True にした．

項目	設定
ミニバッチサイズ	32
エポック数	15
オプティマイザ	RAdam
学習率	0.001
勾配クリッピングのしきい値	10.0
重み減衰の強さ	0.0001
学習率スケジューラ	linear warmup つき cosine annlealing （半周期）
warm up のエポック数	5
warm up 開始時の学習率	0.000001
annealing の終端での学習率	0.0005

音声フレーム特徴量の設定は以下の通りである．標本化周波数は16000Hzにリサンプリングしている（オリジナルは48000Hz）．

項目	設定
標本化周波数	16000Hz
音響特徴量	対数振幅スペクトルおよび位相スペクトル
分析窓	ハン窓
FFTの窓長	512
フレーム長	512
フレームシフト	128

今回はbasic5000の各音声クリップを0.9秒ごとに事前に分割した．分割により生じた0.9秒未満の端数は訓練に利用しない．最終的には19027個の音声クリップに分割された．端数を全て合計しても1分程度であり，訓練に使用せずともその影響は実質的に無視できる．

提案手法については位相復元の「重み行列」の値を制御するハイパパラメタを設定している．論文ではべき乗パラメタ$p$および係数$\gamma_0$である．本実験はそれぞれ2.5と500に設定した．

実験結果

図6に各評価指標の箱ひげ図を示す．提案手法は"TOPR"（Two-stage Online Phare Reconstruction）で示されている．図より，提案手法はLSC以外の2指標で他のオンライン手法を上回ったことが示された．LSCに関してTOPRがSPSIとほぼ同レベルに留まった点やPESQが4.0程度に留まった点は論文の傾向とは大きく異なる点であり，後の考察セクションで詳しく述べたい．

おまけ

提案手法のオフライン版も実装してあるので，信号処理ベースのオフライン手法PGHIとRTISI_LAとで比較してみる．ちなみに提案手法のオフライン版といっても先読みフレーム数を3だけ追加したに過ぎず，入力側FreqConvのチャネル数が7（= 4 + 3）まで増えている点以外はオンライン版と全く変わらない．

客観評価指標の箱ひげ図を図7に示す．

図より，他のオフライン手法はスコアが改善されていることが分かる．RTISI_LA は改善が著しい．しかしながら提案手法はオフライン版に変更してもオンライン版とスコアがほぼ変わらなかった．そこで実装が悪かったのかを検証するために以下の「おまけ」実験を追加で実施した．

おまけ2

RPUおよびwRPU（重みつきRPU）でも差分特徴量を推定するDNNと位相復元アルゴリズムの組み合わせだったが，ここでは提案手法のDNNによって推定された差分特徴量をRPU側でも共通して用いる．よってRPUと提案手法との間で差異は位相復元アルゴリズムのみであり，よりフェアな比較ができる．wRPUを比較に含めたのは，提案手法にも「重み行列」が導入されているからである．

オンライン版DNNで得られた客観評価指標の箱ひげ図を図8に示す．

図より，提案手法はRPUからの音質改善が十分に得られており，位相復元アルゴリズムはRPUから着実に改善されていることが分かる．wRPUはRPUから大きく改善しているのはもちろんだが，TOPRのスコアも上回った．

続いてオフライン版DNNで得られた客観評価指標の箱ひげ図を図9に示す．

上図の結果より，DNNをオフライン設定にすることでRPUのスコアが改善した様子がよく分かる．wRPUは改善が頭打ちになっている．RPUの結果はオフライン版DNNによる差分特徴量の推定精度の改善，ひいては位相スペクトルの推定精度の改善によるものだと考えられる．したがってオフライン設定であってもDNN自体は期待通り動作しており，その実装は首尾よく完成したといえるだろう．したがってまた「おまけ1」や「おまけ2」でTOPRのスコアが頭打ちになっているのはDNNによる推定精度の問題ではなく，位相復元アルゴリズム自体の特性によるものかもしれない．考えられる他の原因については続く「考察」セクションで述べる．

考察

論文では提案手法のPESQスコアが4.4近くも出ているが，今回の実験では4.0程度に留まった．スコアに開きが出た原因を考えてみると，まずは訓練データ量の違いである．論文のLJ speech datasetは合計24時間の英語単一話者の音声データセットだが，basic5000は合計5時間程度の日本語単一話者の音声データセットである．ネットワークのパラメータ数は論文と揃えているが訓練データ量が大きく異なるので，basic5000上で過学習の傾向が強まったのかもしれない．最終的な復元音声の音質が訓練データ量が多いほど高いという仮説は分かりやすいが，それほど多くの訓練用音声を集められない場合にはデータ拡張（data augmentation）を図るほうが良いのかもしれない．もしくは重み減衰をより強くかけるなどのチューニングもまた効果的かもしれない．
次に考えられる原因として，訓練誤差を十分に収束させることができなかったかもしれない，という点である．訓練が不十分であるために，差分特徴量の推定誤差もまた大きくなり，評価指標の値も伸び悩んだというわけである．学習率やエポック数，スケジューリングをもう少し入念にチューニングする必要があったのかもしれない．チューニングには評価用データonomatopee300以外のJSUTコーパスのサブセット，例えばvoiceactress100が使えるだろう．評価用データ上で直接チューニングするのは好ましくないので．
ほかにスコアの開きが出た原因として音声の分析条件の違いが挙げられる．論文の実験ではLJ speech datasetはサンプリング周波数が22050Hzであり，リサンプリングせずそのまま訓練と推論に用いられたが，本実験ではbasic5000オリジナルのサンプリング周波数48000Hzから16000Hzにリサンプリングしている．分析窓長およびシフト長も論文では1024と256であったが本実験ではそれぞれ512と128であった．これら分析条件の違いが最終的な復元音声の音質にどこまで影響するかは未検証である．最初から条件揃えて実験しろよ，と言われればその通りであるが，いずれにせよPESQ計算時には波形データを16000Hzへとリサンプリングする必要はあるので，本実験では特徴抽出の段階からリサンプリングしたのであった．本記事の実験結果にもスコアのリファレンスとして一定の価値があると期待する．
提案手法のハイパパラメタ（重み行列にかかるべき乗指数など）もまたスコアに大きく影響を与えることが分かっている．きちんと探索を行えば，よりよいスコアが得られるようなパラメタの組が見つかるだろうし，Optunaなどのブラックボックス最適化のツールを援用するのも1つの手だろう．しかしながら，ある程度の範囲を人力で探しても，PESQは3.9〜4.0の範囲がほとんどであったし，ESTOIも高々0.98程度だった．
提案手法はオンライン手法ではあるが計算コストが大きい．複素係数の線形方程式を frame-by-frame で何度も解かねばならず，他のオンライン手法と比べても計算時間が延びがちである．アプリケーションに依存するが，実装する際にはさらなる工夫が必要となるかもしれない．
RPUの実験結果から，差分特徴量の推定精度が向上しても復元音声の音質が向上しなかったという点，またすべての実験を通じて提案手法のLSCは論文と異なり全然振るわなかった点にも言及したい．1st stageにおけるDNNの損失関数にはLSC由来のものは含まれておらずcosine lossのみなので，これらに一応の説明はつけられる．つまりDNNは差分特徴量のcosine lossを最小化するように訓練されているが，最終的な復元音声の評価指標，特にLSCを最小化するようには訓練されていない．よってLSCを含めた最終的な音質を評価できる損失関数を新たに導入してDNNをfine-tuningし，さらなる音質改善を図ることは研究の1つの方向性になるだろう．ほか2nd stageの最適化にも音質にかかる制約は含まれていない．位相再構築にかかる全体の処理に通じた損失関数ないし制約条件を新たに導入し，最適化問題を定式化し直すこともまた興味深い方向性である．

おわりに

本記事ではDNNに基づくオンラインの位相復元手法を実装し，その性能をJSUTコーパスを使った客観評価実験により確認した．実際，論文ほどの突き抜けた音質は達成できなかったものの，信号処理ベースのオンライン手法を上回る性能を達成でき，おおむね期待通りの結果が得られたと言えるだろう．その一方で論文で報告されているスコアとの傾向の食い違いも散見されるので，さらなる検証は必要だろう．論文には他にも多くの実験結果が報告されているが，すべてを再現するのは大変なので，この記事の読者に任せよう（他力本願）．

追記

LJ speech dataset上で行った実験結果は続く記事に書いておいた。

tam5917.hatenablog.com