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

音声信号処理の基本

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

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

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

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

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

• 周波数分解能・時間分解能：信号分析における精度のトレードオフである。周波数分解能は周波数の識別精度、時間分解能は時間変化の追跡精度を示す。一方を高めると他方が低下する。

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

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

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

• ノイズ：目的とする信号に混入する不要な成分である。ランダムな性質を持つことが多く、信号対雑音比（SNR）で信号品質を評価する。

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

用語リスト

本節では、音声信号処理の専門用語を分野別に整理する。用語は以下の9分野に分類している。音声信号の基本概念、周波数領域分析、スペクトル特徴量、基本周波数推定、フォルマント分析、時間領域特徴量、リズム・時間構造特徴量、環境音解析、信号処理技術である。各用語に対応するPython関数は「Pythonでの実現」節の表で確認でき、コード例は「Pythonプログラム例」節に示す。

音声信号の基本概念

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

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

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

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

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

周波数領域分析

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

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

• 窓関数：信号を短い区間に分割する際に適用する関数である。ハン窓、ハミング窓、ブラックマン窓などがある。窓関数の選択は周波数分解能と時間分解能のトレードオフに影響する。

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

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

• オーバーラップ率：連続するフレーム間で重複する部分の割合である。通常50%や75%を使用する。オーバーラップにより時間分解能を向上できる。

• FFTサイズ：高速フーリエ変換（FFT）で使用するサンプル数である。2のべき乗（512、1024、2048など）が一般的である。大きいほど周波数分解能が高くなる。

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

• スペクトル包絡：スペクトルの全体的な形状を表す曲線である。声道の共鳴特性を反映し、音色の特徴を表す。

スペクトル特徴量

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

• スペクトル重心：パワースペクトルの重心周波数である。音色の明るさを表し、値が高いほど高周波成分が多い。

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

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

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

基本周波数推定

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

• 基本周波数（F0）：音声の周期的な振動の基本となる周波数である。人間の声では通常、男性で80-200Hz、女性で150-350Hz程度である。

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

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

• 自己相関法：信号と時間シフトした自身との相関を計算し、周期性を検出する手法である。基本周波数推定の古典的手法である。

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

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

• オクターブ誤り：基本周波数推定で、真の基本周波数の整数倍または整数分の1の値を誤って推定する現象である。

フォルマント分析

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

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

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

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

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

• 予測次数：線形予測分析で使用する過去のサンプル数である。通常、サンプリング周波数（kHz単位）の2倍程度に設定する。

• バンド幅：フォルマントの周波数幅である。共鳴のシャープさを表し、減衰の速さに関係する。

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

時間領域特徴量

時間領域特徴量は、信号の波形から直接計算する指標である。周波数分析を必要としないため計算が軽量であり、リアルタイム処理に適する。

• RMSエネルギー：信号の二乗平均平方根である。音声の強度を表す基本的な指標である。

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

• クレストファクター：信号のピーク値とRMS値の比である。波高率を表し、インパルス性の強さを示す。

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

リズム・時間構造特徴量

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

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

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

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

環境音解析

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

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

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

• ノイズの色：ノイズのスペクトル特性を色で表現する概念である。光のスペクトルとの類似性から命名された。ホワイトノイズ（全周波数で均等なパワー）、ピンクノイズ（周波数に反比例するパワー）、ブラウンノイズ（周波数の二乗に反比例するパワー）などがある。

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

信号処理技術

• 正規化：信号の振幅を一定の範囲（通常-1から1）に収める処理である。クリッピング（波形の頭切れ）を防ぐために行う。

• 減衰：信号の振幅が時間とともに減少する現象である。指数減衰が自然音の特徴として広く見られる。

• 共鳴：特定の周波数で振幅が増幅される現象である。声道のフォルマントや楽器の音色形成で重要な役割を果たす。

評価・検証手法

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

• 有効範囲：パラメータが妥当とされる値の範囲である。基本周波数推定では通常50-2000Hzが人間の音声の有効範囲である。

Pythonでの実現

本節では、用語リストの各項目に対応するPython関数を表形式で示す。主にscipy、librosa、numpyの3つのライブラリを使用する。インストールはpip install scipy librosa numpyで行う。コード例は「Pythonプログラム例」節を参照のこと。

