自己相関・相互相関の理論とPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate, correlation_lags

# --- テスト信号の生成 ---
np.random.seed(42)
fs = 1000      # サンプリング周波数 [Hz]
T = 2.0        # 信号長 [s]
t = np.arange(0, T, 1/fs)

# 信号1: 50Hzの正弦波 + ノイズ
x = np.sin(2 * np.pi * 50 * t) + 0.3 * np.random.randn(len(t))

# 信号2: xを100サンプル遅延させたもの + 異なるノイズ
delay_samples = 100
y = np.zeros_like(x)
y[delay_samples:] = x[:-delay_samples] + 0.3 * np.random.randn(len(t) - delay_samples)

# --- 自己相関 ---
def autocorr_fft(x):
    """FFTを用いた高速自己相関計算（正規化あり）"""
    N = len(x)
    # ゼロパディングして循環誤差を防ぐ
    X = np.fft.fft(x, n=2*N)
    # パワースペクトル → 逆FFT
    acf = np.fft.ifft(X * np.conj(X)).real
    # 正の遅延部分のみ取得し、正規化
    acf = acf[:N]
    acf /= acf[0]  # 正規化
    return acf

# --- 相互相関 ---
# scipy.signal.correlate は完全な相関を計算
corr_xy = correlate(x, y, mode='full')
lags = correlation_lags(len(x), len(y), mode='full')

# 正規化
norm = np.sqrt(np.sum(x**2) * np.sum(y**2))
corr_xy_normalized = corr_xy / norm

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

# 元の信号
axes[0].plot(t[:500], x[:500], label='x(t)', alpha=0.8)
axes[0].plot(t[:500], y[:500], label=f'y(t) = x(t - {delay_samples/fs*1000:.0f}ms)',
             alpha=0.8)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Input Signals')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 自己相関
acf = autocorr_fft(x)
lag_times = np.arange(len(acf)) / fs * 1000  # ms単位
axes[1].plot(lag_times, acf)
axes[1].axhline(y=0, color='k', linestyle='--', linewidth=0.8)
axes[1].set_xlabel('Lag [ms]')
axes[1].set_ylabel('Normalized ACF')
axes[1].set_xlim(0, 100)
axes[1].grid(True, alpha=0.3)
# 最初のピーク検出（ラグ0を除く）
first_peak_idx = np.argmax(acf[10:]) + 10
axes[1].set_title(f'Autocorrelation of x(t): argmax peak at '
                   f'{first_peak_idx/fs*1000:.0f} ms → {fs/first_peak_idx:.0f} Hz')
print(f"ACF最初のピーク: lag={first_peak_idx/fs*1000:.1f} ms → "
      f"周波数 = {fs/first_peak_idx:.1f} Hz")

# 相互相関
axes[2].plot(lags / fs * 1000, corr_xy_normalized)
axes[2].axhline(y=0, color='k', linestyle='--', linewidth=0.8)
axes[2].set_xlabel('Lag [ms]')
axes[2].set_ylabel('Normalized CCF')
axes[2].set_xlim(-300, 300)
axes[2].grid(True, alpha=0.3)

# 最大相関のラグを表示
peak_lag = lags[np.argmax(corr_xy_normalized)]
axes[2].set_title(f'Cross-Correlation (scipy raw): peak at lag = {peak_lag/fs*1000:.0f} ms')
print(f"相互相関ピーク: lag={peak_lag/fs*1000:.1f} ms → "
      f"遅延 = {peak_lag} サンプル")

plt.tight_layout()
plt.show()

ACF最初のピーク: lag=40.0 ms → 周波数 = 25.0 Hz
相互相関ピーク: lag=-100.0 ms → 遅延 = -100 サンプル

acf[20] = 0.8309  (真の基本周期のピーク)
acf[40] = 0.8315  (2周期目のピーク。ノイズによりわずかに上回る)

import numpy as np
import matplotlib.pyplot as plt

