このファイルは、３年春学期科目「音声情報処理」の授業において、本来であれば口頭で説明する内容をテキスト(.txt)に書き起こしたものです。

事前に科目のＨＰにアクセスして「講義ノート」の PDF ファイルをダウンロードし、印刷しておいてください。このテキストには講義ノートの記述より多くの情報を含んでいます。講義ノートの行間に、テキストの説明を書きこんで行くことをお勧めします。最後にプリントを綴じれば、将来必要になったときに見返せる自分の「講義ノート」ができあがります。

「学生用ページ」「レポート提出システム」 から、その日の「宿題」レポートを提出してください。 〆切は 次回授業日の前日 23:59 です。

_/_/_/ 第六回　線形予測分析(2) _/_/_/

【講義ノート】音声生成の線形モデル

このページの図は、ケプストラムの回にお見せした「発声機構と信号表現」の図とよく似ています。下の方に音源があって共鳴、放射とつながっています。

今回の図が前とちょっと違うのは、音源の部分を「インパルス」と「生成」に分けて書いているところです。

４つの緑の枠内に下から「インパルス」「生成」「共鳴」「放射」と記入しましょう。

声帯の振動のメカニズムを、声帯という伝達系をインパルスが駆動していると考えることによって、その特性を「e(t)」と「g(t)」に分けます。そして「e(t)」のみを信号源とみて、続く「g(t)」「v(t)」「l(t)」を合わせて伝達系とみます。

e(t) の右の方に「信号源」、l(t),v(t),g(t) の右の方に「伝達系」と記入。

このような分け方をすることで、後に述べる「線形予測」の考え方を導入できるようになります。

さて、この「e(t)」「g(t)」「v(t)」「l(t)」「s(t)」は時間の関数ですが、これらの z-変換であるスペクトル「E(z)」「G(z)」「V(z)」「L(z)」「S(z)」に対しては、次の関係があります。

S(z) = E(z)G(z)V(z)L(z)

フーリエ変換スペクトルのときと同じです。z-変換スペクトルにおいても、総合的スペクトルは、各要素のスペクトルの「積」になります。


【講義ノート】音声生成の線形モデル (2)

各構成要素の z-変換スペクトルを見て行きましょう。

まず、声帯音源を分割した「インパルス列」のスペクトルは、

          ∞                   σ
E(z) = σ Σ (z^{-I})^n = -----------
          n=0              1 - z^{-I}

この式のようになります。

前に勉強したように z^{-1} は 1 標本点の時間遅れに対応しますので、要素「z^{-I}」は、ある瞬間に入力されたインパルスを「I」標本点だけ遅れて出力します。この I 点遅れを「Σ」で繰り返すことによって、I 標本点を周期とするインパルス応答が表現できます。

「(z^{-I})^n)」は、公比「z^{-I}」の等比数列です。よってその和は、等比級数の公式によって、「σ/(1-z^{-I})」と求まります。「σ」は定数です。

音声が周期性を持つとき、その基本周期を「ピッチ」周期と言います。標本化(サンプリング)周期を T とするとき、ピッチ周期 P は I × T で計算できます。

T: サンプリング周期
P: ピッチ周期


【講義ノート】音声生成の線形モデル (3)

続いて、声帯音源を分割した「生成」のスペクトルです。
詳しい説明は省きますが、いろいろな研究の結果、

                 1
G(z) = -----------------------
        (1 - e^{-cT}z^{-1})^2

という式になることが分かっています。

また、口唇における「放射」のモデルは次式のようになります。

L(z) = 1-z^{-1}

これらは、こういうものだ、くらいに思って頂ければよいです。


【講義ノート】音声生成の線形モデル (4)

さて、いよいよ音声認識に直結する、声道形状による「共鳴」のモデルです。

V(z) = ...

ややこしい式ですが、この式が一番重要な式です。式の中に、b_i という係数が出てきますね。これが「ホルマント」周波数と関係します。もう一つ、 c_i という係数が出てきます。これが「帯域幅」という概念と関係し、共鳴の強さを示します。

式の分母に「Π」(ギリシャ大文字「パイ」)という記号があります。これは続く [ ] の値を、 i=1 から k まで全部掛けることを意味する記号です。i がホルマントの番号、k がホルマントの個数に対応します。

この式は複雑すぎて覚えるのは難しいですが、どのようにしてこの式が導かれるかは、知っておく必要があります。次ページ以降で説明していきます。


【講義ノート】音声生成の線形モデル (5)

この図は、声道形状による「共鳴」を z-平面 で理解するための図です。

z-平面とは、複素数 z の実部を横軸に、虚部を縦軸に配した平面です。

