目次
・意気込み
・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$ で割り切れるかは・・・ 頭がいい人に聞いてください