def detect_periodicity(signal, fs, max_lag_sec=None):
    """
    ACFを用いた周期性検出。

    Parameters
    ----------
    signal : np.ndarray
        入力信号
    fs : float
        サンプリング周波数 [Hz]
    max_lag_sec : float, optional
        解析する最大ラグ [s]。省略時は信号長の半分。

    Returns
    -------
    period_sec : float or None
        検出された周期 [s]（検出されない場合は None）
    acf : np.ndarray
        正規化ACF
    lags_sec : np.ndarray
        ラグ [s]
    """
    N = len(signal)
    max_lag = int(max_lag_sec * fs) if max_lag_sec else N // 2

    # FFTベースのACF
    X = np.fft.fft(signal, n=2*N)
    acf = np.fft.ifft(X * np.conj(X)).real[:N]
    acf = acf[:max_lag]
    acf_norm = acf / acf[0]
    lags_sec = np.arange(max_lag) / fs

    # 最初のピーク（1周期分先の極大）を検出
    # ACFはラグ0から減衰して zero_cross[0] で最初にゼロを切り、
    # その後 zero_cross[1] で再び正に戻る（この区間 [0]-[1] は「谷」であり
    # ピークではない点に注意）。真の周期に対応する最初のピークは
    # 谷を抜けた後の区間 zero_cross[1]-zero_cross[2] に現れる。
    zero_cross = np.where(np.diff(np.sign(acf_norm)))[0]
    if len(zero_cross) < 3:
        return None, acf_norm, lags_sec

    search_start = zero_cross[1] + 1
    search_end = zero_cross[2]
    peak_idx = search_start + np.argmax(acf_norm[search_start:search_end])
    period_sec = lags_sec[peak_idx]

    return period_sec, acf_norm, lags_sec


# --- 複数の信号で周期性を比較 ---
fs = 1000
t = np.arange(0, 1.0, 1/fs)
np.random.seed(0)

signals = {
    '50 Hz 正弦波': np.sin(2 * np.pi * 50 * t) + 0.2 * np.random.randn(len(t)),
    '複合正弦波\n(50 + 130 Hz)': (np.sin(2 * np.pi * 50 * t)
                                   + 0.5 * np.sin(2 * np.pi * 130 * t)
                                   + 0.1 * np.random.randn(len(t))),
    'ホワイトノイズ': np.random.randn(len(t)),
}

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

for i, (name, sig) in enumerate(signals.items()):
    period, acf, lags = detect_periodicity(sig, fs, max_lag_sec=0.15)
    if period:
        print(f"{name.strip()}: period={period*1000:.1f} ms → {1/period:.1f} Hz")
    else:
        print(f"{name.strip()}: 周期を検出できず")

    axes[i, 0].plot(t[:300], sig[:300])
    axes[i, 0].set_xlabel('Time [s]')
    axes[i, 0].set_ylabel('Amplitude')
    axes[i, 0].set_title(f'Signal: {name}')
    axes[i, 0].grid(True, alpha=0.3)

    axes[i, 1].plot(lags * 1000, acf)
    axes[i, 1].axhline(y=0, color='k', linestyle='--', linewidth=0.8)
    if period:
        axes[i, 1].axvline(x=period * 1000, color='r', linestyle='--',
                           label=f'Period={period*1000:.1f} ms\n→ {1/period:.1f} Hz')
        axes[i, 1].legend(fontsize=9)
    axes[i, 1].set_xlabel('Lag [ms]')
    axes[i, 1].set_ylabel('Normalized ACF')
    axes[i, 1].set_title('Autocorrelation')
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

50 Hz 正弦波: period=20.0 ms → 50.0 Hz
複合正弦波(50 + 130 Hz): period=22.0 ms → 45.5 Hz
ホワイトノイズ: period=5.0 ms → 200.0 Hz

import numpy as np
import matplotlib.pyplot as plt

def gcc_phat(x, y, fs, max_delay=None):
    """
    GCC-PHAT による遅延推定。

    Parameters
    ----------
    x, y : np.ndarray
        2つの入力信号
    fs : float
        サンプリング周波数 [Hz]
    max_delay : float, optional
        探索する最大遅延 [s]

    Returns
    -------
    delay_sec : float
        推定遅延 [s]（y は x より delay_sec 遅れている）
    cc : np.ndarray
        GCC-PHAT の相関値
    lags_sec : np.ndarray
        ラグ [s]
    """
    N = len(x) + len(y) - 1
    N_fft = 2 ** int(np.ceil(np.log2(N)))

    X = np.fft.fft(x, n=N_fft)
    Y = np.fft.fft(y, n=N_fft)

    # GCC-PHAT: 位相のみ利用（式(14)(16)と揃え、Xを共役にする）
    G = np.conj(X) * Y
    G_phat = G / (np.abs(G) + 1e-10)  # 数値安定化
    cc = np.fft.ifft(G_phat).real

    # lags の計算
    lags = np.fft.fftfreq(N_fft, d=1) * N_fft
    lags = lags.astype(int)

    # fftshift して中心をゼロラグに
    cc_shifted = np.fft.fftshift(cc)
    lags_shifted = np.fft.fftshift(lags)
    lags_sec = lags_shifted / fs

    if max_delay:
        mask = np.abs(lags_sec) <= max_delay
        cc_shifted = cc_shifted[mask]
        lags_sec = lags_sec[mask]

    peak_idx = np.argmax(cc_shifted)
    delay_sec = lags_sec[peak_idx]

    return delay_sec, cc_shifted, lags_sec


