音声信号処理の Python 実現ガイド

第1章 音声信号処理の基本

本章では、音声信号処理の基礎概念を説明する。これらは後続の用語リストとプログラム例を理解する前提知識となる。

• アナログ信号・デジタル信号：アナログ信号は時間的に連続した値をとり、デジタル信号は離散的な時刻に離散的な値をとる。音声処理ではアナログ信号をデジタル信号に変換（A/D変換：Analog to Digital変換）して処理する。

• 周波数：信号の周期的な変化の速さを表す物理量である。単位はHz（ヘルツ）で、1秒あたりの振動回数を示す。

• 時間領域・周波数領域：信号の2つの表現方法である。時間領域は振幅の時間変化を表し、周波数領域は信号に含まれる周波数成分の分布を表す。

• フーリエ変換：時間領域の信号を周波数領域に変換する数学的手法である。信号に含まれる周波数成分とその強度を抽出できる。逆フーリエ変換により周波数領域から時間領域へ戻すことができる。

• 周波数分解能・時間分解能：信号分析における精度の指標であり、両者はトレードオフの関係にある。周波数分解能は周波数の識別精度、時間分解能は時間変化の追跡精度を示す。窓長を長くすると周波数分解能は高くなり、時間分解能は低下する。

• スペクトル：信号の周波数成分とその強度の分布である。横軸に周波数、縦軸に振幅またはパワーをとったグラフで表す。

• 周期性：信号が一定の時間間隔で繰り返すパターンを持つ性質である。音声では有声音が周期性を持つ。

• 定常性・非定常性：信号の統計的性質が時間によって変化するか否かを表す。定常信号は統計的性質が一定であり、非定常信号は時間とともに変化する。音声信号や環境音は非定常信号である。

• ノイズ：目的とする信号に混入する不要な成分である。信号対雑音比（SNR：Signal-to-Noise Ratio）で信号品質を評価する。

• 時間窓：信号を短い区間に分割するための時間範囲である。短時間フーリエ変換では時間窓を移動させながら各区間で周波数分析を行う。

第2章 用語リスト

本章では、音声信号処理の専門用語を分野別に整理する。用語は、音声信号の基本概念、周波数領域分析、スペクトル特徴量、基本周波数推定、フォルマント分析、時間領域特徴量、リズム・時間構造特徴量、環境音解析、信号処理技術、評価・検証手法の10分野に分類している。各用語に対応するPython関数は「Pythonでの実現」章の表で示し、コード例は各演習の章に示す。

音声信号の基本概念

• サンプリング周波数：アナログ信号をデジタル信号に変換する際の1秒あたりのサンプル数である。単位はHz。音声では44.1kHzや48kHzが用いられる。

• ビット深度：各サンプルの振幅を表現するビット数である。16ビットや24ビットが用いられる。ビット深度が大きいほど振幅を細かく表現でき、ダイナミックレンジ（表現可能な音量範囲）が広がる。

• チャンネル数：音声信号の独立した信号系列の数である。モノラル（1チャンネル）、ステレオ（2チャンネル）、マルチチャンネル（5.1チャンネルなど）がある。

• 振幅：音声信号の瞬時的な大きさを表す値である。

• 時間長：音声信号の継続時間である。サンプル数をサンプリング周波数で除算して求める。

周波数領域分析

• スペクトログラム：音声信号の時間-周波数表現である。横軸に時間、縦軸に周波数、色の濃淡でパワーを表示する。

• 短時間フーリエ変換（STFT：Short-Time Fourier Transform）：信号を短い時間窓で区切り、各区間でフーリエ変換を適用する手法である。非定常信号の周波数成分の時間変化を分析する。

• 窓関数：信号を短い区間に分割する際に適用する関数である。ハン窓、ハミング窓、ブラックマン窓などがある。窓関数の選択は、メインローブ幅とサイドローブによるスペクトルの漏れに影響する。

• フレーム長：STFTで信号を分割する各区間の長さ（サンプル数）である。長いほど周波数分解能が高くなるが、時間分解能は低下する。

• ホップ長：STFTで連続するフレーム間の移動量である。フレーム長より短く設定してフレーム間のオーバーラップ（重複）を実現する。

• オーバーラップ率：連続するフレーム間で重複する部分の割合である。50%や75%が用いられる。

• FFTサイズ：高速フーリエ変換（FFT：Fast Fourier Transform）で使用するサンプル数である。2のべき乗（512、1024、2048など）が用いられる。窓長とともに大きくすると周波数分解能が高くなる。窓長を変えずにゼロ埋めで大きくした場合は、周波数軸の点の間隔が細かくなるだけで、分解能そのものは窓長で決まる。

• パワースペクトル密度（PSD：Power Spectral Density）：単位周波数あたりのパワーを表す指標である。

• スペクトル包絡：スペクトルの全体的な形状を表す曲線である。声道の共鳴特性を反映し、音色の特徴を表す。フォルマント推定ではLPCなどで平滑化した包絡を扱う。

スペクトル特徴量

スペクトル特徴量は、音声の音色や質感を数値化する指標である。音声認識や音楽情報検索で音の特徴を定量的に比較する際に用いる。

• スペクトル重心：スペクトルの重心周波数である。値が高いほど高周波成分が多い。

• スペクトル帯域幅：スペクトル重心を中心とした周波数の広がりである。

• スペクトルロールオフ：全パワーの一定割合（85%や95%）が含まれる周波数上限である。高周波成分の広がりを示す。