音声信号の基本概念

周波数領域分析

スペクトル特徴量

基本周波数推定

フォルマント分析

時間領域特徴量

リズム・時間構造特徴量

環境音解析

信号処理技術

評価・検証手法

Pythonプログラム例

本節では、音声信号処理のプログラム例を5つ示す。各プログラムは独立して実行可能であり、用語リストの概念を実際のコードで確認できる。実行にはPython 3.8以降が必要である。グラフ表示にはmatplotlibを使用するため、pip install matplotlibでインストールする。

1. 音声ファイルの入出力

目的：音声ファイルの読み込み・書き出しと基本情報の取得方法を理解する。

処理内容：scipyとlibrosaによるWAVファイル読み込みを比較し、サンプリング周波数、ビット深度、チャンネル数、時間長を表示する。テスト信号として440Hz（A4音）のサイン波を生成し、WAVファイルとして保存する。

注意事項：WAVファイル書き出し時は正規化処理で振幅を-1.0から1.0に制限し、クリッピングを防ぐ。librosa.loadはデフォルトで22050Hzにリサンプリングするため、元のサンプリング周波数を保持するにはsr=Noneを指定する。

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)
    print("=== Scipyでの読み込み結果 ===")
    print(f"サンプリング周波数: {sampling_rate} Hz")
    print(f"チャンネル数: {data.shape[1] if len(data.shape) > 1 else 1}")
    print(f"サンプル数: {len(data)}")
    print(f"時間長: {len(data)/sampling_rate:.2f} 秒")
    return sampling_rate, data

def read_wav_librosa(filename):
    """librosaによるWAVファイル読み込み"""
    data, sampling_rate = librosa.load(filename, sr=None)
    print("\n=== Librosaでの読み込み結果 ===")
    print(f"サンプリング周波数: {sampling_rate} Hz")
    print(f"チャンネル数: {data.shape[1] if len(data.shape) > 1 else 1}")
    print(f"サンプル数: {len(data)}")
    print(f"時間長: {len(data)/sampling_rate:.2f} 秒")
    return sampling_rate, data

def write_wav(filename, sampling_rate, data):
    """WAVファイルの書き出し"""
    normalized_data = np.clip(data, -1.0, 1.0)
    int_data = (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))
    signal = 0.5 * np.sin(2 * np.pi * 440 * t)
    return sampling_rate, signal

# テスト信号の生成と保存
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()

2. スペクトログラムの生成と表示

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

処理内容：チャープ信号（周波数が100Hzから4000Hzへ指数的に変化する信号）を生成し、scipyとlibrosaでスペクトログラムを計算・比較する。窓関数（ハン窓）、フレーム長（2048サンプル）、ホップ長（1024サンプル）の設定例を示す。

注意事項：フレーム長を長くすると周波数分解能は向上するが時間分解能は低下する（トレードオフ）。用途に応じてパラメータを調整する。

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))
    chirp = signal.chirp(t, f0=100, f1=4000, t1=duration, method='exponential')
    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'
    )
    return frequencies, times, 10 * np.log10(Sxx + 1e-10)

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.amplitude_to_db(np.abs(D), 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='gouraud')
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)
plt.ylabel('周波数 [Hz]')
plt.xlabel('時間 [秒]')
plt.title('Librosaによるスペクトログラム')
plt.colorbar(label='振幅 [dB]')

plt.tight_layout()
plt.show()

3. 基本周波数（F0）推定

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

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

