JP4469986B2 - 音響信号分析方法および音響信号合成方法 - Google Patents

Description

本発明は、音声や楽器音などの音響信号を高精度で分析可能な方法および分析された音響信号を所望の目的のために目標とする形態に再合成可能な方法に関する。

音声の有声部や多くの楽器音は、その周波数応答が線スペクトル構造になることが知られている。例えば、有声音声の発話においては、声帯の振動周期や声道の形状は時間と共に変化する。これらの変化は、生成される音声の基本周波数(F0)や振幅包絡の変化として現れる。したがって、有声音声や楽器音を精度良く音響分析するためには、スペクトル全体形状の時間変化パタンを正確に分析する技術が必要である。

この分野における技術開発の進展を考えてみるに、Fant(1960 年) は、声帯振動により生成された音源信号が、声道を通過する際にその形状に応じて変調され、多様な有声音声が生成されるという音声生成のソース・フィルタ理論を提案した。声帯振動による音源信号は、基本周波数（F0）の整数倍の周波数にエネルギーを持つ調波構造を持ち、F0は声帯振動の速さに応じて時々刻々変化する。また、各調波成分が受ける変調の強さは、声道形状の特性を反映して周波数ごとに異なっており、かつ、声道の時間変化に伴って変化する。

このような特性を持つ有声音信号は、瞬時振幅や周波数が滑らかに変化する正弦波（sinusoid）の和として表すことができる。この考え方はSinusoidal Modelingと呼ばれ、現在まで活発な研究が続いている（例えば、非特許文献１）。このSinusoidal Modelingにおいては、有声音信号は次式（数１）のように表される。ここでx(t)は有声音信号、Kはsinusoidの数を表し、ξk(t),ψk(t), ηkはｋ番目のsinusoidの瞬時振幅、瞬時周波数、初期位相にそれぞれ対応する。

上記数１から明らかなように、もし全てのsinusoidのパラメータ（ξk(t),ψk(t), ηk）が与えられれば、そこから有声音信号x(t)を得ることは容易である。しかし、この逆変換、即ち有声音信号x(t)自体から、各sinusoid成分のパラメータを推定する過程は一種の不良設定問題になっているため、何らかの拘束条件を設けなければ解くことができない。

上述の非特許文献１(1986)で、McAulayとQuatieriは、有声音声を一定の時間間隔で周波数分析し、それぞれの分析時刻近傍ではsinusoid成分の瞬時振幅と周波数が定常とみなせるという拘束条件を設定して、パラメータ推定を行った。まず各時刻のスペクトルから調波成分に対応するピーク振幅と周波数を検出し、次に時間的連続性に基づいてそれらを結びつけることで、sinusoid成分の振幅と周波数の軌跡を計算する。彼らは、ひとつの自然発話音声を例に取り、このアイディアが上手く働くことを示しているが、そのシステムには多数のヒューリスティクスが含まれているため、有声音声全般、例えば話者の性が異なる場合にも同等の結果が得られる保証はない。

また、分析時間近傍で振幅と周波数が定常的であるという拘束条件は、これらのパラメータの時間変化が大きい時刻での推定精度の劣化を招く。この問題を解決するために、Hermus(2005)は、ガウス型の窓関数を用いてスペクトルを計算するExponential Sinusoidal Modeling (ESM)を提案した（非特許文献２）。彼はESMのパラメータ推定性能を、聴覚マスキングパターンに基づく精度の指標であるNMR(Noise to Masking Ratio)を用いて評価している。ESMを用いることにより確かにNMRは向上するが、推定されたパラメータから再合成した信号は、元の音声とは大きく異なる。実際、彼が行った聴取実験では、再合成信号と原音は容易に弁別可能（つまり、再合成した信号は、元の音声とは大きく異なっている）という結果が示されている。

Sinusoidal Modelingの長所のひとつとして、信号に含まれるsinusoid成分のパラメータが推定されれば、そこから波形自体を再合成するのが容易であるという点が挙げられる。従って、入力信号と再合成信号の差を計算すれば、パラメータ推定性能を定量的に評価することが可能である。しかしながら、現在までのところ、この様な厳密な指標に基づいてsinusoidal modelingの性能評価を行った研究は報告されていない。
McAulay, R.J. and Quatieri, T.F., "Speech Analysis/Synthesis Based on a Sinusoidal Representation", IEEE trans. on Speech and Signal Processing, ASSP-34(4), p744-754, 1986 年 Hermus, K., "Perceptual audio modeling with exponentially damped sinusoids", Signal Processing 85, p163-176, 2005 年

そこで本発明の課題は、時間変化を伴うsinusoid信号のパラメータ推定精度を大幅に向上でき、それによって、音響分析、さらには、音響再合成の精度を大幅に向上できる、新規な音響分析／合成技術を提供することにある。より具体的には、本明細書にて局所変化率符号化（Local Vector Coding：LVC）と呼ぶ新しい手法を提案し、高精度の音響分析技術を提供し、その分析に基づいて自由に音響を再合成可能な技術を提供することを課題とする。

上記課題を解決するために、本発明に係る音響信号分析方法は、時間変化特性を持つ入力音響信号に対して、局所変化率符号化により第一調波成分の位相関数を推定した後、信号の時間軸をこの位相軸に変換し、変換信号を再び局所変化率符号化により分析することで、入力信号の全成分の瞬時振幅、瞬時周波数および瞬時位相を、再合成可能なパラメータ関数として出力することを特徴とする方法からなる。分析対象となる音響信号としては、音声信号は勿論のこと、楽器音信号も含まれる。すなわち、より正確に言えば、詳細を後述するように、時間変化特性を持つ入力音響信号に対して、入力信号の正弦波成分の瞬時振幅、瞬時周波数、瞬時位相を下記の４段階の計算処理に基づいて、再合成可能なパラメータ関数として出力することを特徴とする音響信号分析方法である。
（１）入力信号〔数１〕（前述）の任意の分析時刻において、正弦波パラメータを〔数２、３〕（後述）のように近似し、パワースペクトルと位相スペクトルから、〔数１２〜１４〕（後述する。特許請求の範囲における〔数４〜６〕に対応している。）のようにこれらのパラメータ候補を定める局所変化率変換処理を行う段階
（２）各分析時刻で得られた瞬時振幅と位相のパラメータ候補を、動的計画法を用いた時間連続性によりグルーピング化し、これらのパラメータを連続関数として出力するパラメータ軌跡推定処理を行う段階
（３）入力信号〔数１〕に含まれる成分間の干渉により、上記（１）と（２）の段階では十分なパラメータ推定精度が得られない場合に、入力信号の基本位相のみを推定した後で、入力信号の時間軸を基本位相軸に変換した信号を生成することで成分間の干渉を抑制する時間−位相変換処理を行う段階
（４）上記（３）の段階で得られた変換信号に、（１）と（２）の段階を適用して正弦波パラメータを推定し、得られた連続関数の基本位相軸を時間軸に逆変換することで入力信号の瞬時振幅と位相を連続関数として出力する時間−位相逆変換処理を行う段階

