ケプストラム分析の理論とPython実装

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

def compute_real_cepstrum(signal, n_fft=None):
    """
    実数ケプストラムの計算。

    Parameters
    ----------
    signal : np.ndarray
        入力信号
    n_fft : int, optional
        FFTの点数（省略時は信号長の次の2の累乗）

    Returns
    -------
    cepstrum : np.ndarray
        実数ケプストラム
    """
    if n_fft is None:
        n_fft = 2 ** int(np.ceil(np.log2(len(signal))))
    # ステップ1: DFT
    X = np.fft.fft(signal, n=n_fft)
    # ステップ2: 振幅スペクトルの対数
    log_spec = np.log(np.abs(X) + 1e-10)
    # ステップ3: 逆DFT（実部のみ取得）
    cepstrum = np.fft.ifft(log_spec).real
    return cepstrum


def formant_ar_coeffs(formants_bw, fs):
    """
    フォルマント周波数・帯域幅から安定な AR フィルタ係数を設計する。

    各フォルマント (F, BW) を2次共振器（複素共役極ペア）としてモデル化し、
    極の半径 r = exp(-pi*BW/fs) を必ず1未満に保つことで、lfilter が発散
    しない安定フィルタを保証する。AR係数を天下り的に与えると容易に不安定
    （極が単位円の外）になりうる点は「実装上のエッジケースと注意点」の節
    で検証する。
    """
    poly = np.array([1.0])
    for F, BW in formants_bw:
        r = np.exp(-np.pi * BW / fs)
        theta = 2 * np.pi * F / fs
        second_order = np.array([1.0, -2 * r * np.cos(theta), r**2])
        poly = np.convolve(poly, second_order)
    return poly


# --- 合成音声信号の生成 ---
fs = 8000
T = 0.05      # 50 ms のフレーム
t = np.arange(0, T, 1/fs)
f0 = 150      # 基本周波数 [Hz]

# 声帯パルス列（グロッタルパルスの近似）
glottal = np.zeros_like(t)
period_samples = int(fs / f0)
for i in range(0, len(t), period_samples):
    if i + 30 < len(t):
        glottal[i:i+30] = np.exp(-np.arange(30) * 0.15)

# 声道フィルタ: フォルマント F1=800Hz(帯域幅80Hz)・F2=1200Hz(帯域幅100Hz) を持つ AR(4) モデル
ar_coeffs = formant_ar_coeffs([(800, 80), (1200, 100)], fs)
print("AR係数:", np.round(ar_coeffs, 4))
print("極の絶対値（全て1未満なら安定）:", np.round(np.abs(np.roots(ar_coeffs)), 4))

voiced = lfilter([1.0], ar_coeffs, glottal)
voiced = voiced / np.max(np.abs(voiced) + 1e-8)

# --- ケプストラム計算 ---
n_fft = 512
cepstrum = compute_real_cepstrum(voiced, n_fft=n_fft)
quefrency = np.arange(n_fft) / fs * 1000   # ms 単位