• スペクトル平坦度：スペクトルの平坦さを表す指標である。幾何平均と算術平均の比で算出する。1に近いほどエネルギーが全周波数に分散し、0に近いほど特定周波数に集中している状態を示す。

基本周波数推定

基本周波数（F0）は音声の高さを決定するパラメータである。音声の抑揚分析や歌声合成で用いる。

• 基本周波数（F0）：音声の周期的な振動の基本となる周波数である。会話における人間の声では、男性で約85-155Hz、女性で約165-255Hz程度である。

• 有声音：声帯の周期的な振動により生成される音声である。母音や有声子音が該当し、明確な基本周波数を持つ。

• 無声音：声帯の振動を伴わず、気流の乱れにより生成される音声である。無声子音（/s/、/k/など）が該当し、基本周波数を持たない。

• 自己相関法：信号と時間シフトした自身との相関を計算し、周期性を検出する手法である。ラグ0近傍ではなく、F0の探索範囲に対応するラグ範囲内の局所ピークを用いる。

• ケプストラム法：対数パワースペクトルに逆フーリエ変換を適用して周期性を検出する手法である。

• PYIN：確率的YIN法である。各時刻での基本周波数の確率分布を推定する手法であり、librosaライブラリで実装されている。

• オクターブ誤り：基本周波数推定で、真の基本周波数の整数倍または整数分の1の値を誤って推定する現象である。探索範囲の指定と時間方向の連続性の確認により軽減できる。

フォルマント分析

フォルマントは母音の音色を特徴づけるパラメータである。音声認識や音声合成で用いる。

• フォルマント：声道の共鳴により生じるスペクトル包絡のピーク周波数である。低い側から第1フォルマント（F1）、第2フォルマント（F2）と呼ぶ。F1は口の開きの程度、F2は舌の前後位置に関係する。

• 声道：声帯から唇までの音声生成経路である。共鳴管として機能し、音源スペクトルを変形させる。

• 線形予測分析（LPC：Linear Predictive Coding）：過去のサンプル値の線形結合で現在のサンプル値を予測するモデルである。声道特性のモデル化やフォルマント推定に用いる。

• LPC係数：線形予測分析で得られる予測係数である。声道の共鳴特性を表す。

• 予測次数：線形予測分析で使用する過去のサンプル数である。実音声では、サンプリング周波数（kHz単位）に2〜4を加えた程度（16kHzなら18前後）に設定するのが目安である。

• バンド幅：フォルマントの周波数幅である。共鳴のシャープさを表し、減衰の速さに関係する。LPCの極の半径を用いて計算できる。

• のこぎり波：基本周波数とその整数倍の倍音を含む周期信号である。有声音の音源モデルとして用いる。

時間領域特徴量

時間領域特徴量は、信号の波形から直接計算する指標である。周波数分析を必要としないため計算負荷が小さい。

• RMSエネルギー：信号の二乗平均平方根（Root Mean Square）である。音声の強度を表す。

• ゼロ交差率（ZCR：Zero-Crossing Rate）：信号が0を横切る回数の割合である。高周波成分の多さを示し、有声音と無声音の判別に用いる。

• クレストファクター：信号のピーク値とRMS値の比である。波高率を表し、インパルス性（突発的な大振幅成分の強さ）を示す。

• 尖度（Kurtosis）：確率分布の尖り具合を表す統計量である。正規分布からの逸脱度を示し、インパルス性の評価に用いる。

リズム・時間構造特徴量

リズム・時間構造特徴量は、音楽のビート検出やテンポ推定に用いる。

• オンセット：音の立ち上がり時刻である。音符の開始時刻や打楽器の打撃時刻に対応する。

• オンセット強度包絡線：各時刻でのオンセットの強さを表す時系列である。スペクトルの変化量から算出する。

• テンポグラム：時間-テンポ表現である。各時刻での周期性の強さをテンポ（BPM：Beats Per Minute、1分あたりの拍数）の関数として表示する。

環境音解析

環境音解析は、音声以外の音響信号を分析する分野である。都市騒音の監視や生態系のモニタリングに応用する。

• 環境音：音声以外の音響信号の総称である。自然音（雨、風、鳥の鳴き声）、人工音（機械音、交通音）、生活音などを含む。

• 非定常性：信号の統計的性質が時間とともに変化する性質である。環境音の多くは非定常な特性を持つ。

• ノイズの色：ノイズのスペクトル特性を色で表現する概念である。ホワイトノイズ（パワースペクトル密度が周波数に依存せず一定）、ピンクノイズ（パワースペクトル密度が周波数に反比例：1/f）、ブラウンノイズ（パワースペクトル密度が周波数の二乗に反比例：1/f²）などがある。

• スペクトル傾き：両対数軸（横軸：周波数の対数、縦軸：PSDの対数）でのパワースペクトル密度の傾きである。ノイズの色を定量化する指標であり、線形回帰で算出する。

信号処理技術

• 正規化：信号の振幅を一定の範囲（-1から1など）に収める処理である。クリッピング（振幅が表現可能範囲を超えた際に頭打ちとなる現象）を防ぐために行う。

• プリエンファシス：音声の高域成分を強調する前処理である。1次差分フィルタ（係数は0.97程度）で実現する。声帯波源に由来するスペクトルの傾斜（高域ほどパワーが小さい）を補正し、線形予測分析によるフォルマント推定の精度を向上させる。

• 減衰：信号の振幅が時間とともに減少する現象である。打撃音などでは指数減衰でモデル化する。

• 共鳴：特定の周波数で振幅が増幅される現象である。声道のフォルマントや楽器の音色形成に関係する。