また、本発明に係る音響信号合成方法は、上記のような音響信号分析方法における各パラメータ関数を用いて音響信号を再合成することを特徴とする方法からなる。再合成は、忠実に元の音響信号を再現する場合は勿論のこと、再合成された音響信号から、または、音響信号を再合成するに際し、音質を他の音質に変換することも可能である（例えば、男性音声から女性音声へのモーフィング等も可能である。

本発明に係る音響信号分析方法および音響信号合成方法によれば、従来技術では得られなかった精度で、音響信号を正確に分析でき、かつ、分析された音響信号を、所望の目的のために目標とする形態に自由に再合成することができる。

以下に、本発明の実施の形態について、本発明の完成に至るまでの開発経過とともに、図面を参照しながら詳細に説明する。
まず、本発明を完成するための研究においては、時間変化を伴うsinusoid信号のパラメータ推定精度を、再合成信号と入力信号の差である残差信号と、入力信号のエネルギー比（S/R: Singnal per Residual ratio）を評価指標とし、sinusoidal modelingの性能を評価することを検討した。以下に、有声音声の分析を主体に説明するが、本発明は楽器音の分析にも適用可能である。

（１）基本周波数(F0)検出（Pitch Determination）
有声音声の基本周波数(F0)は、第一調波成分の周波数に対応し、知覚される声の高さ（ピッチ）を決定付ける主要因である。また、F0は、話者の情動表出や個人性の識別等の手がかりであることが知られており、音声の音響分析において応用性の高い特徴量であると言える。通常sinusoidal modelingでは、有声音信号に含まれる全成分のパラメータを推定する。従って、第一調波成分のパラメータであるF0も不可避的に推定される。本発明では、sinusoidal ModelingのF0検出手法としての側面にも注目する。

音声信号からF0を検出する手法は、pitch determination algorithm(PDA)と呼ばれ、これまで広範囲の研究がなされている。これらを大まかに分類すると、次の３つのグループになる（Hess, 1983）。第一は周波数領域の手法であり、これは音声スペクトルに現れる調波構造を利用してF0を検出する。この手法の代表的なものはケプストラム法である（Noll, 1966）。第二は、時間領域の手法である。これは有声音声の準周期性を利用してF0を検出するもので、自己相関関数法などが代表的である（Sondhi, 1968）。第三は時間−周波数領域の手法で、上の二つの組み合わせに対応する。この手法では、帯域通過フィルタを用いて音声を複数の周波数チャネルに分割し、チャネルごとの信号の周期性を評価することでF0を検出する。

初期のPDAの研究では、入力として単一でクリーンな音声が用いられた。近年は、複数の音声や雑音下でのPDAに関する研究が報告されている（Shimamura and Kobayashi, 2001; Wu, et. al, 2003）。自然発話音声のPDAの性能を定量的に評価するためには、検出誤差を計算するためのリファレンスとなるF0が必要となる。これら最近の研究でも、雑音や他の音声と混合させる前に、クリーンな単一音声に対して、ケプストラムや自己相関を用いてリファレンスF0を計算している。しかし、この計算においては、真のF0の整数倍の周波数が検出されるdouble pitch問題や、F0の変化速度が速い時刻での値のジャンプ等がしばしば観測されるため、実験者が手動でそれらの値を修正する作業が必要になっている（Wu, etl al., 2003）。従って、Hess(1983)がその当時に指摘した事実−”例えクリーンで単一の音声であっても、そこから自動的にF0を検出する手法は未だ完成されていない”−は、それから二十年以上を経た現在においても真実であると言わざるを得ない。

自然発話音声のF0検出を困難にする要因のひとつは、sinusoidal modelingの場合と同様、信号パラメータの時間変化である。F0の時間変化が大きい音声では、スペクトルの調波構造が歪むため、分析時間近傍で定常性を仮定するケプストラムや自己相関関数法では、正確なF0を検出することができない。これに対して、本発明で提案する局所変化率符号化（LVC）は、時間変化を伴う信号を正確に分析できる点が特徴であり、新たなF0検出法としての有効性が期待できる手法である。以下に、既存手法と比較実験を行い、LVCのF0検出性能、その手法の有効性について定量的に評価する。

（２）時間変化（以下、単に「時変」と言うこともある。）sinsuoid信号の性質
有声音声のsinusoid modelingにおいては、振幅や位相などのパラメータの時間変化の処理が重要な問題になってくる。ここでは、まず瞬時振幅と位相が単純な関数で時間変化するsinusoid信号を例に取り、その特性を理論的に説明する。

（２−１） Quadratic-Parameter Sinusoid信号
有声音声に含まれるsinusoid成分のパラメータは、複雑な時間変化パタンを示す。このような成分のひとつに注目し、その瞬時対数振幅A(t)とP(t)を時刻tc近傍でテイラー展開し、二次までの項で打ち切ると次式数２、数３が得られる。

ここでa₀,a₁,a₂,p₀,p₁,p₂は定数であり、a₀が瞬時対数振幅、a₁が振幅の変化速度、a₂が振幅の変化加速度、p₀が瞬時位相、p₁が角周波数、p₂が角周波数の変化速度にそれぞれ対応する。例えば、a₁=a₂=p₂=0の場合は角周波数p₁の純音を意味する。以下では簡単のため、時刻tc =0の場合を考える。これらのパラメータから生成される複素信号s(t)は次式数４で表される。