# --- テスト ---
fs = 16000
t = np.arange(0, 0.5, 1/fs)
np.random.seed(7)

# 音源信号
source = np.random.randn(len(t))
# マイク1: 直接波
mic1 = source.copy() + 0.1 * np.random.randn(len(t))
# マイク2: 実際の遅延（5 ms = 80 サンプル）
true_delay_samples = int(0.005 * fs)   # 80 サンプル
mic2 = np.zeros_like(mic1)
mic2[true_delay_samples:] = source[:-true_delay_samples] + 0.1 * np.random.randn(len(t) - true_delay_samples)

delay_est, cc, lags = gcc_phat(mic1, mic2, fs, max_delay=0.02)
print(f"真の遅延:     {true_delay_samples / fs * 1000:.2f} ms ({true_delay_samples} samples)")
print(f"GCC-PHAT推定: {delay_est * 1000:.2f} ms ({int(delay_est * fs)} samples)")

plt.figure(figsize=(10, 4))
plt.plot(lags * 1000, cc)
plt.axvline(x=delay_est * 1000, color='r', linestyle='--',
            label=f'Estimated delay = {delay_est*1000:.2f} ms')
plt.axvline(x=true_delay_samples / fs * 1000, color='g', linestyle=':',
            label=f'True delay = {true_delay_samples/fs*1000:.2f} ms')
plt.xlabel('Lag [ms]')
plt.ylabel('GCC-PHAT')
plt.title('GCC-PHAT Cross-Correlation for TDOA Estimation')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

真の遅延:     5.00 ms (80 samples)
GCC-PHAT推定: 5.00 ms (80 samples)

真の遅延:       5.00 ms (80 samples)
GCC-PHAT推定(逆規約): -5.00 ms (-80 samples)

import numpy as np
import matplotlib.pyplot as plt

def yin_pitch(signal, fs, fmin=75, fmax=500, threshold=0.1):
    """
    YIN アルゴリズムによるピッチ推定（簡易版）。

    Returns
    -------
    f0 : float or None
        推定F0 [Hz]
    d_norm : np.ndarray
        正規化差分関数
    """
    N = len(signal)
    W = N // 2  # 探索ウィンドウ

    # 差分関数（FFTで高速計算）
    X = np.fft.fft(signal, n=2*N)
    acf = np.fft.ifft(X * np.conj(X)).real[:N]

    d = np.zeros(W)
    d[0] = 0
    for l in range(1, W):
        d[l] = 2 * acf[0] - 2 * acf[l]

    # 累積平均正規化差分（CMND）
    d_norm = np.zeros(W)
    d_norm[0] = 1.0
    cumsum = 0.0
    for l in range(1, W):
        cumsum += d[l]
        d_norm[l] = d[l] / (cumsum / l)

    # ピッチ範囲の制限
    lag_min = int(fs / fmax)
    lag_max = int(fs / fmin)

    # 閾値を下回る最初のディップを検出
    f0 = None
    for l in range(lag_min, min(lag_max, W)):
        if d_norm[l] < threshold:
            # サブサンプル精度のためのパラボリック補間
            if 0 < l < W - 1:
                alpha = d_norm[l-1]
                beta = d_norm[l]
                gamma = d_norm[l+1]
                delta = (alpha - gamma) / (2 * (alpha - 2*beta + gamma))
                l_refined = l + delta
            else:
                l_refined = l
            f0 = fs / l_refined
            break

    return f0, d_norm


# --- テスト：異なるピッチの合成音声 ---
fs = 16000
t = np.arange(0, 0.1, 1/fs)  # 100 ms のフレーム

test_cases = {
    '男性声 (120 Hz)': 120,
    '女性声 (220 Hz)': 220,
    '子供声 (350 Hz)': 350,
}

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