評価・検証手法

• 閾値処理：信号や特徴量が閾値（判定の基準値）を超えるか否かで判定を行う処理である。有声/無声判定やオンセット検出で用いる。

• 探索範囲：推定するパラメータの値域として指定する範囲である。本ガイドの演習では、F0の探索範囲として50-2000Hzを指定する。会話音声を対象とする場合は50-500Hz程度に絞ると、オクターブ誤りを減らせる。

第3章 Pythonでの実現

本章では、用語リストの各項目に対応するPython関数を表形式で示す。scipy、librosa、numpyの3つのライブラリを使用し、グラフ表示にmatplotlibを使用する。インストールはpip install --no-user scipy librosa numpy matplotlibで行う。コード例は各演習の章を参照のこと。librosa.pyin、librosa.lpc のように * の後ろに定義された引数はキーワード引数として渡す必要がある（例：librosa.pyin(y, fmin=..., fmax=..., sr=...)）。グラフの日本語ラベルが文字化けする場合は、plt.rcParams['font.family'] = 'MS Gothic' のようにWindows標準の日本語フォントを指定する。

音声信号の基本概念

周波数領域分析

スペクトル特徴量

基本周波数推定

フォルマント分析

時間領域特徴量

リズム・時間構造特徴量

環境音解析

信号処理技術

評価・検証手法

第4章 演習1：音声ファイルの入出力

テーマ名：音声ファイルの読み込み・書き出しと基本情報の取得。

手順：scipyとlibrosaによるWAVファイル読み込みを比較し、サンプリング周波数、チャンネル数、サンプル数、時間長を表示する。テスト信号として440Hz（A4音）のサイン波を生成し、WAVファイルとして保存する。WAVファイル書き出し時は、ピーク振幅が1を超える場合に信号全体を縮小し、波形形状を保ったまま-1.0から1.0の範囲に収める。次のコードを実行する。

ヒント：librosa.loadはデフォルトで22050Hzにリサンプリングする。元のサンプリング周波数を保持するにはsr=Noneを指定する。また、librosaはデフォルトでモノラル化するため、チャンネル数を保持して読み込むにはmono=Falseを指定する。

考察ポイント：scipyとlibrosaでサンプリング周波数やサンプル数の表示が一致するかを確認し、生成したサイン波の波形が440Hzの周期に対応しているかを読み取る。

import numpy as np
from scipy.io import wavfile
import librosa
import matplotlib.pyplot as plt

def read_wav_scipy(filename):
    """scipyによるWAVファイル読み込み"""
    sampling_rate, data = wavfile.read(filename)

    channel_count = 1 if data.ndim == 1 else data.shape[1]
    sample_count = data.shape[0]

    print("=== Scipyでの読み込み結果 ===")
    print(f"サンプリング周波数: {sampling_rate} Hz")
    print(f"チャンネル数: {channel_count}")
    print(f"サンプル数: {sample_count}")
    print(f"時間長: {sample_count / sampling_rate:.2f} 秒")
    print(f"データ型: {data.dtype}")
    return sampling_rate, data

def read_wav_librosa(filename):
    """librosaによるWAVファイル読み込み"""
    data, sampling_rate = librosa.load(filename, sr=None, mono=False)

    channel_count = 1 if data.ndim == 1 else data.shape[0]
    sample_count = data.shape[-1]

    print("\n=== Librosaでの読み込み結果 ===")
    print(f"サンプリング周波数: {sampling_rate} Hz")
    print(f"チャンネル数: {channel_count}")
    print(f"サンプル数: {sample_count}")
    print(f"時間長: {sample_count / sampling_rate:.2f} 秒")
    print(f"データ型: {data.dtype}")
    return sampling_rate, data

def write_wav(filename, sampling_rate, data):
    """WAVファイルの書き出し"""
    peak = np.max(np.abs(data))
    normalized_data = data / max(1.0, peak)
    int_data = np.round(normalized_data * 32767).astype(np.int16)
    wavfile.write(filename, sampling_rate, int_data)

def generate_test_signal():
    """テスト用の音声信号生成（1秒のサイン波）"""
    duration = 1.0
    sampling_rate = 44100
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal_data = 0.5 * np.sin(2 * np.pi * 440 * t)
    return sampling_rate, signal_data

# テスト信号の生成と保存
sampling_rate, test_signal = generate_test_signal()
write_wav("test_signal.wav", sampling_rate, test_signal)

# 保存した音声ファイルの読み込みと情報表示
sr_scipy, data_scipy = read_wav_scipy("test_signal.wav")
sr_librosa, data_librosa = read_wav_librosa("test_signal.wav")

# 波形プロット
plt.figure(figsize=(12, 4))
time = np.arange(len(test_signal)) / sampling_rate
plt.plot(time, test_signal)
plt.xlabel('時間 [秒]')
plt.ylabel('振幅')
plt.title('テスト信号の波形')
plt.grid(True)
plt.show()

第5章 演習2：スペクトログラムの生成と表示

テーマ名：短時間フーリエ変換（STFT）によるスペクトログラム生成。スペクトログラムにより、周波数成分の時間変化を可視化できる。

手順：チャープ信号（周波数が100Hzから4000Hzへ指数関数的に変化する信号。scipy.signal.chirpではmethod='logarithmic'で指定する）を生成し、scipyとlibrosaでスペクトログラムを計算・比較する。窓関数（ハン窓）、フレーム長（2048サンプル）、ホップ長（1024サンプル）の設定例を確認する。次のコードを実行する。