尚、α₀,α₁,α₂は複素定数である。この複素信号の実部が、成分信号の時刻c近傍の近似に対応する。自然発話された音声の瞬時振幅や位相は、二次関数よりも複雑なパタンを有するが、最も単純な時間変化信号という意味で、この信号の特性を調べることは有効である。以下ではs(t)をQuadratic-Parameter Sinusoid(QPS)と呼ぶ。

（２−２） QPSの周波数応答
QPS信号s(t)の周波数応答S₀(ω)を次式数５、数６で計算する。ここでTは分析時間を表す定数、w(t)はガウス型の窓関数、γは窓関数の時間長を決める正の定数である。

QPSのパラメータのあらゆる範囲に対して、数５を解析的に求めることは不可能である。だが、少なくとも分析時刻の端点ｔ＝±Tにおいて、s(t)*w(t)の振幅包絡がゼロとみなせる場合には、数５を解析的に解くことができる。この条件は、具体的には振幅包絡が、t=-T〜Tの間でピーク値Apeakを持つことと、その値が分析時間の端点の値Aedge=A(±T)に対してMergin(dB)だけ大きいことに対応する。これを満たすためには、パラメータa₁とa₂が次式数７、数８の範囲に存在する必要がある。

上記条件が満たされれば、QPSの周波数応答S₀(ω)は次式数９のように書ける。即ち、S₀(ω)の対数振幅応答と位相応答は、角周波数ωの二次関数で表される。

図１にS₀(ω)の例を示す。図１(a) は、QPSの周波数振幅応答を示している。実線は瞬時対数振幅と瞬時周波数が時間変化するQPSの応答を表す（a₀=a₂=p₀=0, a₁=0.4dB/ms, p₁=100Hz, p₂=1.0Hz/ms）。薄線は純音に対応するQPSの応答を表す（a₀= a₁=a₂=p0= p₂=0, p₁=100Hz）。矢印は、振幅応答のピーク周波数と、実際の瞬時周波数との差（F-shift）に対応する。 (b)は、QPSの周波数位相応答。実線と薄線は(a)と同じ。応答のピーク周波数における位相と、実際の瞬時周波数との差（P-shift）を矢印で表す。 (c)は、パラメータp₂とa₁を変化させた場合のF-shiftの変化を示す。それぞれの線がa₁=0.2, 0.4, 0.8 dB/msのQPSに対応する。 (d)は、(c)と同様にパラメータを変化させた場合のP-shiftを表す。図１(a),(c) において薄線が純音の周波数応答（a₀=a₁=a₂=p₀=p₂=0, p₁=100 Hz）に、また実線は振幅と周波数が線形に変化するQPSの応答(a0 =a2=p0 =0, a1= 0.4 dB/ms, p₁=100 Hz p₂=1.0 Hz/ms)にそれぞれ対応する。これら二つの信号はともに時刻t=0で瞬時周波数と位相がゼロ、瞬時周波数が100Hzになるが、変化成分を含むQPSの振幅応答のピーク周波数は100Hzにはならず、ピークにおける位相もゼロにはならない。この応答ピーク周波数と位相のずれを図１(b),(d) に示す。ここで示したずれの大きさは、数９から一意に導くことが可能である。

この結果は、例えば周波数応答のピーク検出のような手法では、時間変化を含むQPSの瞬時周波数や振幅を正確に分析できないことを示唆する。また、調波構造を持つ有声音声では、各成分の瞬時周波数がF0の整数倍であっても、F0や振幅が時間変化する場合には、スペクトル上の振幅ピークが必ずしもF0の整数倍の周波数に存在しないことを説明できる。

（２−３）局所変化率変換（Local Vector Transform: LVT）
QPSの周波数応答の一般式数９から、その角周波数に対する一次微分S₁(ω)と二次微分S₂(ω)は、次式数１０、数１１のように計算できる。

もし入力信号s(t)から、これら3種類の周波数応答S₀（ω), S₁(ω), S₂(ω)を得ることができれば、任意の角周波数ωにおいて信号パラメータα₀,α₁,α₂を一意に決定可能である。これは次式数１２、数１３、数１４（特許請求の範囲における数４、数５、数６に対応）のように表される。

図２に、有声音声入力に対して、推定されたパラメータを示す。入力は男性話者の/iyoiyo/という発話の最初の母音/i/の部分に対応する。音声データはATRディジタル音声データベースの音素バランス216単語のM107B-0002を用いている。図２(a)においては、実線は入力信号の振幅スペクトルを表す。入力は男性話者の発話/iyoiyo/の最初の母音/i/の時間に対応する。薄線はLocal Vector Transform(LVT)により推定された信号パラメータから求めた振幅応答を表す。推定パラメータは、140Hzの整数倍の周波数応答から計算し、第一成分から第五成分までを図示している。 (b)では、実線は入力音声の位相応答を示す。薄線は、LVTにより推定されたパラメータから算出された位相応答を表す。 (c)は、LVTにより推定された瞬時対数振幅（パラメータa₀）を表す。入力音のエネルギーが十分な領域では、ほぼ同一の値が推定されている。 (d)は、LVTにより推定された瞬時位相（p₀）を表す。 (e)は、推定された瞬時対数振幅の変化速度（a₁）を表す。(f) は、推定された瞬時周波数（p₁）を表す。(g) は、推定された瞬時対数振幅の変化加速度（a₂）を表す。(h) は、推定された瞬時周波数の変化速度（p₂）を表す。図２(a),(b) が入力の振幅応答と位相応答を表し、図２(c) 〜(h) が数１２、数１３、数１４を用いて推定されたパラメータを表す。また、周波数が140Hzの整数倍の推定パラメータを用いて、数９の周波数応答を再合成し、図２(a),(b) の薄線に示している。周波数140Hzの整数倍の近傍では、入力の周波数応答と再合成された周波数応答が、振幅、位相とも良く一致していることが確認できる。これは、自然発話音声のような実際の音声に対しても、QPSに基づくパラメータ推定が有効であることを意味する。本明細書では、このパラメータ推定の演算をLocal Vector Transform(LVT)と呼ぶ。

（３）局所変化率に基づく音声符号化システム(LVC)
上述した局所変化率変換(LVT)に基づいて、有声音声のsinusouidal modelingシステムを構築する。入力となる音声信号から、各成分の瞬時対数振幅と位相を推定する処理を、局所変化率符号化(LVC)と呼ぶ。以下ではLVCシステムの構成について述べる。

