JP7218688B2

JP7218688B2 - 位相推定装置、位相推定方法、およびプログラム

Info

Publication number: JP7218688B2
Application number: JP2019135981A
Authority: JP
Inventors: 悠馬小泉; 義紀升山; 浩平矢田部
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 2019-07-24
Filing date: 2019-07-24
Publication date: 2023-02-07
Anticipated expiration: 2039-07-24
Also published as: JP2021018398A

Description

この発明は、振幅スペクトログラムから位相スペクトログラムを推定する技術に関する。

音声合成や音声強調などの多くの音響信号処理のアプリケーションは、時間領域の観測信号（波形）を短時間フーリエ変換（STFT: Short-Time Fourier Transform）などを利用して時間周波数領域に変換して処理を行う。この際、信号は複素数になるが、複素数の取り扱いの難しさから、多くのアプリケーションは振幅スペクトル（複素数の大きさ）のみを制御し、位相スペクトル（複素数の角度）は制御しない。振幅のみを制御すると、複素スペクトログラムの振幅と位相の一貫性（consistency）が崩れてしまうため、逆変換した際に時間周波数表現（複素スペクトログラム）と対応した時間領域信号（波形）が存在せず、出力音の音質が低下するという問題がある。そのため、振幅スペクトログラムから位相スペクトログラムを推定する手法が必要である。

この問題を定義する。なお、文中で使用する記号「~」「^」は、本来直前の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直前に記載する。数式中においてはこれらの記号は本来の位置、すなわち文字の真上に記述している。例えば、「^φ」は数式中では次式で表される。

まず、時間領域のN点の信号をx=[x₀, x₁, …, x_N-1]^T∈R^Nとする。この信号を窓関数g=[g₀, g₁, …, g_L-1]^T∈R^Lを利用して短時間フーリエ変換した結果は、式（１）のように記述できる。

ここでj=√-1、δはシフト長、τ（=0, 1, …, T-1）は時間のインデックス、ω（=0, 1, …, K-1）は周波数のインデックスである。振幅スペクトルとは|X_ω,τ|のことであり、位相スペクトルとはΦ_ω,τ=Arg(X_ω,τ)∈[-π,π)である。ただし、Arg(X)はXの角度である。また以降、振幅スペクトログラムとは振幅スペクトルを周波数×時間に並べた行列、位相スペクトログラムとは位相スペクトルを周波数×時間に並べた行列、複素スペクトログラムとは複素スペクトルを周波数×時間に並べた行列とする。

この問題を解くために様々な研究がされているが、近年、深層学習を用いて振幅スペクトログラムから直接位相スペクトログラムを推定する手法が提案されている。この方法では、位相スペクトログラムの各要素である位相スペクトルを式（２）で推定する。

ここで、Fはディープニューラルネットワーク（DNN: Deep Neural Network）、θはそのパラメータ、Ψは入力する特徴量である。位相は循環変数であるため、ディープニューラルネットワークの学習でよく利用される最小二乗誤差の最小化によるθの学習は適切ではない。非特許文献１では、循環変数の確率分布の一つであるフォン・ミーゼス分布の尤度最大化に基づき、式（３）のコスト関数を最小化するように学習を行っている。

ただし、Φ^* _ω,τは真の位相である。

S. Takamichi, Y. Saito, N. Takamune, D. Kitamura, and H.Saruwatari, "Phase reconstruction from amplitude spectrograms based on von-Mises-distribution deep neural network," In International Workshop on Acoustic Signal Enhancement (IWAENC), pp. 286-290, September 2018.

短時間フーリエ変換では、データがサンプルシフトすると、振幅スペクトログラムはほとんど変化しないにも関わらず、位相スペクトログラムは大きく変化してしまうという問題がある。そのため、振幅スペクトログラムから位相スペクトログラムを一意に推定することはできない。ゆえに式（２）のように、振幅スペクトログラムから位相スペクトログラムを直接推定することは困難である。

図１にサンプルシフトによる振幅スペクトログラムと位相スペクトログラムの変化の例を示す。図１の一段目は音声を短時間フーリエ変換して得た振幅スペクトログラム、二段目は位相スペクトログラム、三段目は瞬時周波数である。図１の左側は元の音声の分析結果であり、右側は波形を時間方向にサンプルシフトしたものの分析結果である。図１によれば、振幅スペクトログラムはほとんど変化していないにも関わらず、位相スペクトログラム（特に、四角形で囲われた箇所）は大きく変化していることがわかる。

この発明の目的は、上記のような技術的課題に鑑みて、振幅スペクトログラムから位相スペクトログラムを推定することである。

