線形予測符号化（LPC）の理論とPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lfilter


def autocorrelate(x, p):
    """
    ラグ 0 から p までの自己相関を計算する。
    バイアス推定（1/N 正規化なし）を返す。
    """
    N = len(x)
    n_fft = 2 ** int(np.ceil(np.log2(2 * N)))
    X = np.fft.fft(x, n=n_fft)
    r = np.fft.ifft(np.abs(X) ** 2).real[: p + 1]
    return r


def levinson_durbin(r, p):
    """
    Levinson-Durbin アルゴリズムで AR 係数を求める。

    Parameters
    ----------
    r : np.ndarray
        ラグ 0 から p までの自己相関 r[0], r[1], ..., r[p]
    p : int
        モデル次数

    Returns
    -------
    a : np.ndarray
        LPC 係数 [a_1, a_2, ..., a_p]（a_0 = 1 は含まない）
    E : float
        最小残差パワー
    k : np.ndarray
        反射係数（PARCOR）[k_1, ..., k_p]
    """
    a = np.zeros(p)
    k_arr = np.zeros(p)
    E = r[0]
    if E <= 0:
        return a, E, k_arr

    for m in range(p):
        # 反射係数 k_{m+1}
        acc = r[m + 1]
        for j in range(m):
            acc += a[j] * r[m - j]
        k = -acc / E

        # 係数更新（旧係数を退避してから更新）
        a_new = a.copy()
        a_new[m] = k
        for j in range(m):
            a_new[j] = a[j] + k * a[m - 1 - j]
        a = a_new

        # 残差パワー
        E = (1.0 - k * k) * E
        k_arr[m] = k

    return a, E, k_arr


def lpc(signal, p):
    """信号から LPC 係数を直接求める便利関数。"""
    r = autocorrelate(signal, p)
    return levinson_durbin(r, p)


# --- 検証: 既知の AR(4) プロセスを生成して係数を再推定 ---
np.random.seed(42)
true_a = np.array([1.0, -2.6895, 3.6076, -2.4801, 0.8546])  # A(z) = 1 + a_1 z^-1 + ...
print("A(z) の根の絶対値（安定性チェック）:", np.abs(np.roots(true_a)))

N = 4000
e = np.random.randn(N)                   # 白色励振
x = lfilter([1.0], true_a, e)            # AR(4) プロセスを生成

p = 4
a_est, E, k = lpc(x, p)

print("真の AR 係数 (a_1..a_p):", true_a[1:])
print("推定 AR 係数 (a_1..a_p):", a_est)
print(f"反射係数: {k}")
print(f"最小残差パワー: {E:.4f}")

A(z) の根の絶対値（安定性チェック）: [0.9539 0.9539 0.9691 0.9691]
真の AR 係数 (a_1..a_p): [-2.6895  3.6076 -2.4801  0.8546]
推定 AR 係数 (a_1..a_p): [-2.6233  3.4489 -2.3242  0.7926]
反射係数: [-0.754   0.9536 -0.6587  0.7926]
最小残差パワー: 5500.8131

import numpy as np
from scipy.linalg import solve_toeplitz


def lpc_toeplitz(signal, p):
    """
    Toeplitz ソルバを使った LPC 係数の計算（参考実装）。
    Levinson-Durbin と同じ解を返す。
    """
    r = autocorrelate(signal, p)
    # 左辺: Toeplitz(r[0..p-1]) , 右辺: -r[1..p]
    a = solve_toeplitz(r[:p], -r[1 : p + 1])
    E = r[0] + np.dot(a, r[1 : p + 1])
    return a, E


a_ref, E_ref = lpc_toeplitz(x, p)
print("Toeplitz ソルバ係数:", a_ref)
print("Levinson-Durbin との差:", np.max(np.abs(a_ref - a_est)))

Toeplitz ソルバ係数: [-2.6233  3.4489 -2.3242  0.7926]
Levinson-Durbin との差: 2.1316282072803006e-14