（３−１）概要
まず、sinusoidal modelingの一般式数１を次の数１５のように変形する。ここでAk(t)はｋ番目の成分の瞬時対数振幅であり、Pk(t)は瞬時位相に対応する。ちなみに成分kの瞬時角周波数はPk(t)の時間微分により計算できる。

図３に、LVCシステムの計算ブロックを示す。ここで、x(t)：入力信号, P_F0(t) ：第一調波成分の瞬時位相, Ai(t)：第i調波成分の瞬時振幅, Pi(t)：第i調波成分の瞬時位相, y(t)：再合成波形、をそれぞれ表している。システムの入力は有声音声信号x(t)であり、システムの出力は各成分の瞬時対数振幅関数Ak(t)と瞬時位相関数Pk(t)である。一度これらのパラメータ関数が推定できれば、信号波形y(t)は容易に再合成できる。

このパラメータ関数の推定は、二段階の計算により実行する。第一段階目の計算では、入力信号の第一成分の瞬時位相 P_F0(t) を推定する。まず入力信号の各時間フレームにおいて、前述のLVTにより成分パラメータを計算する。図２に示した通り、LVTは任意の時刻の周波数応答からパラメータを推定するが、システム出力である振幅関数や位相関数を得るためには、各時間フレームで計算されたパラメータ値を時間方向に接続する必要がある。この処理は、次のTrajectory Estimationのモジュールで計算する。この時点でx(t)の第一成分の瞬時位相関数PF0(t)が得られる。F0検出タスクにおいては、この P_F0(t) の時間微分を２πで除した値を出力とする。

第二段階目の計算では、まず，入力信号x(t)の時間軸を、第一段階で推定された P_F0(t) ＝φによって位相に変換する（Time to Phase Conversion, 後述する。）。高次の調波成分においては、周波数変化速度が大きいため、ひとつの成分のエネルギーピークが他の成分の周波数領域まで広がってしまう。時間−位相変換は、この成分間干渉を低減するために導入する。変換後の波形ｘ（φ）は、F0の時間変化が正規化され、成分間干渉が小さくなる。

この変換波形ｘ（φ）に、もう一度LVTを実行し高次の調波成分を含めた全ての成分のパラメータを推定する。LVTの出力は、第一段階の場合と同様、Trajectory Estimationモジュールで時間方向に接続する。この時点で得られるのは、位相φに対する瞬時対数振幅Ａｋ（φ）と瞬時位相Ｐｋ（φ）である。最後に、これらの関数の位相軸φを時間tに逆変換し、各成分の瞬時対数振幅関数Ak(t)と位相関数Pk(t)を得る。以下に、各計算ブロックの詳細について述べる。

（３−２） Local Vector Transform
LVCでは入力信号x(t)に対して、一定の時間ステップΔごとに周波数応答を計算する。各分析時刻の近傍では、入力信号の各成分はQPSで近似できると仮定し、その時刻におけるQPSのパラメータを計算する。

まずｎ番目の時間フレームにおいて、入力信号に対して時刻ｔ＝ｎΔ中心とする窓関数w(t)を乗じる。次にx(t)の複素周波数応答X₀(ω)、およびその角周波数に対する一次、二次微分X₁(ω),X₂(ω)を次式数１６、数１７、数１８で計算する。

上式から明らかなように、X₁(ω), X₂(ω)の計算は、実際には信号x(t)に(-jωt),(-jωt)²を乗じた信号を入力とした短時間フーリエ変換（FFT）で実装可能である。よって、各時間フレームに対して３回のFFTを実行することになる。これらの周波数応答が得られれば、あとは数１２〜数１４のLVTで全ての角周波数に対するパラメータが一意に推定できる。

（３−３） Trajectory Estimation：最適パスの決定
LVTによって、各時間フレームの全ての周波数におけるQPSパラメータが計算できる（例えば、図２）。この推定されたパラメータは、その周波数近傍に単一のエネルギーしか存在しない場合は有効であるが、複数の成分が干渉している周波数領域では大きな誤差を含むものとなる。例えば図２の210, 350, 490, 630Hz周辺では、推定されたパラメータ値が大きく変動していることが確認できる。各成分の瞬時対数振幅関数Ak(t)と位相関数Pk(t)を得るためには、LVTで推定されたパラメータの中から、意味のあるものを抽出する必要がある。

この処理は、時間フレームｎにおいて、成分kとして妥当な推定パラメータを出力する角周波数ωBEST(k,n)を定める問題と考えることができる。この計算処理を最適パスの決定と呼ぶ。最適パスを決定するための拘束条件は、(1)LVT出力値の安定性、(2)出力されたパラメータの時間連続性の二つである。

図２に示したように、単一成分のエネルギーが支配的な周波数領域では、LVTの出力パラメータのばらつきは小さい。６つのパラメータa₀〜p₂は、周波数が140, 280, 420, 560, 700 Hz近傍で、ほぼ一定の値となっている。即ち、パラメータ推定値のばらつきが小さい周波数ほど、最適パスとして妥当ということになる。これを評価するために、角周波数ω近傍で推定されたパラメータの標準偏差σ(n,ω)を計算する。

また、時間フレームn、角周波数ωにおいて推定されたパラメータには、時間変化率が含まれる。有声音声の瞬時振幅や位相は連続的に変化するため、異なる時間フレームにおける推定パラメータの時間連続性を比較すれば、それがひとつの成分に起因するものか否かを評価することが可能である。ある点G1= (n1,ω1)において推定されたパラメータの値と、点G₂=(n₂,ω₂)のパラメータの間の時間連続性C₁₂は、次式数１９、数２０を用いて評価する。

ここでg(t)はtの一次関数であり、χ₁₂はG₁とG₂のパラメータ分布から計算される正規分布をg(t)が満たす結合確率に対応する。μdはパラメータ推定値である。この時間連続性の計算は、QPSの６つのパラメータのうち、a₀とp₀を除く４つから計算する。これは位相回転の問題を回避するためである。通常、sinusoid信号の瞬時位相は、時間tに対する単調増加関数になるが、各時間フレームで得られる瞬時位相は±πの範囲に限定されているからである。正しい瞬時位相を得るためには、この値に２πの整数倍を加えた値を用いる必要があるが、この整数値の推定と時間連続性を同時に評価するのは困難である。