上記の課題を解決するために、この発明の一態様の位相推定装置は、音声信号の音響特徴量から瞬時周波数を推定する瞬時周波数推定部と、音声信号の音響特徴量から群遅延を推定する群遅延推定部と、直前の時刻の位相スペクトルおよび瞬時周波数と現在時刻の群遅延とを用いて音声信号の位相スペクトログラムを推定する位相推定部と、を含む。

この発明の位相推定技術によれば、振幅スペクトログラムから位相スペクトログラムを推定することができる。

図１はサンプルシフトによる振幅と位相の変化を例示する図である。図２は位相推定装置の機能構成を例示する図である。図３は位相推定方法の処理手順を例示する図である。図４は発明の実行結果を例示する図である。図５はコンピュータの機能構成を例示する図である。

［発明の概要］
位相スペクトログラムの時間方向の微分（瞬時周波数）と周波数方向の微分（群遅延）は、対数振幅スペクトログラムと関係があることが知られている（参考文献１－３参照）。

〔参考文献１〕F. Auger, E. Chassande-Mottin, and P. Flandrin, "On phase-magnitude relationships in the short-time Fourier transform," IEEE Signal Processing Letters, Vol. 19, No. 5, pp. 267-270, May 2012.
〔参考文献２〕Z. Prusa, P. Balazs, and P. L. Sndergaard, "A noniterative method for reconstruction of phase from STFT magnitude," IEEE/ACM Transactions on Audio, Speech, Language Processing, Vol. 25, No. 5, pp. 1154-1164, May 2017.
〔参考文献３〕K. Yatabe, Y. Masuyama, T. Kusano, and Y. Oikawa, "Representation of complex spectrogram via phase conversion," Acoustic Science & Technology, Vol. 40, No. 3, pp. 170-177, May 2019.

図１に示すように、サンプルシフトした場合に、位相スペクトログラムは大きく変化しているにもかかわらず、瞬時周波数はほとんど変化していない。そこで、本発明では、位相の微分値（すなわち、瞬時周波数と群遅延）をディープニューラルネットワークで推定し、そこから位相を逐次推定する。

瞬時周波数と群遅延から位相スペクトログラムを推定する手法として、参考文献４がある。

〔参考文献４〕D. C. Ghiglia and M. D. Pritt, "Two-dimensional phase unwrapping: theory, algorithms, and software,"Wiley-Interscience, 1998.

ここで、参考文献４に記載された画像における位相復元を、本願発明の課題である音響信号における位相スペクトログラムの推定へ応用することを考える。V∈R^K×T-1を瞬時周波数、U∈R^K-1×Tを群遅延とおいたとき、瞬時周波数と群遅延からの位相スペクトログラムの推定を、式（４）の最小二乗問題として解くことができる。

ここでΦ∈R^K×Tは、Φ_ω,τを並べた行列、||・||_Fはフロベニウスノルム、D_τとD_ωはそれぞれ時間方向と周波数方向の差分作用素である。しかし、式（４）を解くためには、KT×KTの行列の逆行列を求める必要があり、多大な計算コストを要する。

そこで、Φのτ列目のベクトルであるτ番目の時間フレームの位相スペクトル^φ_τを式（５）で再帰的に求める。

ここで、~φ_τ-1=W(^φ_τ-1)はτ-1番目の時間フレームの位相スペクトル、||・||₂はL2ノルムである。ただし、W(φ)は式（６）で定義される。

ここで、

は床関数である。また、v_τとu_τはVとUのτ列目のベクトルであり、ディープニューラルネットワークを用いて推定する。つまり、v_τは式（７）であり、u_τは式（８）である。ψ_τはτ番目の時間フレームの音響特徴量である。

式（５）の解は式（９）で求められる。

ここで、Iは単位行列、D^*はDの随伴作用素である。この計算はK×Kの行列演算のみで求めることができるため、式（４）よりも高速かつ逐次的に推定することができる。

［変形例］
式（９）では、瞬時周波数と群遅延の推定精度が、全ての時間周波数ビンで等しいと仮定した。しかし、振幅が小さい時間周波数ビンでは推定精度が悪い可能性がある。そこで、式（９）をσ_ω,τ∈(0, 1)で重みづけした式（10）で求める方法をとることもできる。

ただし、||x||² _Σ=x^TΣxであり、Σ_IFは対角要素がσ_τ=[σ_0,τ, σ_1,τ, …, σ_K-1,τ]^Tで与えられる対角行列であり、Σ_GD=I-Σ_IFである。式（10）は重み付き最小二乗問題のため、式（９）と同様、解析的に解くことができる。σ_τ=[σ_0,τ, σ_1,τ, …, σ_K-1,τ]^Tは、例えば、式（11）のようにディープニューラルネットワークで推定すればよい。

