線形予測法

線形予測法（せんけいよそくほう、英: linear prediction）は、離散信号の将来の値をそれまでの標本群の線型写像として予測する数学的操作である。

デジタル信号処理では、線形予測法を線形予測符号 (LPC) と呼び、デジタルフィルタのサブセットと見ることができる。（数学の一分野としての）システム分析では、線形予測法は数学的モデルや最適化の一種と見ることができる。

モデル

系列 $x(n)$ に対して $p$ 次の線形予測法で推定した値 ${\widehat {x}}(n)$ は予測係数 $a_{i}$ を用いて次で表される。

{\widehat {x}}(n)=-(a_{1}x(n-1) a_{2}x(n-2) ... a_{p}x(n-p))=-\sum _{i=1}^{p}a_{i}x(n-i)=-{\boldsymbol {a_{1:p}}}^{\top }{\boldsymbol {x_{n-1:n-p}}}

すなわち線形予測法とは $p$ 次の過去系列を用いた線形回帰である。

この誤差は多次元信号においてベクトルノルム $\|.\|$ を用いて次のように定義される。

e(n)=\|x(n)-{\widehat {x}}(n)\|\,

パラメータ推定

線形予測法における予測係数 $a_{i}$ には様々な推定方法が存在する。

最適化においてパラメータ $a_{i}$ の典型的な選択法は、二乗平均平方根基準であり、これを自己相関基準とも呼ぶ。これは、以下の式で得られる二乗誤差 E[e²(n)] の期待値を最小化する手法である。

\sum _{i=1}^{p}a_{i}R(i-j)=-R(j)

ここで 1 ≤ j ≤ p であり、R は信号 x_n の自己相関であり、次のように定義される。

\ R(i)=E\{x(n)x(n-i)\}\,

ここで E は期待値である。多次元の場合、これはL₂ノルムを最小化することに対応する。

上の式を正規方程式または Yule-Walker 方程式と呼ぶ。行列形式でこの方程式を表すと、次のようになる。

Ra=-r,\,

ここで、自己相関行列 R は対称なテプリッツ行列であり、その要素は r_i,j = R(i − j) である。また、ベクトル r は自己相関ベクトル r_j = R(j) であり、ベクトル a は係数ベクトルである。

より汎用的な形式として、次を最小化する方式もある。

e(n)=x(n)-{\widehat {x}}(n)=x(n) \sum _{i=1}^{p}a_{i}x(n-i)=\sum _{i=0}^{p}a_{i}x(n-i)

ここで、係数 $a_{i}$ について $a_{0}=1$ とし、自明な解を防ぐのが一般的である。これにより上述と同じになるが、正規方程式は以下のようになる。

\ Ra=[1,0,...,0]^{\mathrm {T} }

ここで、インテックス i の範囲は 0 から p、R は (p 1) 行 (p 1) 列の行列である。

パラメータの最適化は大きな問題であり、他にも様々な手法が提案されている。

その中でも自己相関手法が最もよく使われており、例えばGSMでの音声符号化に使われている。

行列方程式 Ra = r の解の計算は、比較的時間のかかる処理である。ガウスの消去法を使った解法が最も古くからあるが、R と r の対称性をうまく利用していない。より高速なアルゴリズムとして、1947年に Norman Levinson が考案したレビンソン再帰という再帰的解法がある。その後、Philippe Delsarte らが、これを改良した分割レビンソン再帰というアルゴリズムを発表した。これは、乗除算回数を約半分にしたもので、パラメータベクトルの特殊な対称性をそれぞれの再帰で利用する。

二乗予測誤差

元信号を用いた推定法の1つが予測誤差の二乗を最小化する手法である^[1]。線形予測値 ${\widehat {x_{n}}}$ と真の信号 $x_{n}$ の間の予測誤差 $e_{n}$ は $0$ 次の係数 $a_{0}=1$ を用いて次で表される。

e_{n}=1x_{n}-{\widehat {x_{n}}}=a_{0}x_{n} {\boldsymbol {x_{n-1:n-P}}}{}^{\top }{\boldsymbol {a_{1:P}}}={\boldsymbol {x_{n:n-P}}}{}^{\top }{\boldsymbol {a_{0:P}}}

ここで全時点にわたる予測誤差の二乗和をもって誤差関数 $e=\sum _{n}^{N}{e_{n}}^{2}$ とし、これを最小化して予測係数を推定する。最小値において予測係数の偏微分は $0$ になるので $a_{q}\ (1\leq q\leq P)$ において次が成立する。 ${\partial e \over \partial a_{q}}=\sum _{n}^{N}{\partial \over \partial a_{q}}({\boldsymbol {x_{n:n-P}}}{}^{\top }{\boldsymbol {a_{0:P}}})^{2}=\sum _{n}^{N}2x_{n-q}({\boldsymbol {x_{n:n-P}}}{}^{\top }{\boldsymbol {a_{0:P}}})=2\sum _{n}^{N}\sum _{p=0}^{P}x_{n-q}x_{n-p}a_{p}=2\sum _{p=0}^{P}a_{p}\sum _{n}^{N}x_{n-q}x_{n-p}=0$