ヒント：フレーム長を長くすると周波数分解能は向上するが時間分解能は低下する（トレードオフ）。scipyとlibrosaの表示は、比較しやすいよう最大値を0dBとする相対パワーで表示する。

考察ポイント：スペクトログラム上で周波数が時間とともに上昇する様子を読み取り、scipyとlibrosaの結果で周波数の時間変化の傾向が一致するかを確認する。

import numpy as np
from scipy import signal
import librosa
import librosa.display
import matplotlib.pyplot as plt

def generate_chirp_signal():
    """チャープ信号の生成（周波数が時間とともに指数関数的に変化する信号）"""
    duration = 3.0
    sampling_rate = 44100
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    chirp = signal.chirp(t, f0=100, f1=4000, t1=duration, method='logarithmic')
    return sampling_rate, chirp

def calculate_spectrogram_scipy(audio_data, sampling_rate):
    """Scipyを使用したスペクトログラム計算"""
    nperseg = 2048
    noverlap = nperseg // 2
    frequencies, times, Sxx = signal.spectrogram(
        audio_data,
        fs=sampling_rate,
        window='hann',
        nperseg=nperseg,
        noverlap=noverlap,
        scaling='spectrum'
    )
    Sxx_db = 10 * np.log10(Sxx / (np.max(Sxx) + 1e-12) + 1e-12)
    return frequencies, times, Sxx_db

def calculate_spectrogram_librosa(audio_data, sampling_rate):
    """Librosaを使用したスペクトログラム計算"""
    n_fft = 2048
    hop_length = n_fft // 2
    D = librosa.stft(
        audio_data,
        n_fft=n_fft,
        hop_length=hop_length,
        window='hann',
        center=True
    )
    S_db = librosa.power_to_db(np.abs(D) ** 2, ref=np.max)
    return S_db

# チャープ信号の生成
sampling_rate, chirp = generate_chirp_signal()

# Scipyによるスペクトログラム
freqs, times, spec_scipy = calculate_spectrogram_scipy(chirp, sampling_rate)

# Librosaによるスペクトログラム
spec_librosa = calculate_spectrogram_librosa(chirp, sampling_rate)

# 結果の可視化
plt.figure(figsize=(15, 10))

# 波形
plt.subplot(3, 1, 1)
time = np.arange(len(chirp)) / sampling_rate
plt.plot(time, chirp)
plt.xlabel('時間 [秒]')
plt.ylabel('振幅')
plt.title('チャープ信号の波形')
plt.grid(True)

# Scipyによるスペクトログラム
plt.subplot(3, 1, 2)
plt.pcolormesh(times, freqs, spec_scipy, shading='auto')
plt.ylabel('周波数 [Hz]')
plt.xlabel('時間 [秒]')
plt.title('Scipyによるスペクトログラム')
plt.colorbar(label='相対パワー [dB]')

# Librosaによるスペクトログラム
plt.subplot(3, 1, 3)
librosa.display.specshow(
    spec_librosa,
    sr=sampling_rate,
    x_axis='time',
    y_axis='hz',
    hop_length=1024,
    n_fft=2048
)
plt.ylabel('周波数 [Hz]')
plt.xlabel('時間 [秒]')
plt.title('Librosaによるスペクトログラム')
plt.colorbar(label='相対パワー [dB]')

plt.tight_layout()
plt.show()

第6章 演習3：基本周波数（F0）推定

テーマ名：音声の高さを決定する基本周波数を推定する3つの手法の比較。基本周波数の推定は、音声の抑揚分析や話者識別の基礎となる。

手順：440Hzの基本波とその倍音を含むテスト信号に対し、自己相関法、ケプストラム法、PYIN法を適用する。各手法の推定結果を真のF0と比較する。次のコードを実行する。

ヒント：自己相関法では、ラグ0近傍の値は信号の滑らかさを強く反映するため、そのまま最大値を選ぶと小さいラグを誤って選ぶことがある。本プログラムでは、F0の探索範囲（50-2000Hz）に対応するラグ範囲を計算し、その範囲内の局所ピークを用いる。実際の音声では無声区間でF0が検出されないことがあり、voiced_flag（有声/無声を示すブール配列）で判定する。

考察ポイント：3手法の推定F0が真のF0（440Hz）にどの程度一致するかを比較し、手法ごとの推定のばらつきやオクターブ誤りの有無を読み取る。

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

def generate_test_signal(f0=440, duration=1.0, sr=44100):
    """基本周波数を持つテスト信号の生成"""
    t = np.linspace(0, duration, int(sr * duration), endpoint=False)

    # 基本波と倍音を含む信号
    signal_data = (
        0.5 * np.sin(2 * np.pi * f0 * t)
        + 0.25 * np.sin(2 * np.pi * 2 * f0 * t)
        + 0.125 * np.sin(2 * np.pi * 3 * f0 * t)
    )
    return signal_data, sr

