目次

・意気込み
・LPC分析
　- 簡単な解説
　- 実装
・LSP分析
　- 簡単な解説
　- 実装
・最小位相復元
・おしまい
・おまけ

意気込み

LPC分析はコードいくらでも落ちてるけどLSP分析は全然無いな〜と思い。
メモ代わりに。

あと、最後に、最小位相復元信号についても書く。
WORLDのパワースペクトル包絡から最小位相復元した時間信号を算出する。

LPC分析

簡単な解説

Linear Prediction Coding で LPC。線形予測分析。

Predictive とか Coefficient とかいろいろ書き方あるっぽい。

時間信号 $x_n$ を $x_n = - a_1x_{n-1} - a_2x_{n-2} - \dots - a_Kx_{n-K}$ で表現しようという試み。

このベクトル $\boldsymbol{A} = \{ a_k \}$ を K次の線形予測係数と呼ぶ。
※ $\boldsymbol{A}$ はK+1次のベクトルで、$A_0 = 1$ になる。

例えば、$N = 1000$ の信号の場合でも、$\boldsymbol{x}$ をたったK個の数字で表せちゃう最強の符号化モデルということ。

\begin{align}
L = \sum_n\left(x_n - \sum_k -a_kx_{n-k} \right)^2
\end{align}

を最小化する $\boldsymbol{A}$ を求めればOK.

信号と予測信号の二乗誤差を最小化すりゃOK.

実装

# 自己相関関数   order = K+1
def autocorr(x, order):
  N = len(x)
  r = np.zeros(order)
  for i in range(order):
    r[i] = np.sum(x[0:N-i] * x[i:N])
  
  return r

"""
レビンソンダービン法
LPC係数Aを、A = [A0,A1] -> [A0,A1,A2] -> [A0,A1,A2,A3]
のように、1個ずつ追加していく感じで求める
"""
def levinson(signal, order=10):
  R = autocorr(signal, order+1)  # 自己相関関数
  A = np.ones(2)  # LPC係数 A[0] == 1.
  
  A[1] = -R[1] / R[0]
  e = R[0] + R[1]*A[1]  # 残差
  
  for k in range(2, order+1):
    lam = -np.sum(A * R[k:0:-1]) / e
    U = np.hstack((A, 0))
    V = lam * U[::-1]
    a = U+V
    e = (1 - lam**2) * e

  return A, e

なんでこうなるかは「おまけ〜LPCについて」へ

LPC係数からスペクトル包絡を求めると・・・

N = 1024  # FFTサイズ
A, e = levinson(x, order=32)
H = np.fft.fft(A, N)
Amp = np.abs(H[0:int(N/2)+1])

env = 20*np.log10(np.sqrt(e)) - 20*np.log10(Amp)

なんでこうなるかは「おまけ〜LPCについて」へ

LSP分析

Line Spectral Pairs で LSP。線スペクトル対。

詳細とか数学的なことはよくわからないけど、
LPC係数をもっと素晴らしくした感じ。

K個のLSP式は、下記の2つの多項式の根（解）として定義され、
根（これが各周波数に相当）をLSP周波数と呼ぶ。

\begin{align}
P(z) &= A(z) + z^{-(K+1)}A(z^{-1}) \\
Q(z) &= A(z) - z^{-(K+1)}A(z^{-1})
\end{align}

"""
input : A  LPC係数（K+1 サイズ)
output : LSP分析した周波数
"""
def lpc2lsp(A):
  p = len(A)-1  # #LPC-order
  A1 = np.hstack((A, 0))
  A2 = A1[::-1]
  P1, Q1 = A1 + A2, A1 - A2

  if p%2: # if order is odd
    P_ = P1
    Q_, r = scipy.signal.deconvolve(Q1, [1,0,-1])
  else: # if order is even
    P_, r = scipy.signal.deconvolve(P1, [1, 1])
    Q_, r = scipy.signal.deconvolve(Q1, [1, -1])
  
  zP = np.roots(P_)  # 多項式P=0の解
  zQ = np.roots(Q_)  # 多項式Q=0の解

  omegaP = np.angle(zP[1::2])  # zPは複素数。その偏角が求めたいやつ
  omegaQ = np.angle(zQ[1::2])  # zQは複素数。その変革が求めたいやつ

  lsp = np.sort(np.hstack((-omegaP, -omegaQ)))
  return lsp / 2 / np.pi  # 角周波数を周波数に変換

なんでこうなるかは「おまけ〜LSPについて」へ

最小位相復元

最小位相を復元した実数ケプストラム。

〜手順〜
step 1. $\ \ \boldsymbol{x}_t$ の振幅スペクトル $S(\boldsymbol{\omega})$ のナイキスト成分を復元 → $S'(\boldsymbol{\omega})$
step 2. 対数とって逆フーリエ変換して実数部だけ取得（実数ケプストラム）
step 3. 実数ケプストラムに窓かけてFFTして、expに乗せて、iFFTして実部取得 → $\hat{\boldsymbol{x}}_t$