［学習の実装例］
本発明で利用するディープニューラルネットワークは、F_IFとF_GDの２つである。また変形例としてF_weightを利用することもできる。この学習は何を利用してもよいが、一例として以下の方法を挙げる。

まず、F_IFとF_GDは式（12）（13）のコスト関数の最小化で学習できる。

ここで、F_IF(ψ_τ,θ_IF)_ωとF_GD(ψ_τ,θ_GD)_ωはv_τとu_τのω番目の要素である^V_ω,τと^U_ω,τを意味する。またV^* _ω,τとU^* _ω,τは数値微分で推定された瞬時周波数と群遅延の真値である。

F_weightの学習は、例えば式（14）の一貫性に基づくコスト関数の最小化で学習できる。なお、一貫性については参考文献５，６を参照されたい。

ただし、f⁺は逆STFTであり、^X(θ_weight)のω,τ要素は式（15）で求められる。

ここで、^φ(θ_weight)_ω,τは式（10）で求められる位相の推定値である。

〔参考文献５〕S. O Arik, H. Jun, and G. Diamos, "Fast spectrogram inversion using multi-head convolutional neural networks," IEEE Signal Processing Letters, Vol. 26, No. 1, pp. 94-98, January 2019.
〔参考文献６〕J. Le Roux, N. Ono, and S. Sagayama, "Explicit consistency constraints for STFT spectrograms and their application to phase reconstruction," In ISCA Workshop on Statistical and Perceptual Audition (SAPA), pp. 23-28, September 2008.

以下、この発明の実施の形態について詳細に説明する。なお、図面中において同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。

［実施形態］
実施形態の位相推定装置は、音声信号の振幅スペクトログラムから位相スペクトログラムを推定して出力する信号処理装置である。実施形態の位相推定装置１は、図２に示すように、瞬時周波数推定部１１、群遅延推定部１２、および位相推定部１３を備える。この位相推定装置１が図３に例示する各ステップの処理を行うことにより実施形態の位相推定方法が実現される。

位相推定装置１は、例えば、中央演算処理装置（CPU: Central Processing Unit）、主記憶装置（RAM: Random Access Memory）などを有する公知又は専用のコンピュータに特別なプログラムが読み込まれて構成された特別な装置である。位相推定装置１は、例えば、中央演算処理装置の制御のもとで各処理を実行する。位相推定装置１に入力されたデータや各処理で得られたデータは、例えば、主記憶装置に格納され、主記憶装置に格納されたデータは必要に応じて中央演算処理装置へ読み出されて他の処理に利用される。位相推定装置１は、少なくとも一部が集積回路等のハードウェアによって構成されていてもよい。

以下、図３を参照して、実施形態の位相推定装置１が実行する位相推定方法の処理手続きを説明する。

位相推定装置１には、入力音響信号の各フレームの音響特徴量ψ_ω,0, ψ_ω,1, …, ψ_ω,T-1が入力される。以下では、τ番目の時間フレームの位相スペクトルφ_ω,τを推定する処理を説明する。τ-1番目の時間フレームの音響特徴量ψ_ω,τ-1は瞬時周波数推定部１１へ入力される。τ番目の時間フレームの音響特徴量ψ_ω,τは群遅延推定部１２へ入力される。

ステップＳ１１において、瞬時周波数推定部１１は、入力された音響特徴量ψ_ω,τ-1から学習済みのディープニューラルネットワークF_IFを用いてτ-1番目の時間フレームの瞬時周波数v_ω,τ-1を推定する。具体的には、瞬時周波数推定部１１は、上記式（７）を計算する。瞬時周波数推定部１１は、推定した瞬時周波数v_ω,τ-1を位相推定部１３へ出力する。

ステップＳ１２において、群遅延推定部１２は、入力された音響特徴量ψ_ω,τから学習済みのディープニューラルネットワークF_GDを用いてτ番目の時間フレームの群遅延u_ω,τを推定する。具体的には、群遅延推定部１２は、上記式（８）を計算する。群遅延推定部１２は、推定した群遅延u_ω,τを位相推定部１３へ出力する。

ステップＳ１３において、位相推定部１３は、瞬時周波数v_ω,τ-1と群遅延u_ω,τとτ-1番目の時間フレームについて推定した位相スペクトル~φ_ω,τ-1とを用いて式（５）を解くことでτ番目の時間フレームの位相スペクトル^φ_ω,τを求める。具体的には、位相推定部１３は、式（９）を計算して位相スペクトル^φ_ω,τを得る。変形例に示したように、位相推定部１３は、式（10）を解くことで位相スペクトル^φ_ω,τを求めてもよい。位相推定部１３は、推定した位相スペクトル^φ_ω,τを出力する。同時に、位相推定部１３は、次の時間フレームの位相スペクトル^φ_ω,τ+1を求めるときのために、位相スペクトル^φ_ω,τを図示しない記憶部等に記憶する。