for i, (label, f0_true) in enumerate(test_cases.items()):
    # ハーモニクスを含む合成音声
    np.random.seed(i)
    voice = (np.sin(2 * np.pi * f0_true * t)
             + 0.5 * np.sin(2 * np.pi * 2 * f0_true * t)
             + 0.3 * np.sin(2 * np.pi * 3 * f0_true * t)
             + 0.05 * np.random.randn(len(t)))

    f0_est, d_norm = yin_pitch(voice, fs)
    print(f"{label}: 真のF0={f0_true} Hz, 推定F0={f0_est:.2f} Hz, "
          f"誤差={f0_est - f0_true:+.2f} Hz ({(f0_est - f0_true) / f0_true * 100:+.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'{label}: True F0 = {f0_true} Hz')
    axes[i, 0].grid(True, alpha=0.3)

    lag_ms = np.arange(len(d_norm)) / fs * 1000
    axes[i, 1].plot(lag_ms, d_norm)
    axes[i, 1].axhline(y=0.1, color='r', linestyle='--', label='threshold=0.1')
    if f0_est:
        axes[i, 1].axvline(x=1000/f0_est, color='g', linestyle='--',
                           label=f'Est. F0 = {f0_est:.1f} Hz')
    axes[i, 1].set_xlabel('Lag [ms]')
    axes[i, 1].set_ylabel('Normalized d[l]')
    axes[i, 1].set_title('YIN: Cumulative Mean Normalized Difference')
    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()

男性声 (120 Hz): 真のF0=120 Hz, 推定F0=120.09 Hz, 誤差=+0.09 Hz (+0.07%)
女性声 (220 Hz): 真のF0=220 Hz, 推定F0=218.98 Hz, 誤差=-1.02 Hz (-0.47%)
子供声 (350 Hz): 真のF0=350 Hz, 推定F0=348.34 Hz, 誤差=-1.66 Hz (-0.48%)

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# --- AR(2) プロセスのシミュレーション ---
# x[t] = 0.7 * x[t-1] - 0.3 * x[t-2] + noise
np.random.seed(42)
N = 300
phi = [0.7, -0.3]   # AR係数
noise_std = 1.0

x = np.zeros(N)
for t in range(2, N):
    x[t] = phi[0] * x[t-1] + phi[1] * x[t-2] + np.random.randn() * noise_std

# ACF と PACF の計算
nlags = 20
acf_vals = acf(x, nlags=nlags, fft=True)
pacf_vals = pacf(x, nlags=nlags)

# 信頼区間（近似: ±1.96/√N）
ci = 1.96 / np.sqrt(N)

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

# 時系列プロット
axes[0].plot(x)
axes[0].set_xlabel('Time step')
axes[0].set_ylabel('x[t]')
axes[0].set_title('AR(2) Process: x[t] = 0.7·x[t-1] − 0.3·x[t-2] + noise')
axes[0].grid(True, alpha=0.3)

# ACF プロット
lags = np.arange(nlags + 1)
axes[1].bar(lags, acf_vals, color='steelblue', alpha=0.7)
axes[1].axhline(y=ci, color='r', linestyle='--', linewidth=0.8, label='±95% CI')
axes[1].axhline(y=-ci, color='r', linestyle='--', linewidth=0.8)
axes[1].axhline(y=0, color='k', linewidth=0.5)
axes[1].set_xlabel('Lag')
axes[1].set_ylabel('ACF')
axes[1].set_title('Autocorrelation Function (ACF): Gradual decay → AR process')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# PACF プロット
axes[2].bar(lags, pacf_vals, color='coral', alpha=0.7)
axes[2].axhline(y=ci, color='r', linestyle='--', linewidth=0.8, label='±95% CI')
axes[2].axhline(y=-ci, color='r', linestyle='--', linewidth=0.8)
axes[2].axhline(y=0, color='k', linewidth=0.5)
axes[2].set_xlabel('Lag')
axes[2].set_ylabel('PACF')
axes[2].set_title('Partial ACF (PACF): Cuts off after lag 2 → confirms AR(2)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"PACF(1) = {pacf_vals[1]:.3f} (真値: φ₁ = {phi[0]})")
print(f"PACF(2) = {pacf_vals[2]:.3f} (真値: φ₂ = {phi[1]})")
print(f"PACF(3) = {pacf_vals[3]:.3f} (カットオフ: 0 に近い)")

PACF(1) = 0.486 (真値: φ₁ = 0.7)
PACF(2) = -0.291 (真値: φ₂ = -0.3)
PACF(3) = -0.031 (カットオフ: 0 に近い)