最適パスは時間連続性C₁₂をコスト関数として動的計画法(dynamic programming)により計算する。時間−周波数平面における計算範囲を限定することで、複数の成分の最適パスを計算することも可能である。LVCの第一段階目の計算においては、50〜500 Hzの周波数範囲から第一調波成分の最適パスを算出し、第二段階では成分kに対してk+0.5〜k+1.5（単位は調波番号）の範囲からパスを計算する。

図４に最適パスの計算例を示す。図４は、男性話者の発話/iyoiyo/の、最初の/i/から/y/への音韻変化部において決定された最適パス（音声は図２と同一）を示している。決定されたパス上の(a₁,p₁)を瞬時値とし、その変化率(a₂,p₂)を用いて、各点のパラメータ値を三次元ベクトルで表示している（実線）。各ベクトルの方向が、近傍の時間のベクトルと滑らかに連続していることが確認できる。

（３−４） Trajectory Estimation：軌跡関数の生成
最適パスの計算により、N点（t₀〜tN_N-1 ）におけるQPSパラメータの値が得られている。これらの値に基づいて成分kの瞬時対数振幅Ak(t)と瞬時位相Pk(t)を決定する。まず前節で述べた位相回転の問題を解決する。推定すべき瞬時位相関数Pk(t)は時間に対する単調増加関数であるが、各時間フレームから得られるp0の値は±πに限定されているため、必ずしも単調増加にはならない。これは、時間フレーム間で２π×u（uは正の整数）の位相回転の情報が欠落したために生じる。正しい位相関数を得るためには、各時間フレーム間の回転係数ｕ_n 定める必要がある。

時間フレームnとn+1の間の位相回転係数ｕ_n は、次式数２１で求めることができる。ここでv(t)は瞬時角周波数を表す三次関数で、その係数はp₁(n),p₂(n),p₁(n+1),p₂(n+1)から一意に計算できる。Δは時間フレーム間隔である。

上式で位相回転の問題を解決した後、全時刻のパラメータ推定値から、振幅関数と位相関数を決定する。これらの関数は、ひとつの時間フレーム間では時間tに対する５次関数で表現される。推定されたパラメータ値からこの５次関数の係数は一意に決定可能だが、この手法で推定された関数はしばしば振動的になることがある。これはパラメータ推定時の微妙な誤差が、増幅されるために生じる。

この振動を抑制するために、各５次関数の係数には以下の拘束条件を用いる。(1)計算された各5次関数の接続点において、時間に対する４次微分までが連続であること。(2)計算された５次関数の接続点において、瞬時値、一次、二次微分が、推定パラメータを平均値、最適パス決定で用いた分散値で表される正規分布の結合確率を最大にすること。このような拘束条件を用いて各５次関数の係数を決定すると、上で述べた振動を抑制することが可能である。

図５に最適パスのパラメータ値と分散から計算した瞬時対数振幅関数、瞬時位相関数の例を示す。図５は、軌跡関数の生成例で、瞬時対数振幅関数、位相関数の計算例（男性話者/iyoiyo/）を示している。図５ (a)は、太線が計算された瞬時対数振幅関数を表す。白丸は最適パス上のパラメータ値、エラーバーはその分散に対応する。 (b)は、瞬時位相関数（位相回転の補正後）、 (c)は、瞬時対数振幅の変化速度関数、 (d)は、瞬時周波数関数、 (e)は、瞬時対数振幅の変化加速度関数、 (f)は、瞬時周波数の変化速度関数を示している。軌跡関数の計算においては、上述した時間的な滑らかさも拘束条件として用いているため、各軌跡は推定パラメータとは完全には一致しない。しかし、振動が少なくパラメータ値に近い軌跡が生成できていることが確認できる。

（３−５） Time to Phase Conversion
入力信号の成分間の干渉が小さい場合は、上記の手法により各成分の軌跡関数が簡単に求められる。しかし通常の有声音声では、特に高次の調波成分において成分間の干渉が無視できなくなるため、この様な手法が適用できない。図６にこのような例を示す。図６は、時間／位相変換の例を示しており、 (a)は、男性話者/iyoiyo/の/i/から/y/への遷移部分の波形とスペクトログラムで、エネルギーの強い部分を白で表示してある。F0は140 Hzから170 Hzへ上昇している（図５参照）。 (b)は、(a)の時刻450 msにおける周波数振幅応答を示している。F0の変化速度が大きいため、高次の調波成分（2.0 kHz以上）間で干渉が強くなり、振幅ピークを二次関数で近似できなくなっている。 (c)は、(a)の信号を時間／位相変換した後の波形とスペクトログラムで、スペクトログラムの横軸は、第一調波成分の位相、縦軸は調波番号を表す。F0がほぼ一定であることが確認できる。 (d)は、(c)の時刻に対応する位相における周波数振幅スペクトルで、成分間の干渉が軽減し、高次の調波成分の振幅ピークも明確に現れている。この図６(a) は、男性話者の発話/iyoiyo/の波形とスペクトログラムであり、この図の時間範囲は図５と同一である。時間400 msから500 msにかけてF0が急速に上昇している。図５(f) からF0の変化速度は最大0.8 Hz/msほどであることが分かる。第ｎ調波成分の周波数変化速度はこのｎ倍になるため、成分間の干渉が無視できなくなる。実際、図６(c) のt=450 msにおけるスペクトルでは、2 kHz以上の成分は明確なピークを確認できない。

この問題に対応するため、時間／位相軸変換を導入する。この変換は、干渉の影響の少ない第一調波成分の位相PF_F0(t) をLVTにより推定し、入力信号の時間軸tをこの位相軸に置き換えることで実現される。実際にはx(t)は離散時間でサンプリングされているため、任意の時刻tにおけるx(t)の値はsinc関数で補完する。変換後の信号x(φ) においては、時間軸の代わりに位相(rad)が用いられ、これを周波数分析すれば、角周波数の代わりに調波成分の番号が用いられることになる。これを次式数２２、数２３に示す。

この様な変換により、F0の時間変化に起因する成分間の干渉を低減することが可能である。図６(b) に時間／位相変換を施した波形とスペクトログラムを示す。このスペクトログラムでは、成分周波数（単位は調波番号になる）の時間変化がほぼ無くなっていることが確認できる。また図６(d) は、図６(c) と同じ時刻のスペクトルだが、高次の成分まで明確なピークが確認できる。

