線形予測符号化（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.2, 2.8, -1.8, 0.5])  # A(z) = 1 + a_1 z^-1 + ...
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}")

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)))

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
fs = 8000
ar_true = [1.0, -2.2, 2.8, -1.8, 0.5]
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

# 生スペクトル
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=f"LPC envelope (p=12)")
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.show()

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)")

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])