注意事項：有効範囲（50-2000Hz）のチェックによりオクターブ誤りを軽減できる。実際の音声では無声区間でF0が検出されないことがあり、voiced_flagで判定する。

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))
    # 基本波と倍音を含む信号
    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):
    """自己相関法によるF0推定"""
    f0_values = []
    times = []

    for i in range(0, len(signal_data) - frame_length, hop_length):
        frame = signal_data[i:i + frame_length]
        # 自己相関の計算
        corr = signal.correlate(frame, frame, mode='full')
        corr = corr[len(corr)//2:]

        # 最初のピークを探す（直流成分を除く）
        peak_idx = np.argmax(corr[50:]) + 50
        if corr[peak_idx] > 0.1:  # 有声音の閾値
            f0 = sr / peak_idx
            if 50 <= f0 <= 2000:  # 有効なF0範囲
                f0_values.append(f0)
            else:
                f0_values.append(0)
        else:
            f0_values.append(0)

        times.append(i / sr)

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

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

    for i in range(0, len(signal_data) - frame_length, hop_length):
        frame = signal_data[i:i + frame_length]

        # パワースペクトルの計算
        spectrum = np.fft.fft(frame)
        log_spectrum = np.log(np.abs(spectrum) + 1e-10)

        # ケプストラムの計算
        cepstrum = np.fft.ifft(log_spectrum).real

        # 有効なF0範囲でピーク検出
        start_idx = int(sr / 2000)  # 2000Hz対応
        end_idx = int(sr / 50)      # 50Hz対応
        peak_idx = np.argmax(cepstrum[start_idx:end_idx]) + start_idx

        if peak_idx > 0:
            f0 = sr / peak_idx
            if 50 <= f0 <= 2000:
                f0_values.append(f0)
            else:
                f0_values.append(0)
        else:
            f0_values.append(0)

        times.append(i / sr)

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

def librosa_f0(signal_data, sr):
    """librosaを使用したF0推定（PYIN法）"""
    f0, voiced_flag, voiced_probs = librosa.pyin(
        signal_data,
        fmin=librosa.note_to_hz('C2'),
        fmax=librosa.note_to_hz('C7'),
        sr=sr
    )
    return f0, voiced_flag

# テスト信号の生成と分析
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 = librosa_f0(signal_data, sr)
times_librosa = librosa.times_like(f0_librosa, sr=sr)

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

# 波形
plt.subplot(4, 1, 1)
t = np.linspace(0, len(signal_data)/sr, len(signal_data))
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()

4. フォルマント分析

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

処理内容：母音/a/と/i/を合成し、LPCでフォルマントを推定する。有声音と無声音の違いをスペクトルで比較する。

注意事項：予測次数は通常サンプリング周波数（kHz）の2倍程度（16kHzなら次数12程度）に設定する。次数が小さいとフォルマントを正確に捉えられず、大きすぎるとノイズの影響を受ける。

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))
    f0 = 130  # 基本周波数

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

    # フォルマントフィルタの作成
    vowel = np.zeros_like(t)
    for formant in formants:
        bw = formant * 0.1  # バンド幅
        w0 = 2 * np.pi * formant / sr
        r = np.exp(-np.pi * bw / sr)

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

    return vowel / np.max(np.abs(vowel)), sr

def generate_unvoiced(duration=1.0, sr=16000):
    """無声音の生成"""
    # ランダムノイズを生成
    noise = np.random.randn(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)), sr

def estimate_formants(signal_data, sr, order=12):
    """LPCによるフォルマント推定"""
    # LPC係数の計算
    a = librosa.lpc(signal_data, order=order)

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

    # 角周波数から周波数へ変換
    angles = np.arctan2(np.imag(roots), np.real(roots))
    freqs = angles * (sr / (2 * np.pi))

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

    # 周波数でソート
    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)

# 結果の可視化
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))
plt.vlines(freqs_a[:3], plt.ylim()[0], plt.ylim()[1],
          colors='r', linestyles='--', label='推定フォルマント')
plt.vlines(formants_a, plt.ylim()[0], plt.ylim()[1],
          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))
plt.vlines(freqs_i[:3], plt.ylim()[0], plt.ylim()[1],
          colors='r', linestyles='--', label='推定フォルマント')
plt.vlines(formants_i, plt.ylim()[0], plt.ylim()[1],
          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))
plt.title('無声音のスペクトル')
plt.xlabel('周波数 [Hz]')
plt.ylabel('パワー [dB]')

plt.tight_layout()
plt.show()

5. 環境音の解析

目的：非定常な環境音を合成し、スペクトル特徴量・時間領域特徴量・リズム特徴量で解析する。環境音解析は、騒音監視や異常音検知に応用できる。