このようにして得られた瞬時振幅関数Ak( φ) や瞬時位相関数Pk( φ) の位相軸は、次式数２４、数２５のようにP _F0(t) を用いて容易に時間軸に逆変換することが可能である。

LVCの出力は、N点×K成分のQPSパラメータ（6×N×K）と、N点の基本位相パラメータ（3N）の和である、3N(2K+1)個の値で表現できる。

（４）合成音声を用いた性能評価実験
次に、LVCシステムの性能評価について説明する。まず定量的な評価を行うため、入力信号はパラメータが既知である合成音声を用いる。評価項目はF0検出誤差と、パラメータ推定性能とする。

（４−１）入力信号
評価用の入力信号として、合計400個の合成二連母音を作成する。刺激のパラメータは、母音の組み合わせ（25通り）、F0の変化パタン（８通り）、話者（男性、女性）の３種類である。入力信号のF0は次式数２６で表される。ここで、FCとFDはF0の平均的な大きさと変化幅に対応する定数であり、男性話者を模擬する信号の場合は FC = 100 Hz, FD = 10 Hzとし、女性話者ではFC = 200 Hz, FD = 20 Hzとする。定数Lは信号の持続時間で全ての信号で300msである。

持続時間中のF0の変化パタンは、パラメータFPにより決まる。FPは0から1.75まで0.25刻みで８通り用意する。例えばFp=0の場合は、F0は持続時間中に滑らかに減少し、FP=0.5では時刻150 msで頂点を持つ山型の変化パタンを示す。図７(a) に代表的なF0の変化パタンを示す。図７(b) は、フォルマント周波数で、男性話者の二連母音/ia/のF1-F4 を示しており、 (c)は、各調波成分の瞬時振幅で、２次〜15次の成分を表示している。

二連母音は、日本語５母音二つの組合せである。各母音に対応する周波数応答は５フォルマントのcascade-Klatt音声合成器で計算する（Klatt, 1980）。母音のフォルマント周波数は、男性話者を模擬するものと、女性話者を模擬するものの二通り用意する。表１に入力信号のフォルマント周波数を示す。尚、これらの値は、成人男女各一名の自然発話音声から抽出した値に基づいている。

通常Cascade-Klatt音声合成器はサンプリング周波数10kHzで実装されるが、これは声道長が約17cmの男性話者を近似している。一般に女性の声道長は男性よりも短いので、女性話者を模擬するためにはKlatt合成器のサンプリング周波数を上げる必要がある。ここでは女性話者のサンプリング周波数は12 kHzに設定する。これは声道長にして約14cmに対応する。

またKlatt合成器のナイキスト周波数以上の周波数応答H(f)は、次式数２７、数２８を用いて決定する。ここでK(f)はKlatt合成器の周波数応答、F5は第５フォルマント周波数を表す。H(f)は第５フォルマント周波数（男性：4.5 kHz、女性：5.4 kHz）以上の周波数で、振幅応答が-48 dB/octで減衰し、位相応答はゼロになることを意味する。この処理により、Klatt合成器のナイキスト周波数以上の領域で、連続的な応答を計算することが可能になる。図８に計算されたH(f)の例を示す。図８は、合成母音の振幅包絡の例を示しており、(a) は、母音/i/の振幅包絡で、線の太さが話者に対応し、太線が男性、細線が女性を表す。(b) は、母音/o/の振幅包絡を示している。

フォルマント周波数は、先行母音から後続母音へ滑らかに変化する。図７(b) に二連母音/ia/のフォルマント周波数を示す。

これらの情報に基づいて、入力信号は次式数２９、数３０により作成する。成分数Kは男性・女性どちらの話者を模擬する場合も80である。信号のサンプリング周波数は48kHz、持続時間は300 msとする。図７(c) に各調波成分の振幅パタンの例を示す。

（４−２） F0検出性能
音声のF0検出（pitch determination）に関しては、既に多くの研究がなされており、様々なアルゴリズムが提案されている。F0検出アルゴリズムの性能を評価する手法として、Rabiner(1977)は次の誤差率e(%)を提案している（数３１）。本明細書でも、この誤差率を用いて性能を評価する。

また、既存手法との性能差を評価するために、代表的なF0検出手法であるケプストラム法（Noll, 1966）と比較する。ケプストラムによるF0検出は、次のように実現する。まずLVCと同一の時間フレームに対して40 msのハミング窓をかけ、8192点のFFTを用いてパワースペクトルを計算する。次に各周波数成分のパワーを対数に変換したスペクトルを入力として、もう一度FFTを実行することでケプストラムを得る。このケプストラム上でエネルギーがピークを持つquefrencyが、信号の基本周期（1/F0）に対応する。

ここでは、ケプストラムのF0検出精度を向上させるために、予めF0の存在する範囲を限定している。Fc=100 Hzの入力に対しては、quefrencyのピークを求めるF0の探索範囲は80〜120 Hzとし、Fc=200 Hzの入力に対しては、160 Hzから240 Hzを探索範囲とした。尚、この条件はケプストラム法の検出性能を向上させるために用いており、LVCには適用していない。LVCは全ての入力信号に対して同一のパラメータでF0検出を行う。

図９に、LVCとケプストラム法のF0検出誤差を表す。この図９から入力信号の種類に関わらず、LVCのF0検出誤差は常にケプストラム法より小さいことが分かる。全入力に対する平均誤差は、LVCで0.000650%、ケプストラム法で0.163%である。この結果は、パラメータの時間変化を伴う有声音声のF0検出に対して、LVCシステムが有効であることを支持するものである。

図９はF0検出性能の評価を示しており、横軸は合成音声の母音の種類、縦軸はF0検出誤差（％)を示す。LVCのF0推定誤差を○と●で、cepstrum法（ケプストラム法）による検出誤差を◇と◆で表す。エラーバーはF0変化パタンによる標準偏差に対応する。

（４−３）パラメータ推定性能
ここではLVCのシステム全体としての性能を評価する。入力音のパラメータが既知であるので、各成分ごとにパラメータの推定精度を算出することも可能だが、入力パラメータが未知である場合と対応付ける為、推定されたパラメータから再合成された信号と入力信号の比較によってシステム全体の性能を評価する。具体的には入力信号x(t)と、推定されたパラメータから再合成した波形y(t)を用いて、次式数３２で表す入力/残差信号比（Signal to residual power ratio: S /R）を計算する。S/Rは入力信号と再合成信号の残差が小さいほど値が大きくなる。再合成信号y(t)は全ての成分の和から計算されるため、一部の成分のパラメータ推定における誤差が、そのままS/Rに反映される。高いS/Rを得るためには、全ての成分のパラメータ推定で誤差が小さくなっている必要がある。

