PythonでLPC分析・LSP分析・最小位相復元とか

# 目次
・[意気込み](#意気込み)
・[LPC分析](#lpc分析)
　- [簡単な解説](#簡単な解説)
　- [実装](#実装)
・[LSP分析](#lsp分析)
　- [簡単な解説](#簡単な解説)
　- [実装](#実装)
・[最小位相復元](#最小位相復元)
・[おしまい](#おしまい)
・[おまけ](#おまけ)


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

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

# LPC分析
### 簡単な解説
Linear Prediction Coding で LPC。線形予測分析。
> Predictive とか Coefficient とかいろいろ書き方あるっぽい。

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

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

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

<br>

```math
\begin{align}
L = \sum_n\left(x_n - \sum_k -a_kx_{n-k} \right)^2
\end{align}
```
を最小化する $\boldsymbol{A}$ を求めればOK.
>信号と予測信号の二乗誤差を最小化すりゃOK.


### 実装
```python
# 自己相関関数   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
```

```python
"""
レビンソンダービン法
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について) へ

<br>

LPC係数からスペクトル包絡を求めると・・・
```python
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について」](#おまけlpcについて) へ

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/156754/eed96ed4-9b6b-e8f8-a86e-52e7e1af44d3.png)



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

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

K個のLSP式は、下記の2つの多項式の根（解）として定義され、
根（これが各周波数に相当）をLSP周波数と呼ぶ。
```math
\begin{align}
P(z) &= A(z) + z^{-(K+1)}A(z^{-1}) \\
Q(z) &= A(z) - z^{-(K+1)}A(z^{-1})
\end{align}
```

<br>

```python
"""
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について」](#おまけlspについて) へ








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

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


### 実装
```python
"""
 入力 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スペクトル同士ではほぼ一致する

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/156754/4bdfa6e1-501d-1e2f-71ee-c15354273b61.png)




<br>
<br>
ということで、、、

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

素晴らしい。




<br>
<br>

実際にWORLDとコラボ
```python
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の対数パワースペクトルがコレ↓
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/156754/3ed3c478-8144-f137-be59-3ed224335364.png)




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

# おまけ
### おまけ〜LPCについて
レビンソンダービンの解法は死ぬほどネットに転がってるので割愛。

<br>
<br>

自己回帰モデルの伝達関数は、周波数領域では
```math
\begin{align}
H(z) &= \frac{\sigma}{A(z)} \\
\log|H(z)| &= \log|\sigma| - \log|A(z)|
\end{align}
```
らしい。

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

よく、
```python
H = scipy.signal.freqz(1, A)
H = scipy.signal.freqz(np.sqrt(e), A)
```
っていうのを見かけるけど、、
前者はさっきの式の分子が1、後者は今説明したやつと等価。

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




### おまけ〜LSPについて
```math
\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$ で割り切れる（へぇ〜）。
そこで、
```python
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_を係数とする多項式の根を
求めれば良いので、
```python
zP = np.roots(P_)  # 多項式P=0の解
zQ = np.roots(Q_)  # 多項式Q=0の解
```
になる。

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

```python
P_ = P1
Q_, r = scipy.signal.deconvolve(Q1, [1,0,-1])
```
で、
・1行目は、係数P1のまま
・2行目の返り直Q_は、係数Q1を、$1z^2 + 0z - 1$ で割った商の係数を表す。
>なんで $z^2 - 1$ で割り切れるかは・・・ 頭がいい人に聞いてください