def autocorrelation_f0(signal_data, sr, frame_length=2048, hop_length=512,
                       fmin=50, fmax=2000):
    """自己相関法によるF0推定"""
    f0_values = []
    times = []

    min_lag = max(1, int(np.floor(sr / fmax)))
    max_lag = min(frame_length - 1, int(np.ceil(sr / fmin)))
    window = signal.windows.hann(frame_length, sym=False)

    for i in range(0, len(signal_data) - frame_length + 1, hop_length):
        frame = signal_data[i:i + frame_length]
        frame = frame - np.mean(frame)
        frame = frame * window

        # 自己相関の計算
        corr = signal.correlate(frame, frame, mode='full', method='fft')
        corr = corr[len(corr) // 2:]

        # ラグ0の値で正規化
        corr = corr / (corr[0] + 1e-10)

        # F0の探索範囲に対応するラグ範囲で局所ピークを探索
        search = corr[min_lag:max_lag + 1]
        peaks, _ = signal.find_peaks(search, height=0.1)

        if len(peaks) > 0:
            # ラグ0に最も近い有意な局所ピークを基本周期とする
            peak_idx = peaks[0] + min_lag
            f0 = sr / peak_idx
            if fmin <= f0 <= fmax:
                f0_values.append(f0)
            else:
                f0_values.append(np.nan)
        else:
            f0_values.append(np.nan)

        times.append((i + frame_length / 2) / sr)

    return np.array(f0_values), np.array(times)

def cepstrum_f0(signal_data, sr, frame_length=2048, hop_length=512,
                fmin=50, fmax=2000):
    """ケプストラム法によるF0推定"""
    f0_values = []
    times = []

    qmin = max(1, int(np.floor(sr / fmax)))
    qmax = min(frame_length // 2, int(np.ceil(sr / fmin)))
    window = signal.windows.hann(frame_length, sym=False)

    for i in range(0, len(signal_data) - frame_length + 1, hop_length):
        frame = signal_data[i:i + frame_length]
        frame = frame - np.mean(frame)
        frame = frame * window

        # 対数パワースペクトルの計算
        spectrum = np.fft.rfft(frame)
        log_power_spectrum = np.log(np.abs(spectrum) ** 2 + 1e-10)

        # 実ケプストラムの計算
        cepstrum = np.fft.irfft(log_power_spectrum, n=frame_length)

        # 探索範囲でピーク検出
        search = cepstrum[qmin:qmax + 1]
        peaks, _ = signal.find_peaks(search)

        if len(peaks) > 0:
            peak_idx = peaks[np.argmax(search[peaks])] + qmin
            f0 = sr / peak_idx
            if fmin <= f0 <= fmax:
                f0_values.append(f0)
            else:
                f0_values.append(np.nan)
        else:
            f0_values.append(np.nan)

        times.append((i + frame_length / 2) / sr)

    return np.array(f0_values), np.array(times)

def librosa_f0(signal_data, sr, frame_length=2048, hop_length=512,
               fmin=50, fmax=2000):
    """librosaを使用したF0推定（PYIN法）"""
    f0, voiced_flag, voiced_probs = librosa.pyin(
        signal_data,
        fmin=fmin,
        fmax=fmax,
        sr=sr,
        frame_length=frame_length,
        hop_length=hop_length
    )
    return f0, voiced_flag, frame_length, hop_length

# テスト信号の生成と分析
signal_data, sr = generate_test_signal(f0=440)

# 自己相関法によるF0推定
f0_autocorr, times_autocorr = autocorrelation_f0(signal_data, sr)

# ケプストラム法によるF0推定
f0_cepstrum, times_cepstrum = cepstrum_f0(signal_data, sr)

# librosaによるF0推定
f0_librosa, voiced_flag, frame_length, hop_length = librosa_f0(signal_data, sr)
times_librosa = librosa.times_like(
    f0_librosa,
    sr=sr,
    hop_length=hop_length,
    n_fft=frame_length
)

# 結果の可視化
plt.figure(figsize=(12, 10))

# 波形
plt.subplot(4, 1, 1)
t = np.arange(len(signal_data)) / sr
plt.plot(t, signal_data)
plt.xlabel('時間 [秒]')
plt.ylabel('振幅')
plt.title('テスト信号の波形')
plt.grid(True)

# 自己相関法によるF0
plt.subplot(4, 1, 2)
plt.plot(times_autocorr, f0_autocorr, 'o-', label='推定F0')
plt.axhline(y=440, color='r', linestyle='--', label='真のF0')
plt.xlabel('時間 [秒]')
plt.ylabel('周波数 [Hz]')
plt.title('自己相関法によるF0推定結果')
plt.legend()
plt.grid(True)

# ケプストラム法によるF0
plt.subplot(4, 1, 3)
plt.plot(times_cepstrum, f0_cepstrum, 'o-', label='推定F0')
plt.axhline(y=440, color='r', linestyle='--', label='真のF0')
plt.xlabel('時間 [秒]')
plt.ylabel('周波数 [Hz]')
plt.title('ケプストラム法によるF0推定結果')
plt.legend()
plt.grid(True)

# librosaによるF0
plt.subplot(4, 1, 4)
plt.plot(
    times_librosa[voiced_flag],
    f0_librosa[voiced_flag],
    'o-',
    label='推定F0（有声区間）'
)
plt.axhline(y=440, color='r', linestyle='--', label='真のF0')
plt.xlabel('時間 [秒]')
plt.ylabel('周波数 [Hz]')
plt.title('PYIN（librosa）によるF0推定結果')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

第7章 演習4：フォルマント分析

テーマ名：母音の音色を特徴づけるフォルマントの線形予測分析（LPC）による推定。フォルマントの違いにより、「あ」「い」「う」などの母音を区別できる。

手順：母音/a/と/i/を合成し、プリエンファシス後にLPCでフォルマントを推定する。母音合成では、各フォルマントを表す2次共鳴フィルタを直列に接続し、声道の全極フィルタを模擬する。LPC係数の根のうち、実根、負の虚部を持つ根、バンド幅の広い根はフォルマントに対応しないため除外する。有声音と無声音の違いをスペクトルで比較する。次のコードを実行する。

ヒント：実音声での予測次数の目安はサンプリング周波数（kHz単位）に2〜4を加えた程度（16kHzなら18前後）である。本演習の合成母音はフォルマントが3つ（極が6個）であるため、次数12で推定できる。次数が小さいとフォルマントを正確に捉えられず、大きすぎるとノイズの影響を受ける。バンド幅は、極の半径をrとして -(sr / π) log(r) で計算する。

考察ポイント：推定フォルマントが真のフォルマント（F1、F2、F3）にどの程度一致するかを読み取り、母音/a/と/i/でF1・F2の位置がどう異なるか、無声音のスペクトルが有声音と比べてどのような形状かを確認する。

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

def generate_vowel(formants, duration=1.0, sr=16000):
    """フォルマントを持つ母音の生成"""
    t = np.linspace(0, duration, int(sr * duration), endpoint=False)
    f0 = 130  # 基本周波数

    # 声帯波源（のこぎり波）
    source = signal.sawtooth(2 * np.pi * f0 * t)

    # フォルマントフィルタを直列に適用して声道を模擬
    vowel = source.copy()
    for formant in formants:
        bw = formant * 0.1  # バンド幅
        w0 = 2 * np.pi * formant / sr
        r = np.exp(-np.pi * bw / sr)

        # 共鳴フィルタ（2次IIRフィルタ）
        a = [1, -2 * r * np.cos(w0), r ** 2]
        b = [1 - r]
        vowel = signal.lfilter(b, a, vowel)

    return vowel / (np.max(np.abs(vowel)) + 1e-10), sr

def generate_unvoiced(duration=1.0, sr=16000, seed=1):
    """無声音の生成"""
    rng = np.random.default_rng(seed)

    # ランダムノイズを生成
    noise = rng.standard_normal(int(sr * duration))

    # 高域通過フィルタで無声子音を模擬
    sos = signal.butter(4, 2000, 'hp', fs=sr, output='sos')
    unvoiced = signal.sosfilt(sos, noise)

    return unvoiced / (np.max(np.abs(unvoiced)) + 1e-10), sr

def estimate_formants(signal_data, sr, order=12):
    """LPCによるフォルマント推定"""
    # プリエンファシス（声帯波源によるスペクトル傾斜を補正）
    emphasized = signal.lfilter([1, -0.97], [1], signal_data)

    # 中央の短時間フレームを用いて定常区間を分析
    frame_length = int(0.05 * sr)  # 50 ms
    start = (len(emphasized) - frame_length) // 2
    frame = emphasized[start:start + frame_length]
    frame = frame - np.mean(frame)
    frame = frame * signal.windows.hamming(frame_length, sym=False)

    # LPC係数の計算
    a = librosa.lpc(y=frame, order=order)

    # 根を求める
    roots = np.roots(a)
    roots = roots[np.imag(roots) > 0]  # 正の虚部を持つ根のみ

    # 角周波数から周波数へ変換
    freqs = np.angle(roots) * (sr / (2 * np.pi))

    # バンド幅の計算
    bandwidth = -(sr / np.pi) * np.log(np.abs(roots))

    # フォルマント候補の選別
    keep = (
        (freqs > 50)
        & (freqs < sr / 2)
        & (bandwidth > 0)
        & (bandwidth < 400)
    )
    freqs = freqs[keep]
    bandwidth = bandwidth[keep]

    # 周波数でソート
    sorted_idx = np.argsort(freqs)
    freqs = freqs[sorted_idx]
    bandwidth = bandwidth[sorted_idx]

    return freqs, bandwidth

# 母音/a/のフォルマント
formants_a = [730, 1090, 2440]
# 母音/i/のフォルマント
formants_i = [270, 2290, 3010]

# 母音の生成
vowel_a, sr = generate_vowel(formants_a)
vowel_i, _ = generate_vowel(formants_i)

# 無声音の生成
unvoiced, _ = generate_unvoiced()

# フォルマント推定
freqs_a, bw_a = estimate_formants(vowel_a, sr)
freqs_i, bw_i = estimate_formants(vowel_i, sr)

print("母音/a/ 推定フォルマント [Hz]:", np.round(freqs_a[:3], 1))
print("母音/i/ 推定フォルマント [Hz]:", np.round(freqs_i[:3], 1))

# 結果の可視化
plt.figure(figsize=(15, 12))

# 母音/a/の波形とスペクトル
plt.subplot(3, 2, 1)
plt.plot(vowel_a)
plt.title('母音/a/の波形（有声音）')
plt.xlabel('サンプル')
plt.ylabel('振幅')

plt.subplot(3, 2, 2)
freq, spec = signal.welch(vowel_a, sr, nperseg=2048)
plt.plot(freq, 10 * np.log10(spec + 1e-12))
ymin, ymax = plt.ylim()
plt.vlines(freqs_a[:3], ymin, ymax,
           colors='r', linestyles='--', label='推定フォルマント')
plt.vlines(formants_a, ymin, ymax,
           colors='g', linestyles=':', label='真のフォルマント')
plt.title('母音/a/のスペクトル')
plt.xlabel('周波数 [Hz]')
plt.ylabel('パワー [dB]')
plt.legend()

# 母音/i/の波形とスペクトル
plt.subplot(3, 2, 3)
plt.plot(vowel_i)
plt.title('母音/i/の波形（有声音）')
plt.xlabel('サンプル')
plt.ylabel('振幅')

plt.subplot(3, 2, 4)
freq, spec = signal.welch(vowel_i, sr, nperseg=2048)
plt.plot(freq, 10 * np.log10(spec + 1e-12))
ymin, ymax = plt.ylim()
plt.vlines(freqs_i[:3], ymin, ymax,
           colors='r', linestyles='--', label='推定フォルマント')
plt.vlines(formants_i, ymin, ymax,
           colors='g', linestyles=':', label='真のフォルマント')
plt.title('母音/i/のスペクトル')
plt.xlabel('周波数 [Hz]')
plt.ylabel('パワー [dB]')
plt.legend()

# 無声音の波形とスペクトル
plt.subplot(3, 2, 5)
plt.plot(unvoiced)
plt.title('無声音の波形')
plt.xlabel('サンプル')
plt.ylabel('振幅')

plt.subplot(3, 2, 6)
freq, spec = signal.welch(unvoiced, sr, nperseg=2048)
plt.plot(freq, 10 * np.log10(spec + 1e-12))
plt.title('無声音のスペクトル')
plt.xlabel('周波数 [Hz]')
plt.ylabel('パワー [dB]')

plt.tight_layout()
plt.show()

第8章 演習5：環境音の解析

テーマ名：非定常な環境音の合成と、スペクトル特徴量・時間領域特徴量・リズム特徴量による解析。環境音解析は、騒音監視や異常音検知に応用できる。

手順：雨音を大中小の雨滴（ポアソン分布に従う発生数）でシミュレートし、スペクトル重心、RMSエネルギー、ゼロ交差率、テンポグラム、ノイズの色（スペクトル傾き）を計算・可視化する。乱数を用いるため、結果を再現できるよう乱数シードを固定する。次のコードを実行する。

ヒント：環境音は非定常であるため、短時間分析で時間変化を捉える。スペクトル傾きが0に近い場合はホワイトノイズ、約-1の場合はピンクノイズ、約-2の場合はブラウンノイズに近い特性を示す。本演習の雨音は、雨滴ごとの指数減衰と低域通過型の1次IIRフィルタの影響で、負のスペクトル傾きを持つ音として観察される。

考察ポイント：クレストファクターが大きいほどインパルス性が高い。尖度が正規分布（値0）より大きければ突発的な音が含まれる。これらの値から、合成した雨音がどのノイズの色に近く、どの程度インパルス性を持つかを読み取る。

import numpy as np
import librosa
import librosa.display
import scipy.stats
from scipy.signal import welch
from scipy import signal
import matplotlib.pyplot as plt

def compute_spectral_features(y, sr, n_fft=2048, hop_length=512):
    """スペクトル特徴量の計算"""
    spectral_centroid = librosa.feature.spectral_centroid(
        y=y, sr=sr, n_fft=n_fft, hop_length=hop_length
    )[0]

    spectral_bandwidth = librosa.feature.spectral_bandwidth(
        y=y, sr=sr, n_fft=n_fft, hop_length=hop_length
    )[0]

    spectral_rolloff = librosa.feature.spectral_rolloff(
        y=y, sr=sr, n_fft=n_fft, hop_length=hop_length
    )[0]

    spectral_flatness = librosa.feature.spectral_flatness(
        y=y, n_fft=n_fft, hop_length=hop_length
    )[0]

    return {
        'centroid': spectral_centroid,
        'bandwidth': spectral_bandwidth,
        'rolloff': spectral_rolloff,
        'flatness': spectral_flatness
    }

def compute_temporal_features(y, sr, frame_length=2048, hop_length=512):
    """時間領域特徴量の計算"""
    # RMSエネルギー
    rms = librosa.feature.rms(
        y=y, frame_length=frame_length, hop_length=hop_length
    )[0]

    # ゼロ交差率
    zero_crossing_rate = librosa.feature.zero_crossing_rate(
        y=y, frame_length=frame_length, hop_length=hop_length
    )[0]

    # 波形の統計量
    rms_val = np.sqrt(np.mean(y ** 2)) + 1e-10
    crest_factor = np.max(np.abs(y)) / rms_val
    kurtosis = scipy.stats.kurtosis(y, fisher=True)

    return {
        'rms': rms,
        'zcr': zero_crossing_rate,
        'crest_factor': crest_factor,
        'kurtosis': kurtosis
    }

def compute_rhythm_features(y, sr, hop_length=512):
    """リズム特徴量の計算"""
    # オンセット強度包絡線
    onset_env = librosa.onset.onset_strength(y=y, sr=sr, hop_length=hop_length)

    # テンポグラム
    tempogram = librosa.feature.tempogram(
        onset_envelope=onset_env,
        sr=sr,
        hop_length=hop_length
    )

    # オンセット検出
    onset_frames = librosa.onset.onset_detect(
        onset_envelope=onset_env,
        sr=sr,
        hop_length=hop_length
    )

    return {
        'tempogram': tempogram,
        'onset_frames': onset_frames
    }

def analyze_noise_color(y, sr, nperseg=4096, f_low=100, f_high=10000):
    """ノイズの色分析（パワースペクトル密度の傾き）"""
    freqs, psd = welch(y, sr, nperseg=nperseg)

    # 低すぎる周波数とナイキスト周波数付近を避け、両対数軸で傾きを求める
    band = (freqs >= f_low) & (freqs <= min(f_high, sr * 0.45))
    log_freqs = np.log10(freqs[band])
    log_psd = np.log10(psd[band] + 1e-20)

    # 線形回帰で傾きを計算
    slope, intercept, _, _, _ = scipy.stats.linregress(log_freqs, log_psd)
    fitted_psd = 10 ** (intercept + slope * log_freqs)

    return slope, freqs[1:], psd[1:], freqs[band], fitted_psd

def simulate_rain(duration=5.0, sr=44100, seed=0):
    """雨音のシミュレーション"""
    rng = np.random.default_rng(seed)
    rain = np.zeros(int(sr * duration))

    # 雨滴の種類（大中小）ごとのパラメータ
    drop_types = {
        'small': {
            'freq_range': (2000, 4000),
            'duration_range': (0.01, 0.03),
            'amplitude_range': (0.1, 0.3),
            'rate': 100  # 1秒あたりの平均発生数
        },
        'medium': {
            'freq_range': (1000, 2000),
            'duration_range': (0.03, 0.06),
            'amplitude_range': (0.3, 0.6),
            'rate': 50
        },
        'large': {
            'freq_range': (500, 1000),
            'duration_range': (0.06, 0.1),
            'amplitude_range': (0.6, 1.0),
            'rate': 20
        }
    }

    # 各種類の雨滴を生成
    for drop_type, params in drop_types.items():
        # 雨滴の数（ポアソン分布に従う）
        n_drops = rng.poisson(duration * params['rate'])

        for _ in range(n_drops):
            # 雨滴の発生時刻
            drop_time = rng.random() * duration
            idx = int(drop_time * sr)

            # 雨滴のパラメータ
            freq = rng.uniform(*params['freq_range'])
            drop_duration = rng.uniform(*params['duration_range'])
            amplitude = rng.uniform(*params['amplitude_range'])

            # 雨滴の波形生成
            n_samples = int(drop_duration * sr)
            t_drop = np.arange(n_samples) / sr
            envelope = np.exp(-t_drop * 30)  # 指数減衰
            drop = amplitude * envelope * np.sin(2 * np.pi * freq * t_drop)

            # 1次IIRフィルタ（低域通過）で音色を変化
            drop = signal.lfilter([0.01], [1, -0.99], drop)

            # 雨音に加算
            end = min(idx + n_samples, len(rain))
            rain[idx:end] += drop[:end - idx]

    # 振幅の正規化
    rain = rain / (np.max(np.abs(rain)) + 1e-10)

    return rain, sr

# 解析の実行
y, sr = simulate_rain()

# パラメータ
n_fft = 2048
hop_length = 512

# 各種特徴量の計算
spectral_features = compute_spectral_features(y, sr, n_fft=n_fft, hop_length=hop_length)
temporal_features = compute_temporal_features(y, sr, frame_length=n_fft, hop_length=hop_length)
rhythm_features = compute_rhythm_features(y, sr, hop_length=hop_length)
noise_slope, freqs_psd, psd, freqs_fit, psd_fit = analyze_noise_color(y, sr)

# 結果の可視化
plt.figure(figsize=(15, 10))

# 波形
plt.subplot(3, 2, 1)
plt.plot(np.arange(len(y)) / sr, y)
plt.title('雨音波形')
plt.xlabel('時間 [秒]')
plt.ylabel('振幅')

# スペクトログラム
plt.subplot(3, 2, 2)
D = librosa.power_to_db(
    np.abs(librosa.stft(y, n_fft=n_fft, hop_length=hop_length)) ** 2,
    ref=np.max
)
librosa.display.specshow(
    D,
    sr=sr,
    hop_length=hop_length,
    n_fft=n_fft,
    x_axis='time',
    y_axis='log'
)
plt.colorbar(format='%+2.0f dB')
plt.title('スペクトログラム')

# スペクトル特徴量の時間変化
plt.subplot(3, 2, 3)
times = librosa.times_like(
    spectral_features['centroid'],
    sr=sr,
    hop_length=hop_length,
    n_fft=n_fft
)
plt.plot(times, spectral_features['flatness'], label='平坦度')
plt.plot(times, spectral_features['centroid'] / (sr / 2), label='重心（ナイキスト正規化）')
plt.legend()
plt.title('スペクトル特徴量の時間変化')
plt.xlabel('時間 [秒]')

# パワースペクトル密度（ノイズの色分析）
plt.subplot(3, 2, 4)
plt.loglog(freqs_psd, psd, label='PSD')
plt.loglog(freqs_fit, psd_fit, 'r--', label='回帰直線')
plt.title(f'パワースペクトル密度（傾き: {noise_slope:.2f}）')
plt.xlabel('周波数 [Hz]')
plt.ylabel('PSD')
plt.legend()

# テンポグラム
plt.subplot(3, 2, 5)
librosa.display.specshow(
    rhythm_features['tempogram'],
    sr=sr,
    hop_length=hop_length,
    x_axis='time',
    y_axis='tempo'
)
plt.title('テンポグラム')
plt.xlabel('時間 [秒]')
plt.ylabel('テンポ [BPM]')

# エネルギーとゼロ交差率
plt.subplot(3, 2, 6)
times_rms = librosa.times_like(
    temporal_features['rms'],
    sr=sr,
    hop_length=hop_length,
    n_fft=n_fft
)
plt.plot(times_rms, temporal_features['rms'], label='RMS')
plt.plot(times_rms, temporal_features['zcr'], label='ZCR')
plt.legend()
plt.title('時間領域特徴量')
plt.xlabel('時間 [秒]')

plt.tight_layout()
plt.show()

# 統計値の表示
print(f"クレストファクター: {temporal_features['crest_factor']:.2f}")
print(f"尖度: {temporal_features['kurtosis']:.2f}")
print(f"ノイズの色（スペクトル傾き）: {noise_slope:.2f}")