図１０に各母音ごとのS/Rを示す。S/Rは母音の種類によって変化し、特に入力合成音声のフォルマント周波数が変化しない場合（５種類）と、変化する場合（20種類）の差が大きい。前者の平均S/RはFc=100 Hzで65.9 dB、Fc=200 Hzで69.1 dBと非常に高いが、後者の平均S/RはFc=100 Hzで38.4 dB、Fc=200 Hzで43.5dBである。また全サンプルの平均S/Rは46.3dBで、S/Rが30 dBを下回るサンプルは存在しない。

（５）自然音声を用いた性能評価実験
合成音声を用いた性能表実験により、LVCはF0や振幅包絡が時間変化する有声音声を高い精度で分析できることが示された。ここでは、自然発話音声を入力とした場合の性能を評価する。自然発話音声のF0は未知のため、評価項目はS/Rのみとする。

（５−１）入力信号
入力音声はATRディジタル音声データベースの音素バランス216単語より抽出する。データベースに含まれる単語の中から、全ての音韻が有声音であるサンプルを探し、表２で示す12単語を選定する。各単語音声の開始、終了時刻は、データベースに付属するラベル情報のうち、音響イベント層ラベルにより決定する。男性17名と女性17名（M101-M117, F101-F117）が静かな部屋で発話した、これら12単語を性能評価の入力信号として用いる。合計サンプル数は384である。

（５−２）パラメータ推定性能
合成音声の場合と同様に、自然発話音声の各成分の瞬時振幅関数と瞬時位相関数をLVCにより推定し、得られたパラメータから信号を再合成する。図１１(a)に入力信号と再合成、残差信号の例を示す。入力は男性話者が発話した/yumoa/であり、S/Rは27.6dBである。この図から再合成信号は、ほぼ入力信号と一致していることが確認できる。また図１１(c)に、これらの信号のスペクトログラムを示す。入力音声のF0と振幅包絡は、激しい時間変化を有する。また残差信号のスペクトル（図１１(c)下段）は、入力の振幅包絡のエネルギーが強い時間＝周波数領域でエネルギーが残っていることを示している。図１１は、自然音声のLVC分析例を示しており、(a)は、男性話者の発話/yumoa/の例で、上段が入力信号、中段がLVCの再合成信号、下段が残差信号にそれぞれ対応する。入力と残差信号の比(S/R)は27.6dBである。 (b)は、女性話者の発話/yumoa/の例で、図の形式は(a)と同じである。S/R=30.4 dB である。(c)は、(a)の信号のスペクトログラムで、上段が入力信号、中段が再合成信号、下段が残差信号を表す。 (d)は、(b)のスペクトログラムで、図の形式は(c)と同一である。

図１１(b)と(d)には、女性話者の/yumoa/の入力、再合成、残差信号を同じ形式で示す。これも男性話者の場合とほぼ同様の結果が得られていることが確認できる。特に、残差信号のスペクトルの時刻0-200ms、周波数500Hz周辺のエネルギーは重要である。この残差エネルギーは、入力信号の第一調波成分と第二調波成分の間に存在している。即ち、有声音信号が振幅と周波数が連続的に変化するsinusoid信号の和として表現できるするLVCの仮定の枠外にあるエネルギーということになる。

Fantの音源フィルタ理論によれば、有声音声の発話時には声帯振動が主な音源ではあるが、同時に乱流雑音もわずかながら発生していると考えられる。この雑音音源も、声帯振動音源と同じ共振フィルタと同じ声道フィルタで変調されるため、入力信号の振幅包絡の強い領域で残差エネルギーが強くなるのである。声帯振動に起因する音源信号は周波数とともに減衰することが知られているが、乱流雑音音源にはそのような特性はない。よって、残差信号の大きさは高周波数になるほど大きくなることが予測できる。図１２に、その他の発話の入力信号と残差信号の例を示す。特に図１２(c)の結果は、上で述べた仮説を支持するものである。図１２は、入力信号スペクトル（上段）と、残差信号スペクトル（下段）の比較を示したものであり、 (a)は、男性話者の/nyuIN/の例で、S/R=31.2 dB、 (b)は、女性話者の/nyuiN/の例で、S/R=33.2 dB、 (c)は、男性話者の/reNai/の例で、S/R=27.6dB、 (d)は、女性話者の/reNai/の例で、S/R=30.5 dB の場合を、それぞれ示している。

図１３に各単語の平均S/Rを示す。男性話者の発話単語に対する平均S/Rは21.2ｄB、女性話者の平均S/Rは24.0dB、全サンプルの平均S/Rは22.6dBである。

（６） LVCに基づく音声変調
以上の実験により、LVCが自然発話された有声音声の成分パラメータを高い精度で推定できることを示した。このようなパラメータで有声音声を表現することにより、音声処理の観点から様々な応用が期待できる。ここではその一例として、音源フィルタ理論を応用した音声の変調手法について説明する。

（６−１）発話速度の変調
LVCの出力は、音声の各成分kに対する瞬時対数振幅関数Ak(t)と位相関数Pk(t)である。Fantの音源フィルタ理論を適用するために、声道フィルタに起因する振幅包絡E(t,ω)を次式数３３で計算する。ここでInterpolateは、ある時刻の全ての成分の瞬時角周波数と瞬時振幅値を用いて、任意の角周波数における瞬時振幅値を得る関数とする。これは、例えばωに最も近い二つの成分の瞬時振幅の線形補完などで実現できる。

この表現に基づいて、音声の発話速度を任意の速さに変調する手法を考える。単純なサンプリング周波数を変えて発話速度を変更する手法では、元の音声とは異なる話者のような変調音声が生成される。これを回避するためには、包絡Eが角周波数軸に対して変化しないことと、F0の範囲が元の音声と同一であることが必要となる。具体的には、発話速度をＭ_T 倍した変調音声の瞬時振幅A’(t)と位相P’(t)は次式数３４、数３５で与えられる。