p_max = 20
Es_ar4 = []
for pp in range(1, p_max + 1):
    _, Epp, _ = lpc(x, pp)
    Es_ar4.append(Epp)

for pp, e1 in enumerate(Es_ar4, start=1):
    print(f"p={pp:2d}  E_p={e1:12.4f}")

p= 1  E_p=  288414.5319
p= 2  E_p=   26140.3432
p= 3  E_p=   14797.5360
p= 4  E_p=    5500.8131   <- 真の次数
p= 5  E_p=    5250.2230
p= 6  E_p=    5244.4969
p= 7  E_p=    5244.1644
p= 8  E_p=    5242.0517
p=10  E_p=    5236.2604
p=12  E_p=    5234.8054
p=20  E_p=    5232.3890

import matplotlib.pyplot as plt


def synth_vowel_inline(fs=8000, duration=0.04, f0=140, ar=None, noise=0.02, seed=0):
    np.random.seed(seed)
    t = np.arange(0, duration, 1 / fs)
    glottal = np.zeros_like(t)
    period = int(fs / f0)
    for k in range(0, len(t), period):
        L = min(30, len(t) - k)
        glottal[k : k + L] = np.exp(-np.arange(L) * 0.15)
    voice = lfilter([1.0], ar, glottal)
    voice = voice / (np.max(np.abs(voice)) + 1e-8)
    voice += noise * np.random.randn(len(voice))
    return voice


ar_true = [1.0, -2.6895, 3.6076, -2.4801, 0.8546]
vowel = synth_vowel_inline(ar=ar_true)

Es_vowel = []
for pp in range(1, p_max + 1):
    _, Epp, _ = lpc(vowel * np.hamming(len(vowel)), pp)
    Es_vowel.append(Epp)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
axes[0].plot(range(1, p_max + 1), Es_ar4, marker="o", color="steelblue")
axes[0].axvline(4, color="gray", linestyle=":", label="true order p=4")
axes[0].set_xlabel("Model order p")
axes[0].set_ylabel("Residual power $E_p$")
axes[0].set_title("AR(4) signal")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(range(1, p_max + 1), Es_vowel, marker="o", color="crimson")
axes[1].axvline(12, color="gray", linestyle=":", label="used order p=12")
axes[1].set_xlabel("Model order p")
axes[1].set_ylabel("Residual power $E_p$")
axes[1].set_title("Synthetic vowel")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("lpc_order_selection.png", dpi=150, bbox_inches="tight")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lfilter, freqz


def synth_vowel(fs=8000, duration=0.04, f0=140, ar=None, noise=0.02, seed=0):
    """声帯パルス励振 + AR フィルタの簡易合成音声。"""
    np.random.seed(seed)
    t = np.arange(0, duration, 1 / fs)
    glottal = np.zeros_like(t)
    period = int(fs / f0)
    for k in range(0, len(t), period):
        L = min(30, len(t) - k)
        glottal[k : k + L] = np.exp(-np.arange(L) * 0.15)
    voice = lfilter([1.0], ar, glottal)
    voice = voice / (np.max(np.abs(voice)) + 1e-8)
    voice += noise * np.random.randn(len(voice))
    return t, voice


# 母音「あ」相当の AR(4): フォルマント F1≈800, F2≈1200 Hz（A(z) の根は単位円内で安定）
fs = 8000
ar_true = [1.0, -2.6895, 3.6076, -2.4801, 0.8546]
t, voice = synth_vowel(fs=fs, duration=0.04, f0=140, ar=ar_true, noise=0.02)

# LPC 推定（窓掛けでスペクトルリーク低減）
window = np.hamming(len(voice))
a_est, E, _ = lpc(voice * window, p=12)
A = np.concatenate(([1.0], a_est))   # A(z) = 1 + a_1 z^-1 + ... + a_p z^-p
print("推定LPC係数 (p=12):", a_est)
print("残差パワー E:", E)