実装

"""
 入力 y : フレームtにおける振幅スペクトル
 出力 x : 最小位相復元した信号 
"""
def rceps(y):
  N = (len(y)-1)*2  # fft-size

  # ナイキスト周波数成分の振幅を復元
  y_recon = np.hstack((y, y[np.arange(len(y)-2, 0, -1)]))   # 1.
  
  ceps = np.fft.ifft(np.log(y_recon)).real  # 2.
  win = np.concatenate(([1.], 2.0*np.ones(int(N/2)-1), [1.], np.zeros(int(N/2) -1))) 

  x = np.fft.ifft(np.exp(np.fft.fft(win * ceps))).real  # 3.
  return x

ちなみにFFT($x$) と FFT($\hat{x}$) の振幅スペクトルは一致する。 LPC($x$) と LPC($\hat{x}$) は一致しないが、LPCスペクトル同士ではほぼ一致する

ということで、、、

CheapTrickでゲットしたスペクトル包絡を使って（最小位相復元した信号で）LPCやらLSP分析ができる。

素晴らしい。

実際にWORLDとコラボ

I = 300 # 3000番目のフレームでデモ
x, fs = soundfile('音声パス')
f0, time_axis = pyworld.harvest(x, fs)
sp = pw.cheaptrick(x, f0, time_axis, fs)
ap = sp[I] ** 0.5  # cheaptrickはパワースペクトルなので振幅スペクトルに

# WORLDの振幅スペクトルでLPC分析
x_hat = rceps(ap) # まずは時間領域の信号に変換（最小位相復元信号）
lpc, e = levinson(x_hat, order=32)
Afft = np.fft.fft(lpc, N)[0:int(N/2)+1]
logP = 20*np.log10(e**0.5) - 20 * np.log10(np.abs(Afft)+1e-10)

logP が最終的なLPCパワースペクトル。係数を10にすればLPC振幅スペクトル。お好きに。

WORLDの対数パワースペクトルとLPCの対数パワースペクトルがコレ↓

おしまい

Deep Learningも良いけどこういうのも良いよね
数学的な理論はくっそムズイです

おまけ

おまけ〜LPCについて

レビンソンダービンの解法は死ぬほどネットに転がってるので割愛。

自己回帰モデルの伝達関数は、周波数領域では

\begin{align}
H(z) &= \frac{\sigma}{A(z)} \\
\log|H(z)| &= \log|\sigma| - \log|A(z)|
\end{align}

らしい。

$\sigma$ は残差の標準偏差で、レビンソンダービン法で求めた e が残差分散だから、第一項は平方根とってあげればいい。
第二項はLPC係数の振幅スペクトルを求めればOK。

A, e = levinson(x, order=32)
H = np.fft.fft(A, N)
Amp = np.abs(H[0:int(N/2)+1])

env = 20*np.log10(np.sqrt(e)) - 20*np.log10(Amp)

よく、

H = scipy.signal.freqz(1, A)
H = scipy.signal.freqz(np.sqrt(e), A)

っていうのを見かけるけど、、
前者はさっきの式の分子が1、後者は今説明したやつと等価。

いろんな資料を見たけど
分子が残差だったり $\sigma$ だったりの場合も1の場合もよく見かける。
よくわからんｗ

おまけ〜LSPについて

\begin{align}
P(z) &= A(z) + z^{-(K+1)}A(z^{-1}) \\
Q(z) &= A(z) - z^{-(K+1)}A(z^{-1})
\end{align}

について、
次数が偶数の時、
$P(z)$ は $z + 1$ で割り切れて、$Q(z)$ は $z - 1$ で割り切れる（へぇ〜）。
そこで、

P_, r = scipy.signal.deconvolve(P1, [1, 1])  
Q_, r = scipy.signal.deconvolve(Q1, [1, -1])

X, r = scipy.signal.deconvolve(A, B)は、係数Aの多項式 ÷ 係数Bの多項式の商の係数X、余りの係数r っていう感じ

・1行目の返り直P_は、係数P1を、$1z + 1$ で割った商の係数を表す。
・2行目の返り直Q_は、係数Q1を、$1z - 1$ で割った商の係数を表す。
だから、
・$P(z)$ の根（$P(z) = 0$ の解）は、P_を係数とする多項式の根を
・$Q(z)$ の根（$Q(z) = 0$ の解）は、Q_を係数とする多項式の根を
求めれば良いので、

zP = np.roots(P_)  # 多項式P=0の解
zQ = np.roots(Q_)  # 多項式Q=0の解

になる。

次数が奇数の時、 $Q(z)$ は $z^2 - 1$ で割り切れる。だから、

P_ = P1
Q_, r = scipy.signal.deconvolve(Q1, [1,0,-1])

で、
・1行目は、係数P1のまま
・2行目の返り直Q_は、係数Q1を、$1z^2 + 0z - 1$ で割った商の係数を表す。

なんで $z^2 - 1$ で割り切れるかは・・・頭がいい人に聞いてください