上記の変換により、任意の発話速度で同一話者と感じられる滑らかな音声される。

（６−２）振幅包絡とF0の変調
音声に含まれる話者の個人性は声道フィルタに関係する振幅包落と、声帯振動と関係するF0に強く影響される。例えば声道長の短い女性の音声を変調して、声道長の長い男性の音声を生成するためには、包絡関数の周波数軸を変化させれば良い。また声の高さを変化させるなら、瞬時位相関数を変更すれば良い。具体的には、話者の声道長をMV倍、基本周波数をMP倍に変調した音声の瞬時振幅A’(t)と位相P’(t)は次式数３６、数３７のように書ける。

上記変換により、男性話者の発話から女性話者と感じられる変調音声を生成することができる。無論、その逆も可能である。

（６−３）音声モーフィング
上述した振幅包絡とF0変調手法を応用すると、ある音声から別の音声への連続的なモーフィングが可能となる。この手法の特徴は、モーフィングの途中の音声でも自然な発話のように聞こえる点である。具体的には、二つの音声の包絡関数Ｅ_A ，Ｅ_Bと位相関数Ｐ_A ，Ｐ_B、及び話者の声道長の比R₁₂と二つの音声の混合比Ｍ_X が与えられている場合に、モーフィング音声の振幅A’(t)と位相P’(t)は次式数３８、数３９で計算できる。

図１４にモーフィングの例を示す。図１４は、発話/warauoNna/ に対するモーフィングの例を示しており、 (a)は、女性話者のスペクトログラム (b)は、(a)と(f)を混合比20%でモーフィングした変調音声のスペクトログラム、 (c)は、混合比40% 、(d)は、混合比60% 、(e)は、混合比80% 、(f)は、男性話者のスペクトログラム、をそれぞれ示している。図１４(a)が女性の発話、図１４(f)が男性の発話であり、混合比を変えていくとF0と振幅包絡が滑らかに変化することが確認できる。

このように、本発明においては、既存の技術では困難であった時間変化を伴う有声音声の正確な分析が可能となった。これにより、本発明に係る技術を用いて、多様な応用が期待できる。本願では応用のひとつとして音声変調を例に挙げた。また、自然音声の時間変化を詳細に分析する手法は、今後の音声生成研究（音声合成研究）の発展に大いに寄与できると考えられる。

例えば、本発明による音声信号合成までの概略ステップを表したシステムの構成例をず１５に示す。このように、本発明により、位相スペクトルの時間軸変換まで実行することにより、極めて高精度の音声信号分析を行うことが可能になり、かつ、分析された信号を用いて、実質的に自由に精度良く音声信号を再現、合成できるようになる。

なお、本願では、音声を主体に説明してきたが、楽器音についても、同様の適用、展開が可能であることは言うまでもない。

QPSの周波数応答を示す特性図である。 LVTの計算例を示す特性図である。 LVCシステム構成例を示すブロック図である。 最適パス決定の例を示す特性図である。 軌跡関数の生成例を示す特性図である。 時間／位相変換の例を示す特性図である。 性能評価の入力信号の例を示す特性図である。 合成母音の振幅包絡の例を示す特性図である。 F0検出性能の例を示す特性図である。 パラメータ推定精度の例を示す特性図である。 自然音声のLVC分析例１の例を示す特性図である。 自然音声のLVC分析例２の例を示す特性図である。 自然音声のパラメータ推定性能の例を示す特性図である。 音声モーフィングの例を示す特性図である。 本発明に係るシステムの構成例を示す説明図である。

Claims (5)

1. 時間変化特性を持つ入力音響信号に対して、入力信号の正弦波成分の瞬時振幅、瞬時周波数、瞬時位相を下記の４段階の計算処理に基づいて、再合成可能なパラメータ関数として出力することを特徴とする音響信号分析方法。
   （１）入力信号〔数１〕の任意の分析時刻において、正弦波パラメータを〔数２、３〕のように近似し、パワースペクトルと位相スペクトルから、〔数４〜６〕のようにこれらのパラメータ候補を定める局所変化率変換処理を行う段階
   （２）各分析時刻で得られた瞬時振幅と位相のパラメータ候補を、動的計画法を用いた時間連続性によりグルーピング化し、これらのパラメータを連続関数として出力するパラメータ軌跡推定処理を行う段階
   （３）入力信号〔数１〕に含まれる成分間の干渉により、上記（１）と（２）の段階では十分なパラメータ推定精度が得られない場合に、入力信号の基本位相のみを推定した後で、入力信号の時間軸を基本位相軸に変換した信号を生成することで成分間の干渉を抑制する時間−位相変換処理を行う段階
   （４）上記（３）の段階で得られた変換信号に、（１）と（２）の段階を適用して正弦波パラメータを推定し、得られた連続関数の基本位相軸を時間軸に逆変換することで入力信号の瞬時振幅と位相を連続関数として出力する時間−位相逆変換処理を行う段階
   

   

   

   

   

   

   

   ここで、x(t)は有声音信号、Kはsinusoidの数を表し、ξk(t),ψk(t), ηkはｋ番目の正弦波（sinusoid）の瞬時振幅、瞬時周波数、初期位相にそれぞれ対応する。
   τは変数、t,tcは時刻、A(t)とP(t)は時刻tc近傍での瞬時対数振幅と瞬時位相、a ₀ ,a ₁ ,a ₂ ,p ₀ ,p ₁ ,p ₂ は定数であり、a ₀ が瞬時対数振幅、a ₁ が振幅の変化速度、a ₂ が振幅の変化加速度、p ₀ が瞬時位相、p ₁ が角周波数、p ₂ が角周波数の変化速度にそれぞれ対応する。
   S ₀ , S ₁ , S ₂ は任意の角周波数ωにおける３種類の周波数応答、α ₀ ,α ₁ ,α ₂ は任意の角周波数ωにおける対応する３種類の信号パラメータである。
2. 請求項１に記載の音響信号分析方法における各パラメータ関数を用いて音響信号を再合成することを特徴とする音響信号合成方法。
3. 再合成された音響信号から、または、音響信号を再合成するに際し、音質を他の音質に変換する、請求項２に記載の音響信号合成方法。
4. 音響信号が音声信号からなる、請求項１〜３の方法。
5. 音響信号が楽器音信号からなる、請求項１〜３の方法。

Images