処理内容：雨音を大中小の雨滴（ポアソン分布に従う）でシミュレートし、スペクトル重心、RMSエネルギー、ゼロ交差率、テンポグラム、ノイズの色（スペクトル傾き）を計算・可視化する。

注意事項：環境音は非定常であるため、短時間分析で時間変化を捉える。スペクトル傾きが約-1ならピンクノイズ、約-2ならブラウンノイズに近い特性を示す。

出力の解釈：クレストファクターが大きいほどインパルス性が高い。尖度が正規分布（値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):
    """スペクトル特徴量の計算"""
    if np.all(y == 0):
        return {
            'centroid': np.zeros(1),
            'bandwidth': np.zeros(1),
            'rolloff': np.zeros(1),
            'flatness': np.zeros(1)
        }

    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):
    """時間領域特徴量の計算"""
    if np.all(y == 0):
        return {
            'rms': np.zeros(1),
            'zcr': np.zeros(1),
            'crest_factor': 0,
            'kurtosis': 0
        }

    # 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):
    """リズム特徴量の計算"""
    # テンポグラム
    onset_env = librosa.onset.onset_strength(y=y, sr=sr)
    tempogram = librosa.feature.tempogram(
        onset_envelope=onset_env,
        sr=sr
    )

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

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

def analyze_noise_color(y, sr):
    """ノイズの色分析（パワースペクトル密度の傾き）"""
    freqs, psd = welch(y, sr, nperseg=2048)
    # 対数変換
    log_freqs = np.log10(freqs[1:])  # 0Hzを除外
    log_psd = np.log10(psd[1:])
    # 線形回帰で傾きを計算
    slope, _, _, _, _ = scipy.stats.linregress(log_freqs, log_psd)
    return slope, freqs[1:], psd[1:]

def simulate_rain(duration=5.0, sr=44100):
    """雨音のシミュレーション"""
    t = np.linspace(0, duration, int(sr * duration))
    rain = np.zeros_like(t)

    # 雨滴の種類（大小）を考慮
    drop_types = {
        'small': {
            'freq_range': (2000, 4000),
            'duration_range': (0.01, 0.03),
            'amplitude_range': (0.1, 0.3)
        },
        'medium': {
            'freq_range': (1000, 2000),
            'duration_range': (0.03, 0.06),
            'amplitude_range': (0.3, 0.6)
        },
        'large': {
            'freq_range': (500, 1000),
            'duration_range': (0.06, 0.1),
            'amplitude_range': (0.6, 1.0)
        }
    }

    # 各種類の雨滴を生成
    for drop_type, params in drop_types.items():
        # 雨滴の数（ポアソン分布に従う）
        n_drops = np.random.poisson(duration *
                 (100 if drop_type == 'small' else
                  50 if drop_type == 'medium' else 20))

        for _ in range(n_drops):
            # 雨滴の時刻
            t_drop = np.random.random() * duration
            idx = int(t_drop * sr)

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

            # 雨滴の波形生成
            n_samples = int(drop_duration * sr)
            if idx + n_samples >= len(rain):
                continue

            # 減衰正弦波として雨滴を生成
            t_drop = np.arange(n_samples) / sr
            envelope = np.exp(-t_drop * 30)  # 指数減衰
            drop = amplitude * envelope * np.sin(2 * np.pi * freq * t_drop)

            # ランダムなフィルタリングで音色を変化
            drop = signal.lfilter([1], [1, -0.99], drop)

            rain[idx:idx+n_samples] += drop

    # クリッピングを防ぎつつ正規化
    rain = rain / (np.max(np.abs(rain)) + 1e-10)

    return rain, sr

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

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

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

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

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

# スペクトル特徴量の時間変化
plt.subplot(3, 2, 3)
times = np.linspace(0, len(y)/sr, len(spectral_features['centroid']))
plt.plot(times, librosa.amplitude_to_db(spectral_features['flatness']),
         label='平坦度')
plt.plot(times, spectral_features['centroid']/sr, label='重心')
plt.legend()
plt.title('スペクトル特徴量の時間変化')
plt.xlabel('時間 [秒]')

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

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

# エネルギーとゼロ交差率
plt.subplot(3, 2, 6)
times_rms = np.linspace(0, len(y)/sr, len(temporal_features['rms']))
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}")