# 対数スペクトル
X = np.fft.fft(voiced, n=n_fft)
log_spec = np.log(np.abs(X[:n_fft//2 + 1]) + 1e-10)
freqs = np.fft.rfftfreq(n_fft, 1/fs)

# ケフレンシー1〜20msの範囲でのピーク上位5件
idx_range = np.where((quefrency[:n_fft//2] > 1) & (quefrency[:n_fft//2] < 20))[0]
top5 = idx_range[np.argsort(cepstrum[idx_range])[::-1][:5]]
for idx in top5:
    print(f"q={quefrency[idx]:.3f} ms, cepstrum={cepstrum[idx]:.5f}")

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# 時間波形
axes[0].plot(t * 1000, voiced)
axes[0].set_xlabel('Time [ms]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title(f'Synthetic Voice: F0={f0} Hz, Formants F1≈800Hz, F2≈1200Hz')
axes[0].grid(True, alpha=0.3)

# 対数振幅スペクトル
axes[1].plot(freqs, log_spec)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Log Magnitude')
axes[1].set_title('Log Amplitude Spectrum')
axes[1].set_xlim(0, 4000)
axes[1].grid(True, alpha=0.3)

# ケプストラム
axes[2].plot(quefrency[:n_fft//2], cepstrum[:n_fft//2])
axes[2].axvline(x=1000/f0, color='r', linestyle='--',
                label=f'1/F0 = {1000/f0:.1f} ms')
axes[2].set_xlabel('Quefrency [ms]')
axes[2].set_ylabel('Cepstral Amplitude')
axes[2].set_title('Real Cepstrum: Pitch peak visible at 1/F0')
axes[2].set_xlim(0, 20)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

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

def lifter_cepstrum(cepstrum, cutoff_ms, fs, n_fft):
    """
    ローパスリフタリングによるスペクトル包絡抽出。

    Parameters
    ----------
    cepstrum : np.ndarray
        入力ケプストラム（長さ n_fft）
    cutoff_ms : float
        ローパスリフタのカットオフ [ms]
    fs : float
        サンプリング周波数
    n_fft : int
        FFTの点数

    Returns
    -------
    envelope_log : np.ndarray
        対数スペクトル包絡（長さ n_fft//2 + 1）
    """
    cutoff_samples = int(cutoff_ms * fs / 1000)

    # ローパスリフタ
    lifter = np.zeros(n_fft)
    lifter[:cutoff_samples] = 1.0
    lifter[-cutoff_samples+1:] = 1.0  # 対称成分

    # リフタリング
    liftered = cepstrum * lifter

    # 逆DFT → 対数スペクトル包絡
    envelope_log = np.fft.fft(liftered, n=n_fft).real
    return envelope_log[:n_fft//2 + 1]


# --- 上のコードに続けて実行 ---
cutoff_list = [2.0, 3.5, 5.0]   # ms 単位のカットオフ

fig, axes = plt.subplots(len(cutoff_list), 1, figsize=(12, 10))

for ax, cutoff_ms in zip(axes, cutoff_list):
    envelope = lifter_cepstrum(cepstrum, cutoff_ms, fs, n_fft)

    ax.plot(freqs, log_spec, alpha=0.5, label='Log spectrum', color='steelblue')
    ax.plot(freqs, envelope, linewidth=2, label=f'Envelope (cutoff={cutoff_ms} ms)',
            color='red')
    ax.set_xlabel('Frequency [Hz]')
    ax.set_ylabel('Log Magnitude')
    ax.set_title(f'Spectral Envelope Extraction via Low-pass Liftering ({cutoff_ms} ms)')
    ax.set_xlim(0, 4000)
    ax.legend()
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

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

def formant_ar_coeffs(formants_bw, fs):
    """フォルマント周波数・帯域幅から安定な AR フィルタ係数を設計する（既出のヘルパー関数）。"""
    poly = np.array([1.0])
    for F, BW in formants_bw:
        r = np.exp(-np.pi * BW / fs)
        theta = 2 * np.pi * F / fs
        second_order = np.array([1.0, -2 * r * np.cos(theta), r**2])
        poly = np.convolve(poly, second_order)
    return poly


def cepstrum_pitch_detection(signal, fs, fmin=60, fmax=600):
    """
    ケプストラムによるピッチ（F0）検出。

    Parameters
    ----------
    signal : np.ndarray
        入力音声フレーム
    fs : float
        サンプリング周波数 [Hz]
    fmin, fmax : float
        探索するF0の範囲 [Hz]

    Returns
    -------
    f0 : float or None
        推定F0 [Hz]
    cepstrum : np.ndarray
        ケプストラム
    quefrency_ms : np.ndarray
        ケフレンシー [ms]
    """
    n_fft = 2 ** int(np.ceil(np.log2(len(signal))))

    # ケプストラム計算
    X = np.fft.fft(signal * np.hanning(len(signal)), n=n_fft)
    log_spec = np.log(np.abs(X) + 1e-10)
    cepstrum = np.fft.ifft(log_spec).real

    # ケフレンシーの範囲を F0 範囲に変換
    q_min = int(fs / fmax)
    q_max = int(fs / fmin)
    q_max = min(q_max, n_fft // 2)

    # ピーク検出
    pitch_region = cepstrum[q_min:q_max]
    peak_idx = np.argmax(pitch_region) + q_min
    f0 = fs / peak_idx

    quefrency_ms = np.arange(n_fft) / fs * 1000

    return f0, cepstrum, quefrency_ms


# --- 異なるF0での推定テスト ---
fs = 16000
T = 0.04  # 40 ms フレーム
t = np.arange(0, T, 1/fs)

test_f0s = [120, 180, 250]  # Hz

# 声道フィルタ: フォルマント F1=700Hz(帯域幅80Hz)・F2=1300Hz(帯域幅120Hz) を持つ安定な AR(4) モデル
ar = formant_ar_coeffs([(700, 80), (1300, 120)], fs)

fig, axes = plt.subplots(len(test_f0s), 2, figsize=(14, 10))

for i, f0_true in enumerate(test_f0s):
    np.random.seed(i)
    # 合成音声生成
    period = int(fs / f0_true)
    glottal = np.zeros_like(t)
    for k in range(0, len(t), period):
        if k + 20 < len(t):
            glottal[k:k+20] = np.exp(-np.arange(20) * 0.2)

    voice = lfilter([1.0], ar, glottal)
    voice = voice / (np.max(np.abs(voice)) + 1e-8)
    voice += 0.02 * np.random.randn(len(t))

    f0_est, cep, q_ms = cepstrum_pitch_detection(voice, fs)
    n_fft = len(cep)
    err_pct = abs(f0_est - f0_true) / f0_true * 100
    print(f"True F0={f0_true} Hz -> Estimated F0={f0_est:.2f} Hz (error {err_pct:.2f}%)")

    axes[i, 0].plot(t * 1000, voice)
    axes[i, 0].set_xlabel('Time [ms]')
    axes[i, 0].set_ylabel('Amplitude')
    axes[i, 0].set_title(f'Voice Frame: True F0 = {f0_true} Hz')
    axes[i, 0].grid(True, alpha=0.3)

    axes[i, 1].plot(q_ms[:n_fft//4], cep[:n_fft//4])
    axes[i, 1].axvline(x=1000/f0_true, color='g', linestyle=':',
                       label=f'True 1/F0 = {1000/f0_true:.1f} ms')
    axes[i, 1].axvline(x=1000/f0_est, color='r', linestyle='--',
                       label=f'Estimated F0 = {f0_est:.1f} Hz')
    axes[i, 1].set_xlabel('Quefrency [ms]')
    axes[i, 1].set_ylabel('Cepstral Amplitude')
    axes[i, 1].set_title(f'Cepstrum Pitch Detection: Est={f0_est:.1f} Hz')
    axes[i, 1].set_xlim(0, 30)
    axes[i, 1].legend(fontsize=9)
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def hz_to_mel(hz):
    """Hzをメルに変換"""
    return 2595 * np.log10(1 + hz / 700)

def mel_to_hz(mel):
    """メルをHzに変換"""
    return 700 * (10 ** (mel / 2595) - 1)

def compute_mfcc(signal, fs, n_mfcc=13, n_filters=26, n_fft=512,
                 hop_length=160, fmin=80, fmax=None):
    """
    MFCCの計算。

    Parameters
    ----------
    signal : np.ndarray
        入力音声信号
    fs : float
        サンプリング周波数 [Hz]
    n_mfcc : int
        MFCC次元数
    n_filters : int
        メルフィルタバンクの数
    n_fft : int
        FFTの点数
    hop_length : int
        フレームシフト（サンプル数）
    fmin, fmax : float
        メルフィルタバンクの周波数範囲 [Hz]

    Returns
    -------
    mfcc : np.ndarray
        MFCC行列 (n_mfcc, n_frames)
    """
    if fmax is None:
        fmax = fs / 2

    N = len(signal)
    n_frames = 1 + (N - n_fft) // hop_length

    # --- プリエンファシス ---
    pre_emphasis = 0.97
    signal_pe = np.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

    # --- メルフィルタバンクの構築 ---
    mel_min = hz_to_mel(fmin)
    mel_max = hz_to_mel(fmax)
    mel_points = np.linspace(mel_min, mel_max, n_filters + 2)
    hz_points = mel_to_hz(mel_points)
    bin_points = np.floor((n_fft + 1) * hz_points / fs).astype(int)

    filter_bank = np.zeros((n_filters, n_fft // 2 + 1))
    for m in range(1, n_filters + 1):
        for k in range(bin_points[m-1], bin_points[m]):
            filter_bank[m-1, k] = (k - bin_points[m-1]) / (bin_points[m] - bin_points[m-1])
        for k in range(bin_points[m], bin_points[m+1] + 1):
            filter_bank[m-1, k] = (bin_points[m+1] - k) / (bin_points[m+1] - bin_points[m])

    # --- フレームごとの処理 ---
    window = np.hamming(n_fft)
    mfcc = np.zeros((n_mfcc, n_frames))

    for i in range(n_frames):
        start = i * hop_length
        frame = signal_pe[start:start + n_fft]
        if len(frame) < n_fft:
            frame = np.pad(frame, (0, n_fft - len(frame)))

        # パワースペクトル
        power_spec = np.abs(np.fft.rfft(frame * window)) ** 2

        # メルフィルタバンク適用
        mel_energy = np.dot(filter_bank, power_spec)
        log_mel_energy = np.log(mel_energy + 1e-10)

        # DCT（MFCCの計算）
        K = n_filters
        for m in range(n_mfcc):
            mfcc[m, i] = np.sum(
                log_mel_energy * np.cos(m * (np.arange(K) + 0.5) * np.pi / K)
            )

    return mfcc


# --- テスト信号でMFCCを計算 ---
fs = 16000
duration = 1.5
t = np.arange(0, duration, 1/fs)
np.random.seed(42)

# 母音的な信号（フォルマント変化あり）
from scipy.signal import lfilter, butter

def make_vowel(t, fs, f0, f1, f2, noise_std=0.01):
    """フォルマントを持つ合成音声"""
    period = int(fs / f0)
    glottal = np.zeros(len(t))
    for k in range(0, len(t), period):
        if k + 25 < len(t):
            glottal[k:k+25] = np.exp(-np.arange(25) * 0.15)
    # バンドパスフィルタでフォルマント模擬
    def bp(f_c, bw, sig):
        b, a = butter(2, [f_c - bw/2, f_c + bw/2], btype='band', fs=fs)
        return lfilter(b, a, sig)
    v = bp(f1, 200, glottal) + 0.6 * bp(f2, 200, glottal)
    v /= np.max(np.abs(v)) + 1e-8
    v += noise_std * np.random.randn(len(v))
    return v

# 3段階で変化する信号
t1 = t[t < 0.5]
t2 = t[(t >= 0.5) & (t < 1.0)]
t3 = t[t >= 1.0]
v1 = make_vowel(t1, fs, f0=120, f1=730, f2=1090)   # 「あ」
v2 = make_vowel(t2, fs, f0=120, f1=270, f2=2290)    # 「い」
v3 = make_vowel(t3, fs, f0=120, f1=400, f2=800)     # 「う」
speech = np.concatenate([v1, v2, v3])

# MFCC計算
mfcc = compute_mfcc(speech, fs, n_mfcc=13, n_fft=512, hop_length=160)
n_frames = mfcc.shape[1]
time_axis = np.arange(n_frames) * 160 / fs
print("MFCC shape:", mfcc.shape)
print(f"MFCC value range: [{mfcc.min():.2f}, {mfcc.max():.2f}]")

# プロット
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

axes[0].plot(t[:len(speech)], speech, linewidth=0.5)
axes[0].axvline(x=0.5, color='r', linestyle='--', alpha=0.7)
axes[0].axvline(x=1.0, color='r', linestyle='--', alpha=0.7)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Speech Signal: Vowel "a"→"i"→"u" Transition')
axes[0].grid(True, alpha=0.3)

img = axes[1].imshow(mfcc, aspect='auto', origin='lower',
                     extent=[0, time_axis[-1], 0, 13],
                     cmap='RdBu_r', interpolation='nearest')
axes[1].axvline(x=0.5, color='y', linestyle='--', alpha=0.8)
axes[1].axvline(x=1.0, color='y', linestyle='--', alpha=0.8)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('MFCC coefficient index')
axes[1].set_title('MFCC: Changes with vowel transitions visible')
plt.colorbar(img, ax=axes[1], label='Coefficient value')

# 各時刻のMFCCベクトル比較
frame_a = int(0.25 * fs / 160)
frame_i = int(0.75 * fs / 160)
frame_u = int(1.25 * fs / 160)
d_ai = np.linalg.norm(mfcc[:, frame_a] - mfcc[:, frame_i])
d_au = np.linalg.norm(mfcc[:, frame_a] - mfcc[:, frame_u])
d_iu = np.linalg.norm(mfcc[:, frame_i] - mfcc[:, frame_u])
print(f"Euclidean distance a-i: {d_ai:.2f}, a-u: {d_au:.2f}, i-u: {d_iu:.2f}")
axes[2].plot(range(1, 14), mfcc[:, frame_a], 'o-', label='"a" (t=0.25s)')
axes[2].plot(range(1, 14), mfcc[:, frame_i], 's-', label='"i" (t=0.75s)')
axes[2].plot(range(1, 14), mfcc[:, frame_u], '^-', label='"u" (t=1.25s)')
axes[2].set_xlabel('MFCC index')
axes[2].set_ylabel('Coefficient value')
axes[2].set_title('MFCC vectors for different vowels: distinguishable patterns')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()