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

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

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

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

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

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