z-平面中の単位円の円周上を、複素数 z が移動することによって、スペクトルの横軸である「周波数」が変化します。原点(0, 0)と円周上の複素数 z を結ぶ線が、実軸となす角度が周波数に対応します。

z=(1, 0) のとき、角度は 0 になり、周波数 0Hz(直流)に対応します。
z が、単位円周上を左回りに移動することによって、対応周波数が増加していき、
半周して z=(-1, 0) のとき、標本化周波数の半分の周波数に対応します。

さて、この z-平面の中で、一つのホルマントを実軸に対称な２つの「共役」複素数で表します。図中の赤い点です。式で書くと

z  = exp( -cT + j bT )
z* = exp( -cT - j bT )

これらをこの共鳴特性の「極」と言います。

一つのホルマントについて２つずつの複素数が対応します。例えば、第三ホルマントまで考えるのであれば、赤い点は６個になるわけです。


【講義ノート】極とスペクトル

共鳴によって生じる、音声スペクトル上のエネルギー集中のことを「ホルマント」というのでした。

z-平面中の単位円の円周上を、複素数 z が (1, 0) から (-1, 0) まで、左回りに半周するとき、z の値の変化によって伝達関数 X(z) の値も変化していきます。

z に対応する周波数を横軸にとって、X(z) の値を縦軸にグラフを描くとグラフが描けますね。そうです、それが z-変換 のスペクトルです。

次式のような、伝達系 X(z) を考えます。

                               1