［実験結果］
図４に実施形態の位相推定装置による実行結果の例を示す。図４の一段目はスペクトログラム、二段目は位相スペクトログラム、三段目は瞬時周波数、四段目は群遅延である。図４の二～四段目の左側は真値であり、右側は推定値である。図４の結果から、瞬時周波数と群遅延は、推定値が真値とほぼ一致しており、位相スペクトログラムも精度よく推定できていることがわかる。

以上、この発明の実施の形態について説明したが、具体的な構成は、これらの実施の形態に限られるものではなく、この発明の趣旨を逸脱しない範囲で適宜設計の変更等があっても、この発明に含まれることはいうまでもない。実施の形態において説明した各種の処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。

［プログラム、記録媒体］
上記実施形態で説明した各装置における各種の処理機能をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムを図５に示すコンピュータの記憶部１０２０に読み込ませ、制御部１０１０、入力部１０３０、出力部１０４０などに動作させることにより、上記各装置における各種の処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。

また、このプログラムの流通は、例えば、そのプログラムを記録したDVD、CD-ROM等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるASP（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

１位相推定装置
１１瞬時周波数推定部
１２群遅延推定部
１３位相推定部

Claims

推定すべき位相スペクトログラムをΦとして、Φをτ番目の時間フレームに関する情報ごとに分割した位相情報＾φ _τ を再帰的に求めることでΦを決定する位相推定装置であって、
音声信号のτ－１番目の時間フレームの音響特徴量から瞬時周波数ｖ _τ－１をディープニューラルネットワークで推定する瞬時周波数推定部と、
上記音声信号のτ番目の時間フレームの音響特徴量から群遅延ｕ _τをディープニューラルネットワークで推定する群遅延推定部と、
＾φ _τ と＾φ _τ－１から求めた瞬時周波数とｖ _τ－１の差分、および＾φ _τ から求めた群遅延とｕ _τ の差分を最小化する＾φ _τ を決定する位相推定部
を含む位相推定装置。
音声信号の音響特徴量から瞬時周波数を推定する瞬時周波数推定部と、
上記音声信号の音響特徴量から群遅延を推定する群遅延推定部と、
直前の時刻の位相スペクトルおよび上記瞬時周波数と、現在時刻の上記群遅延とを用いて上記音声信号の位相スペクトログラムを推定する位相推定部と、を含み、
上記位相推定部は、~φ_τ-1を時刻τ-1の位相スペクトルとし、v_τ-1を時刻τ-1の瞬時周波数とし、u_τを時刻τの群遅延とし、D_ωを周波数方向の差分作用素とし、||・||₂をL2ノルムとし、次式を解くことで上記音声信号の位相スペクトル^φ_τを求める、

位相推定装置。
音声信号の音響特徴量から瞬時周波数を推定する瞬時周波数推定部と、
上記音声信号の音響特徴量から群遅延を推定する群遅延推定部と、
直前の時刻の位相スペクトルおよび上記瞬時周波数と、現在時刻の上記群遅延とを用いて上記音声信号の位相スペクトログラムを推定する位相推定部と、を含み、
上記位相推定部は、~φ_τ-1を時刻τ-1の位相スペクトルとし、v_τ-1を時刻τ-1の瞬時周波数とし、u_τを時刻τの群遅延とし、D_ωを周波数方向の差分作用素とし、Σ_IFを対角要素が予め学習した重みσ_τ=[σ_0,τ, σ_1,τ, …, σ_K-1,τ]^Tで与えられる対角行列とし、Iを単位行列とし、Σ_GD=I-Σ_IFとし、||・||² _Σ=・^TΣ・とし、次式を解くことで上記音声信号の位相スペクトル^φ_τを求める、

位相推定装置。
推定すべき位相スペクトログラムをΦとして、Φをτ番目の時間フレームに関する情報ごとに分割した位相情報＾φ _τ を再帰的に求めることでΦを決定する位相推定方法であって、
瞬時周波数推定部が、音声信号のτ－１番目の時間フレームの音響特徴量から瞬時周波数ｖ _τ－１をディープニューラルネットワークで推定し、
群遅延推定部が、上記音声信号のτ番目の時間フレームの音響特徴量から群遅延ｕ _τをディープニューラルネットワークで推定し、
位相推定部が、＾φ _τ と＾φ _τ－１から求めた瞬時周波数とｖ _τ－１の差分、および＾φ _τ から求めた群遅延とｕ _τ の差分を最小化する＾φ _τ を決定する、
位相推定方法。
請求項１から３のいずれかに記載の位相推定装置としてコンピュータを機能させるためのプログラム。