# 生スペクトル
n_fft = 1024
X = np.fft.rfft(voice * window, n=n_fft)
freqs = np.fft.rfftfreq(n_fft, d=1 / fs)
log_spec = 20 * np.log10(np.abs(X) + 1e-10)

# LPC スペクトル: H(e^{jw}) = sqrt(E) / A(e^{jw})
w, H = freqz([np.sqrt(E)], A, worN=n_fft // 2 + 1, fs=fs)
lpc_spec_db = 20 * np.log10(np.abs(H) + 1e-10)

fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(freqs, log_spec, color="steelblue", alpha=0.6, label="Log spectrum (windowed)")
ax.plot(w, lpc_spec_db, color="red", linewidth=2, label="LPC envelope (p=12)")
ax.set_ylim(-40, 40)
ax.axvline(820.5, color="gray", linestyle=":", linewidth=1)
ax.text(835, 35, "F1", color="gray", fontsize=10)
ax.axvline(1197.5, color="gray", linestyle=":", linewidth=1)
ax.text(1212, 35, "F2", color="gray", fontsize=10)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("Magnitude [dB]")
ax.set_title("LPC Spectral Envelope vs. Raw Spectrum")
ax.set_xlim(0, fs / 2)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("lpc_spectral_envelope.png", dpi=150, bbox_inches="tight")
plt.show()

推定LPC係数 (p=12): [-1.7831  0.9418  0.4716 -0.2637 -0.4911  0.2977  0.022   0.1982 -0.252
 -0.1395  0.2205 -0.0617]
残差パワー E: 0.4598260612610578

import numpy as np


def lpc_formants(a, fs, max_bw=400.0, min_freq=90.0):
    """
    LPC 係数からフォルマント（周波数と帯域幅）を抽出する。

    Parameters
    ----------
    a : np.ndarray
        LPC 係数 [a_1, ..., a_p]（a_0 = 1 は含まない）
    fs : float
        サンプリング周波数 [Hz]
    max_bw : float
        フォルマントとみなす帯域幅の上限 [Hz]
    min_freq : float
        フォルマントとみなす最低周波数 [Hz]

    Returns
    -------
    formants : list of (freq_Hz, bandwidth_Hz)
    """
    A = np.concatenate(([1.0], a))
    roots = np.roots(A)
    # 上半平面の根のみ（共役対の重複を避ける）
    roots = roots[np.imag(roots) > 0]

    formants = []
    for r in roots:
        theta = np.arctan2(np.imag(r), np.real(r))
        magnitude = np.abs(r)
        f = theta / (2 * np.pi) * fs
        bw = -np.log(max(magnitude, 1e-10)) / np.pi * fs
        if f >= min_freq and bw <= max_bw:
            formants.append((f, bw))

    return sorted(formants, key=lambda x: x[0])


formants = lpc_formants(a_est, fs)
for i, (f, bw) in enumerate(formants[:4], start=1):
    print(f"F{i}: {f:7.1f} Hz   (BW = {bw:5.1f} Hz)")

F1:   820.5 Hz   (BW =  98.1 Hz)
F2:  1197.5 Hz   (BW = 164.6 Hz)

import numpy as np


def lpc_to_cepstrum(a, n_cep):
    """
    LPC 係数 [a_1, ..., a_p] から LPC ケプストラム c[1..n_cep] を計算する。
    """
    p = len(a)
    c = np.zeros(n_cep)
    for m in range(1, n_cep + 1):
        if m <= p:
            c[m - 1] = -a[m - 1]
        # 再帰項
        s = 0.0
        for k_ in range(1, m):
            if (m - k_) <= p:
                s += (k_ / m) * c[k_ - 1] * a[m - k_ - 1]
        c[m - 1] -= s
    return c


cep = lpc_to_cepstrum(a_est, n_cep=20)
print("LPC ケプストラム c[1..10]:", cep[:10])

LPC ケプストラム c[1..10]: [ 1.78314  0.64802 -0.26099 -0.60069 -0.24648 -0.02235  0.25715  0.11109
  0.06295  0.1004 ]