X(z) = -------------------------------------------------
         k
        Π　(1- (z_i)×(z^{-1}) (1- (z_i)* × (z^{-1}))
         i

複素数 z を (1, 0) から (-1, 0) まで、円周に沿って動かしていくとき、z が複素数 z_i (赤い点)の近くを通過すると、(1- (z_i)×(z^{-1}) の値が 0 に近づき、分母が 0 に近づくので、X(z) の値は無限大に近づきます。言い換えると、スペクトル X(z) は z_i の近くで「ピーク」を持つことになります。

伝達系にこのような式を用いることによって、設定した特定の複素数の近くでスペクトルにピークを作ることができることが分かりました。スペクトルのピークとは、その周波数におけるエネルギー集中です。「スペクトル上のエネルギー集中」すなわち「ホルマント」ですね。

伝達系を z^{-1} の有理多項式、すなわち (z^{-1} の多項式)/(z^{-1} の多項式) の形で表現するとき、分母多項式=0 とおいた方程式の「根」を「極」(pole) と呼びます。

以上をまとめると、結局『ホルマントは伝達系の「極」として表現される』ことが分かりました。


【講義ノート】音声生成の線形モデル (6)

以上述べてきた、生成 G(z), 共鳴 V(z), 放射 L(z) を掛け合わせて、伝達系の総合特性を求めると、この式のようになります。

G(z)V(z)L(z) = ...

研究者の知見により、cT << 1 であることが分かっていますので、分子の「(1-z^{-1})」と分母の「1-e^{-cT} z^{-1})^2」のうち 1個が約分できます。結局、総合特性は、この下の式のように求められます。

                                       1
 =  --------------------------------------------------------------------
                         k
    (1-e^{-cT} z^{-1})  Π　(1- (z_i)×(z^{-1}) (1- (z_i)* × (z^{-1}))
                         i

この式で重要なのは、分子の多項式が「1」になったということです。伝達特性が、分子が 1 である有理多項式で表現されるとき、その伝達系は「全極型」のフィルタになると言います。


【講義ノート】音声生成の線形モデル (7)

この分母の z の多項式を A(z) と書くとき、全極型音声生成のモデルは、
              1
S(z) = E(z) ------
             A(z)
となります。音源の駆動インパルス列 E(z) に、音声生成フィルタ 1/A(z) を掛けると、音声のスペクトル S(z) が生成できます。

この A(z) を「逆フィルタ」と呼びます。
             1
A(z) = --------------
        G(z)V(z)L(z)

何故逆フィルタというかというと、式を変形して、

E(z) = S(z)A(z)

と書くと、音声スペクトル S(z) に、フィルタ A(z) を掛けると、「逆に」音源の駆動インパルス列 E(z) に戻るという式が書けるからです。


【講義ノート】全極型音声生成フィルタ

「z^{-1}」は 1 標本点の時間遅れに相当します。よって、入力波形 x_n 、出力波形 y_n に対し次の関係がなりたちます。

x_n = y_n × ( z^{-1} の多項式 )

この信号の流れを図で表すとこのようになります。図中の「z^{-1}」は 1 標本点の遅れ、「+」は加算を意味しています。「a_i」は多項式の係数です。

コンピュータで時間遅れを実現するのは簡単です。ある時刻の値をメモリしておいて、次の(1 標本周期後の)時刻の計算に使えばよい。すると、図のような信号処理は、記憶と加算だけで実行できることになります。

このような演算によって、コンピュータによって伝達系のシミュレーションを行うことができます。伝達系の z-変換 さえ知っていれば、その入出力関係をコンピュータで計算することができるわけです。このような仕掛けを「ディジタルフィルタ」と呼びます。

この図は「全極型」ディジタルフィルタの信号の流れを示したものです。

ある時刻に入力された信号は、そのまま出力される部分もありますが、遅延を与えられ、ある係数を掛けられて、入力にフィードバックします。この結果、ループができて、ある一瞬に信号が入力すればそれがグルグル廻って、減衰しながらも無限に出力が続きます。

一瞬の「インパルス」信号の入力に対して、どのような信号系列が出力されるかを、その伝達系の「インパルス応答」といいます。

全極型のディジタルフィルタは、無限にインパルス応答が続く特性を持っています。
このような特性を IIR (Infinite Impulse Response) と言います。

全極型音声生成フィルタ 1/A(z) は、IIR フィルタの一種です。


【講義ノート】逆フィルタ

同じように、逆フィルタの信号の流れを見てみたいとおもいます。逆フィルタはこのような z-変換 で表されるのでした。

                           k
A(z) = (1-e^{-cT} z^{-1})  Π　(1- (z_i)×(z^{-1}) (1- (z_i)* × (z^{-1}))
                           i

これは z^{-1} の多項式です。式を展開して次のように書くことができます。

       M
A(z) = Σ　a_i × z^{-i}
       i=0

ここで M は多項式の次数で、ホルマントの数を k として、高々 2k + 1 です。

M が有限の値であることに注意してください。このフィルタにインパルス信号が入力されると、M 標本点より後に応答が出力されることはありません。

すなわち、逆フィルタは、有限のインパルス応答を持ちます。
このような特性を FIR (Finite Impulse Response) と言います。

逆フィルタ A(z) は、FIR フィルタの一種です。

ピンクの枠に「インパルスレスポンス長 M の FIR フィルタ」
とメモしておきましょう。


【講義ノート】逆フィルタ (2)

逆フィルタでは、入力波形 x_n 、出力波形 y_n に対し次の関係がなりたちます。

y_n = x_n × ( z^{-1} の多項式 )

この信号の流れを図で示します。 IIR の場合と違ってループができていないことが分かります。

逆フィルタ A(z) は、インパルス応答長 M の FIR フィルタになります。


【講義ノート】線形予測分析

ようやく準備が整いました、いよいよ線形予測によるスペクトル分析の話を始めましょう。

線形予測分析は、まずは音声の符号化に多用されています。
LPC: Linear Predictive Coding
という言葉を覚えておきましょう。

線形予測分析は、スペクトル解析の一手法です。音声の高圧縮符号化、認識、合成に使われます。また、地震波による地層解析、音響的故障診断など非常に応用範囲が広いです。

・スペクトル解析の一手法
　- 音声の高圧縮符号化、認識、合成
　- 地震波による地層解析、音響的故障診断

基本的な考え方としては、伝達系を全極型 IIR フィルタで表現し、音源としてインパルス列を仮定し、観測波形から逆フィルタの係数を求めるという考え方です。

・伝達系を全極型 IIR フィルタで表現
　- 音源としてインパルス列を仮定し、
　　観測波形から逆フィルタの係数を求める


【講義ノート】線形予測分析 (2)

今、言った考え方を、音声生成のモデルにあてはめると、このようになります。

声帯音源を、インパルス列と生成に分け、「生成」「共鳴」「放射」を伝達系として一つにまとめます。

e(t) の右の方に「インパルス列を仮定」と記入。それから、
s(t) の右の方に「観測された音声波形」と記入。

このように考えれば、問題を『「インパルス列」を観測された「波形」に変換する「伝達系」の逆フィルタ A(z) を求める問題』と整理することができます。

l(t),v(t),g(t) の右の方に『インパルス列を、観測された波形に変換する「伝達系」の逆フィルタ A(z)  を求める』と記入。ちゃんと意味を考えながら書くこと。


【講義ノート】線形予測分析 (3)

スペクトル領域における逆フィルタの式を再掲します。

E(z) = S(z)A(z)

音声スペクトル S(z) にフィルタ A(z) を掛けると、音源の駆動インパルス列 E(z) に戻るという式です。

スペクトル領域の掛け算は、時間領域では畳み込み和になります。すなわち、

        M
e(n) = Σ　a_i × s(n-i)
       i=0

ここで a_i は IIR フィルタの係数です。「線形予測係数」と呼びます。

i=0 のときだけ Σ の外に出します。a_0 = 1 なので、

              M
e(n) = s(n) + Σ　a_i × s(n-i)
              i=1
ここで、過去 M-1 点の音声波形より、現在の波形 s~(n) を予測します。
           M
s~(n) =  - Σ　a_i × s(n-i)
           i=1
過去 M-1 点の値に、重み係数 a_i を掛けて加算すること、すなわち線形演算によって、現在の値 s~(n) を予測します。この計算を「線形予測」といいます。

結局、

e(n) = s(n) - s~(n)

と整理されます。予測された音声波形の値と観測された音声波形の値の差が e(t) になるということです。


【講義ノート】線形予測分析 (4)

e(n) がインパルス列であることを思い出しましょう。インパルス列はインパルスのある時刻以外は値が 0、すなわちほとんどの n で e(n)=0 となります。

よって、インパルスを避けるように n_0, n_1 を設定し、二乗誤差「α」を最小化するように、線形予測係数 a_i を求めればよいのです。
      n_1
α =  Σ　　e^2(n)
      n=n_0
講義ノートのように式を変形し、結局、
       M   M
α = 　Σ　Σ　a_i × c_{i,j} × a_j
       i=0 j=0
という式が得られます。ここで重要なのが、新しく導入された
          n_1
c_{i,j} = Σ　　s(n-i)×s(n-j)
          n=n_0
という量です。

この量を音声波形 s(n) に対する「共分散」係数と言います。式から明らかに、この量は音声波形 s(n) だけあれば、計算することができます。


【講義ノート】線形予測分析 (5)

さて、二乗誤差αを最小化する、線形予測係数 a_i を求めましょう。

       M   M
α = 　Σ　Σ　a_i × c_{i,j} × a_j
       i=0 j=0

を、各 a_k (k=1, ... ,M) で偏微分して 0 とおきます。

        M
0  =  2 Σ　a_i × c_{i,k}
        i=0

a_0 = 1 は定数なので、未知数は a_k (k=1, ... ,M) の M個、
偏微分するのも a_k (k=1, ... ,M) の M個なので、結局、
M元の連立一次方程式が M個得られることになります。

M
Σ　a_i × c_{i,k} = - c_{0,k}　(k=1, ... ,M)
i=0

後はこの M個の式からなる M元連立一次方程式を解けばよいのです。

連立方程式の解として「線形予測係数」が求まります。

このように、線形予測に基づくスペクトル分析とは、音声波形から共分散係数を求め、共分散係数を含む連立方程式を解くことによって、波形予測の二乗誤差が最小になるような、線形予測係数を求めることです。

この線形予測係数を使って音声生成ディジタルフィルタを作り、そのスペクトルを計算します。それが、入力された音声波形に対するスペクトル分析結果です。


【講義ノート】線形予測分析 (6)

線形予測分析には、共分散係数の取り扱い方によって大きく分けて２種類があります。

一つは「共分散法」と呼ばれる方法で、M元連立一次方程式M個をそのまま代数的に解きます。この方法は、周波数の分解精度が高いのですが、若干安定性に欠け、推定された線形予測係数が変動しやすい欠点があります。

共分散法は、地震波による地層解析や、音響的故障診断に用いられています。

もう一つは「自己相関法」と呼ばれる方法です。共分散係数の計算をするときに、時間ずれ i, j を２自由度で扱わずに、その差「l」のみで平均し、近似する方法です。この方法は、周波数の分解精度が若干落ちるのですが、安定性が高いという大きな利点があります。

音声信号処理に用いられている線形予測分析は、だいたいこの自己相関法です。


【講義ノート】線形予測によるスペクトル概形分析

線形予測によるスペクトル分析が、ホルマント判定にどのくらい有効か、実例を観察してみましょう。

これは、母音 /o/ のスペクトル概形の例です。。

図の横軸が周波数。0 Hz 〜 6000 Hz の範囲が表示されています。縦軸は各周波数成分の強さ。

青い線がフーリエ変換を用いて求めたスペクトルの形状です。

緑の線がケプストラムとリフタを用いて、スペクトルの形状です。青い線のフーリエ変換スペクトルと比べて、緑の線は形状が滑らかで、ホルマントを判断しやすい。これは前に見た通りです。

さて、赤い線が線形予測によるスペクトル概形です。緑の線よりさらに形状を把握しやすく、ホルマントの位置が明確に抽出されていることがわかります。

線形予測によるスペクトル分析の威力を実感して頂けましたでしょうか？

図の右側に２つ黄色い枠があります。各々、
音源の周期成分を除いた「スペクトル概形」を求めることができる
ホルマント周波数を判定しやすい
と記入しておいてください。

--