ここで $\sum _{n}^{N}{x_{n-q}x_{n-p}}$ がラグ $\left\vert p-q\right\vert$ の自己相関関数になっていることから $r_{\left\vert p-q\right\vert }$ と表記すると $\sum _{p=0}^{P}r_{\left\vert p-q\right\vert }a_{p}=0$ となり、ベクトル表記すると ${\boldsymbol {r_{0-q:P-q}}}{}^{\top }{\boldsymbol {a_{0:P}}}=0$ となる。 $1\leq q\leq P$ 全てを行列に集約すると以下になる。

${\begin{pmatrix}r_{1}&r_{0}&r_{1}&...&r_{P-1}\\r_{2}&r_{1}&r_{0}&...&r_{P-2}\\...&...&...&...&...\\r_{P}&r_{P-1}&r_{P-2}&...&r_{0}\\\end{pmatrix}}{\begin{pmatrix}1\\a_{1}\\a_{2}\\...\\a_{P}\\\end{pmatrix}}={\boldsymbol {0}}$

すなわちユールウォーカー方程式に帰着する。この $p$ 元連立一次方程式を解くことで予測係数が求まる。解くにあたってガウスの消去法などの汎用解法が利用できるが、この行列はテプリッツ行列を導くからLevinson再帰を用いた高速解法が存在する。実務的にはこれがLPCでよく利用される。

スペクトル変換

信号の符号化（エンコード）では元信号を用いて予測係数を得ることが有用であるが、信号生成タスクでは元信号を利用することができない。そのような場合に利用可能な、元信号を直接用いない推定法の1つがスペクトルの変換である。この手法ではスペクトルを自己相関関数へ変換しユールウォーカー方程式へ持ち込みこれを解くことで予測係数を推定する。

この手法の中心にある原理はウィーナー＝ヒンチンの定理である。この定理はパワースペクトル密度のフーリエ逆変換が自己相関関数となることを示している。

ゆえにメルスペクトログラムやMFCCをパワースペクトル密度へ変換し、それにフーリエ逆変換を適用、得られた自己相関でテプリッツ行列を構成しユールウォーカー方程式を解くことで予測係数が得られる^[2]。

係数表現

線形予測における係数は複数の形式で表現できる^[3]。以下は代表的な係数表現である。

線形予測係数（英: linear predictive coefficients; LP coefficients）: 定義式における $a_{i}$ 。ノイズに対する脆弱性が知られる
line spectral frequencies; LSF
reflection coefficients; RC
autocorrelations; AC
ログ面積比（英: log area ratios; LAR）
arcsine of reflection coefficients; ASRC
impulse responses of LP synthesis filter; IR

係数表現によってノイズ耐性や計算量の特性が異なる。例えば音声符号化における線形予測（線形予測符号化）ではアナログ回線に由来するノイズへ耐性を持たせるために、LP coefficients 以外の係数表現がしばしば用いられてきた。

脚注

[脚注の使い方]

^ "最小二乗誤差推定による定式化" p.7 of 亀岡. (2014). 応用音響学第3回. 東京大学.
^ "The prediction coefficients $a_{k}$ are computed by first converting ... Bark-frequency cepstrum into a linear-frequency power spectral density ... then converted to an autocorrelation using an inverse FFT. From the auto-correlation, the Levinson-Durbin algorithm is used to compute the predictor" Valin & Skoglund. (2018). LPCNet: Improving Neural Speech Synthesis Through Linear Prediction. arxiv.
^ "Linear predictive coefficients (LP coefficients) have other representations: line spectral frequencies (LSF) ... etc" Tamanna Islam. (2000). Interpolation of Linear Prediction Coefficients for Speech Coding. Master thesis of Engineering. McGill University. p.20.

参考文献

G. U. Yule. On a method of investigating periodicities in disturbed series, with special reference to wolfer’s sunspot numbers. Phil. Trans. Roy. Soc., 226-A:267–298, 1927.
J. Makhoul. Linear prediction: A tutorial review. Proceedings of the IEEE, 63 (5):561–580, April 1975.
M. H. Hayes. Statistical Digital Signal Processing and Modeling. J. Wiley & Sons, Inc., New York, 1996.

外部リンク

PLP and RASTA (and MFCC, and inversion) in Matlab

[1] "最小二乗誤差推定による定式化" p.7 of 亀岡. (2014). 応用音響学第3回. 東京大学.

[2] "The prediction coefficients $a_{k}$ are computed by first converting ... Bark-frequency cepstrum into a linear-frequency power spectral density ... then converted to an autocorrelation using an inverse FFT. From the auto-correlation, the Levinson-Durbin algorithm is used to compute the predictor" Valin & Skoglund. (2018). LPCNet: Improving Neural Speech Synthesis Through Linear Prediction. arxiv.

[3] "Linear predictive coefficients (LP coefficients) have other representations: line spectral frequencies (LSF) ... etc" Tamanna Islam. (2000). Interpolation of Linear Prediction Coefficients for Speech Coding. Master thesis of Engineering. McGill University. p.20.

[1]

[2]

[3]

モデル