librosa による音声スペクトログラムと CREPE によるF0時間変化（ソースコードと実行結果）

Python開発環境，ライブラリ類

ここでは、最低限の事前準備について説明する。機械学習や深層学習を行う場合は、NVIDIA CUDA、Visual Studio、Cursorなどを追加でインストールすると便利である。これらについては別ページ https://www.kkaneko.jp/cc/dev/aiassist.htmlで詳しく解説しているので、必要に応じて参照してください。

Python 3.12 のインストール

インストール済みの場合は実行不要。

管理者権限でコマンドプロンプトを起動（手順：Windowsキーまたはスタートメニュー > cmd と入力 > 右クリック > 「管理者として実行」）し、以下を実行する。管理者権限は、wingetの--scope machineオプションでシステム全体にソフトウェアをインストールするために必要である。

REM Python をシステム領域にインストール
winget install --scope machine --id Python.Python.3.12 -e --silent
REM Python のパス設定
set "PYTHON_PATH=C:\Program Files\Python312"
set "PYTHON_SCRIPTS_PATH=C:\Program Files\Python312\Scripts"
echo "%PATH%" | find /i "%PYTHON_PATH%" >nul
if errorlevel 1 setx PATH "%PATH%;%PYTHON_PATH%" /M >nul
echo "%PATH%" | find /i "%PYTHON_SCRIPTS_PATH%" >nul
if errorlevel 1 setx PATH "%PATH%;%PYTHON_SCRIPTS_PATH%" /M >nul

【関連する外部ページ】

Python の公式ページ: https://www.python.org/

AI エディタ Windsurf のインストール

Pythonプログラムの編集・実行には、AI エディタの利用を推奨する。ここでは，Windsurfのインストールを説明する。

管理者権限でコマンドプロンプトを起動（手順：Windowsキーまたはスタートメニュー > cmd と入力 > 右クリック > 「管理者として実行」）し、以下を実行して、Windsurfをシステム全体にインストールする。管理者権限は、wingetの--scope machineオプションでシステム全体にソフトウェアをインストールするために必要となる。

winget install --scope machine Codeium.Windsurf -e --silent

【関連する外部ページ】

Windsurf の公式ページ: https://windsurf.com/

必要なライブラリのインストール

コマンドプロンプトを管理者として実行（手順：Windowsキーまたはスタートメニュー > cmd と入力 > 右クリック > 「管理者として実行」）し、以下を実行する

pip install crepe librosa matplotlib numpy scipy sounddevice japanize-matplotlib tensorflow

librosa による音声スペクトログラムと CREPE によるF0時間変化プログラム

概要

このプログラムは音声から周波数スペクトログラムを計算し、深層学習モデルCREPEを用いて基本周波数（F0）を推定する。音声の物理的特性を数値化し、時間変化を可視化する。

主要技術

• CREPE（Convolutional Representation for Pitch Estimation）

  深層畳み込みニューラルネットワークを用いたピッチ推定技術である[1]。音声信号から直接基本周波数を推定し、従来の信号処理手法と比較して雑音環境下でも高い精度を維持する。1024サンプルの音声フレームを入力とし、360次元の確率分布として出力する。

• STFT（Short-Time Fourier Transform）

  音声信号を短時間窓で区切り、各窓でフーリエ変換を適用することで時間-周波数表現を得る手法である[2]。librosaライブラリの実装を使用し、2048点のFFTと512サンプルのホップ長で計算している。

参考文献

[1] Kim, J. W., Salamon, J., Li, P., & Bello, J. P. (2018). CREPE: A Convolutional Representation for Pitch Estimation. In 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (pp. 161-165). IEEE.

[2] McFee, B., Raffel, C., Liang, D., Ellis, D. P., McVicar, M., Battenberg, E., & Nieto, O. (2015). librosa: Audio and music signal analysis in python. In Proceedings of the 14th python in science conference (Vol. 8, pp. 18-25).

# プログラム名: librosa による音声スペクトログラムと CREPE による基本周波数（F0）時間変化解析プログラム
# 特徴技術名: CREPE (Convolutional Representation for Pitch Estimation)
# 出典: Kim, J. W., Salamon, J., Li, P., & Bello, J. P. (2018). CREPE: A convolutional representation for pitch estimation. In 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (pp. 161-165). IEEE. doi:10.1109/ICASSP.2018.8461329
# 特徴機能: 時間領域音声波形からの直接的ピッチ推定。深層畳み込みニューラルネットワークが音声波形を直接入力として受け取り、前処理なしに高精度な基本周波数推定を実現。従来の周波数領域解析手法と比較してノイズ耐性が高く、複雑な音響環境でも安定した性能を発揮
# 学習済みモデル: CREPE tiny/small/medium/large/full モデル（容量：tiny-0.9MB、small-1.8MB、medium-3.3MB、large-5.9MB、full-12MB）。GitHub公式リポジトリ（https://github.com/marl/crepe）から提供。TensorFlow SavedModel形式で、初回実行時にPython crepeライブラリを通じて自動ダウンロード
# 方式設計:
#   - 関連利用技術: librosa（音声読み込み・STFT・スペクトログラム変換）、matplotlib（可視化）、numpy（数値計算）、scipy（信号処理・フィルタリング）、sounddevice（音声再生）、tkinter（ファイル選択GUI）、urllib（サンプルファイルダウンロード）
#   - 入力と出力: 入力: 音声ファイル（ユーザは「0:音声ファイル，1:サンプル音声」のメニューで選択．0:音声ファイルの場合はtkinterでファイル選択．1の場合はhttps://github.com/librosa/librosa-test-data/raw/main/test1_44100.wavを使用）、出力: スペクトログラムとF0軌跡の重ね合わせ表示、F0統計情報のコンソール出力、1秒間隔でのF0値表示、result.txtファイルへの結果保存、音声再生機能
#   - 処理手順: 1.音声ファイル読み込み（librosa.load） 2.ハイパスフィルタ適用（scipy.butter, 50Hz以下除去） 3.音声信号正規化 4.STFT計算（librosa.stft） 5.スペクトログラム変換（amplitude_to_db） 6.CREPE F0推定（crepe.predict, Viterbiデコーディング適用） 7.信頼度フィルタリング 8.メディアンフィルタによるスムージング 9.F0統計計算 10.可視化・出力・保存
#   - 前処理、後処理: 前処理: バターワースハイパスフィルタ（5次、50Hz以下除去）による低域ノイズ除去、最大振幅による正規化、後処理: 5点メディアンフィルタによるF0軌跡平滑化、信頼度0.5未満の区間をNaN設定による無音区間処理
#   - 追加処理: 統計情報算出（平均・標準偏差・最小・最大・中央値・範囲）、1秒間隔サンプリングによる時系列F0値出力、結果の自動ファイル保存、処理日時記録、音声再生機能
#   - 調整を必要とする設定値: CREPE_MODEL_CAPACITY（モデルサイズ選択、精度と計算時間のトレードオフ）、CREPE_STEP_SIZE（時間解像度、デフォルト10ms）、HOP_LENGTH（STFT解像度、デフォルト512サンプル）、CONF_THRESH（無音判定閾値、デフォルト0.5）、MED_FILTER_SIZE（スムージング強度、デフォルト5）、HIGHPASS_CUTOFF（ノイズ除去周波数、デフォルト50Hz）
# 将来方策: 音声長に応じたHOP_LENGTHの動的調整機能実装（短時間音声では高解像度、長時間音声では処理効率を優先した自動パラメータ調整）
# その他の重要事項: Windows環境でのNumba/librosaキャッシュディレクトリ問題対応のため環境変数設定実装。CREPEモデルは初回のみダウンロード、以降はローカルキャッシュ使用。処理結果は日時付きでresult.txtに自動保存
# 前準備: pip install crepe librosa matplotlib numpy scipy sounddevice japanize-matplotlib tensorflow

import os
# Windowsでのキャッシュディレクトリ設定
os.environ['NUMBA_CACHE_DIR'] = os.path.join(os.path.expanduser('~'), '.numba_cache')
os.environ['LIBROSA_CACHE_DIR'] = os.path.join(os.path.expanduser('~'), '.librosa_cache')

import crepe
import librosa
import matplotlib.pyplot as plt
import japanize_matplotlib
import numpy as np
from scipy.signal import medfilt, butter, filtfilt
import tkinter as tk
from tkinter import filedialog
import urllib.request
import sounddevice as sd
from datetime import datetime

# 設定値（必要に応じて調整してください）
CREPE_MODEL_CAPACITY = 'tiny'  # CREPEモデルサイズ: 'tiny', 'small', 'medium', 'large', 'full'
CREPE_STEP_SIZE = 10  # CREPEの時間ステップ（ミリ秒）
HOP_LENGTH = 512  # 時間解像度（サンプル数）：小さいほど時間解像度が高い
N_FFT = 2048  # FFTウィンドウサイズ
WINDOW_TYPE = 'hann'  # 窓関数の種類
CONF_THRESH = 0.5  # 無音判定の信頼度閾値（0.0-1.0）
MED_FILTER_SIZE = 5  # メディアンフィルタのカーネルサイズ
SAMPLE_AUDIO_URL = 'https://github.com/librosa/librosa-test-data/raw/main/test1_44100.wav'
SPECTROGRAM_MAX_FREQ = 2000  # スペクトログラム表示の最大周波数（Hz）
F0_MAX_FREQ = 800  # F0表示の最大周波数（Hz）
HIGHPASS_CUTOFF = 50  # ハイパスフィルタのカットオフ周波数（Hz）
FILTER_ORDER = 5  # バターワースフィルタの次数
RESULT_FILENAME = 'result.txt'  # 結果保存ファイル名

# 表示関連定数
FIGURE_SIZE_WIDTH = 12
FIGURE_SIZE_HEIGHT = 6
LINE_WIDTH = 4
# 統計表示関連定数
PRINT_INTERVAL_SEC = 1  # 1秒間隔表示用

# ファイル選択関連
AUDIO_FILE_TYPES = [('Audio Files', '*.wav *.mp3 *.flac *.m4a *.ogg')]
FILE_DIALOG_TITLE = '音声ファイルを選択'
TEMP_FILENAME = 'sample_audio.wav'

# プログラムの説明表示
print('=== 音声スペクトログラムとF0時間変化表示プログラム ===')
print('このプログラムは音声ファイルから以下の処理を行います：')
print('1. スペクトログラム（周波数成分の時間変化）の計算')
print('2. 基本周波数（F0/ピッチ）の推定')
print('3. スペクトログラムとF0の重ね合わせ表示')
print('4. F0統計情報の計算と表示')
print('5. 結果をresult.txtに自動保存')
print()

# 音声ファイル選択
print('0: 音声ファイル')
print('1: サンプル音声')

choice = input('選択: ')
temp_file = None
result_data = []  # 結果保存用

try:
    if choice == '0':
        root = tk.Tk()
        root.withdraw()
        audio_path = filedialog.askopenfilename(
            title=FILE_DIALOG_TITLE,
            filetypes=AUDIO_FILE_TYPES
        )
        root.destroy()
        if not audio_path:
            print('ファイルが選択されませんでした。')
            exit()
        print(f'選択されたファイル: {audio_path}')
    elif choice == '1':
        # サンプル音声ダウンロード
        url = SAMPLE_AUDIO_URL
        filename = TEMP_FILENAME
        try:
            print('サンプル音声をダウンロードしています...')
            urllib.request.urlretrieve(url, filename)
            temp_file = filename
            audio_path = filename
            print('サンプル音声のダウンロードが完了しました。')
        except Exception as e:
            print(f'音声のダウンロードに失敗しました: {url}')
            print(f'エラー: {e}')
            exit()
    else:
        print('無効な選択です')
        exit()

    # 結果データの記録開始
    result_data.append(f'処理日時: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}')
    result_data.append(f'音声ファイル: {audio_path}')
    result_data.append('')

    # 音声ファイル読み込み
    print('音声ファイル読み込み中...')
    print(f'ファイルパス: {audio_path}')

    try:
        audio, sr = librosa.load(audio_path, sr=None, mono=True)
        print(f'読み込み成功: サンプリングレート={sr}Hz, 長さ={len(audio)/sr:.1f}秒')
        result_data.append(f'サンプリングレート: {sr} Hz')
        result_data.append(f'音声長: {len(audio)/sr:.1f} 秒')
        result_data.append('')
    except Exception as e:
        print(f'音声ファイルの読み込みに失敗しました: {e}')
        exit()

    # 音声信号の前処理
    print('音声信号前処理中...')

    # ハイパスフィルタの適用（50Hz以下をカット）
    nyquist = sr / 2
    if HIGHPASS_CUTOFF >= nyquist:
        print(f'警告: ハイパスフィルタのカットオフ周波数({HIGHPASS_CUTOFF}Hz)がナイキスト周波数({nyquist}Hz)以上です。フィルタを適用しません。')
        audio_filt = audio
    else:
        normal_cutoff = HIGHPASS_CUTOFF / nyquist
        b, a = butter(FILTER_ORDER, normal_cutoff, btype='high', analog=False)
        audio_filt = filtfilt(b, a, audio)

    # 正規化
    if np.max(np.abs(audio_filt)) > 0:
        audio_norm = audio_filt / np.max(np.abs(audio_filt))
    else:
        print('警告: 音声信号の振幅が0です。')
        audio_norm = audio_filt

    # スペクトログラム計算
    print('スペクトログラム計算中...')

    # librosaによるSTFT計算
    D = librosa.stft(audio_norm, n_fft=N_FFT, hop_length=HOP_LENGTH, window=WINDOW_TYPE)
    S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)

    # CREPEによるF0推定
    print('F0推定中...')
    print(f'音声データ: shape={audio_norm.shape}, sr={sr}')
    try:
        time, frequency, confidence, activation = crepe.predict(
            audio_norm, sr, viterbi=True,
            model_capacity=CREPE_MODEL_CAPACITY,
            step_size=CREPE_STEP_SIZE
        )
        print('F0推定完了')
    except Exception as e:
        print(f'CREPE F0推定に失敗しました: {e}')
        exit()

    # 無音区間のF0をNaNに設定（グラフで不連続表示）
    frequency[confidence < CONF_THRESH] = np.nan

    # F0のスムージング（メディアンフィルタ、NaNは保持）
    freq_smooth = frequency.copy()
    valid_mask = ~np.isnan(frequency)
    if np.any(valid_mask):
        # 有効な値のみでメディアンフィルタを適用
        valid_indices = np.where(valid_mask)[0]
        for i in valid_indices:
            start = max(0, i - MED_FILTER_SIZE // 2)
            end = min(len(frequency), i + MED_FILTER_SIZE // 2 + 1)
            window_data = frequency[start:end]
            valid_window = window_data[~np.isnan(window_data)]
            if len(valid_window) > 0:
                freq_smooth[i] = np.median(valid_window)

    # F0統計情報の計算（簡素化）
    valid_f0 = freq_smooth[~np.isnan(freq_smooth)]
    if len(valid_f0) > 0:
        stats = {
            '平均F0': np.mean(valid_f0),
            '標準偏差': np.std(valid_f0),
            '最小F0': np.min(valid_f0),
            '最大F0': np.max(valid_f0),
            '中央値': np.median(valid_f0)
        }
        stats['F0範囲'] = stats['最大F0'] - stats['最小F0']
    else:
        stats = {key: 0 for key in ['平均F0', '標準偏差', '最小F0', '最大F0', '中央値', 'F0範囲']}

    # グラフ描画
    fig, ax1 = plt.subplots(figsize=(FIGURE_SIZE_WIDTH, FIGURE_SIZE_HEIGHT))

    # スペクトログラム表示
    freqs = librosa.fft_frequencies(sr=sr, n_fft=N_FFT)
    times = librosa.frames_to_time(np.arange(S_db.shape[1]), sr=sr, hop_length=HOP_LENGTH)
    img = ax1.pcolormesh(times, freqs, S_db, shading='auto', cmap='viridis')
    ax1.set_ylim(0, SPECTROGRAM_MAX_FREQ)

    # F0グラフの重ね合わせ（水平線分として表示）
    ax2 = ax1.twinx()

    # 各測定点で水平線分を描画
    step_size_sec = CREPE_STEP_SIZE / 1000.0  # ミリ秒を秒に変換
    for i in range(len(time)):
        if not np.isnan(freq_smooth[i]):
            # 各測定点を中心とした水平線分
            x_start = time[i] - step_size_sec / 2
            x_end = time[i] + step_size_sec / 2
            ax2.plot([x_start, x_end], [freq_smooth[i], freq_smooth[i]],
                    'r-', linewidth=LINE_WIDTH)

    ax2.set_ylabel('Frequency (Hz)', color='r')
    ax2.tick_params(axis='y', labelcolor='r')
    ax2.set_ylim(0, F0_MAX_FREQ)

    # タイトルとラベル
    ax1.set_xlabel('Time (s)')
    ax1.set_ylabel('Frequency (Hz)')
    plt.title('Spectrogram with F0 Contour')

    # カラーバー追加
    cbar = plt.colorbar(img, ax=ax1, format='%+2.0f dB')

    # F0統計情報の出力と記録（簡素化）
    duration = len(audio) / sr
    print(f'\n音声長: {duration:.1f}秒')
    print('\nF0統計情報:')

    result_data.append('F0統計情報:')
    for key, value in stats.items():
        print(f'{key}: {value:.1f} Hz')
        result_data.append(f'{key}: {value:.1f} Hz')
    result_data.append('')

    # 1秒ごとのF0値出力と記録
    print('\n1秒ごとのF0値:')
    result_data.append('1秒ごとのF0値:')
    for t in range(0, int(duration) + 1, PRINT_INTERVAL_SEC):
        if t < len(time):
            idx = np.argmin(np.abs(time - t))
            if not np.isnan(freq_smooth[idx]):
                print(f'{t:3d}秒: {freq_smooth[idx]:.1f} Hz')
                result_data.append(f'{t:3d}秒: {freq_smooth[idx]:.1f} Hz')
            else:
                print(f'{t:3d}秒: 無音')
                result_data.append(f'{t:3d}秒: 無音')

    # 結果をファイルに保存
    try:
        with open(RESULT_FILENAME, 'w', encoding='utf-8') as f:
            f.write('\n'.join(result_data))
        print(f'\n結果を {RESULT_FILENAME} に保存しました')
    except Exception as e:
        print(f'結果ファイルの保存に失敗しました: {e}')

    print('保存内容:')
    print('- 処理日時と音声ファイル情報')
    print('- F0統計情報（平均、標準偏差、最小、最大、中央値、範囲）')
    print('- 1秒ごとのF0値')

    plt.tight_layout()
    plt.show()

    # 音声再生
    print('\n音声を再生しますか？')
    print('y: 再生')
    print('n: 終了')
    play_choice = input('選択: ')

    if play_choice.lower() == 'y':
        print('音声再生中... (Ctrl+Cで中断)')
        try:
            sd.play(audio, sr)
            sd.wait()
            print('再生完了')
        except KeyboardInterrupt:
            sd.stop()
            print('\n再生を中断しました')
        except Exception as e:
            print(f'音声再生エラー: {e}')

except Exception as e:
    print(f'予期しないエラーが発生しました: {e}')
finally:
    # 一時ファイルの削除
    if temp_file and os.path.exists(temp_file):
        try:
            os.remove(temp_file)
            print('一時ファイルを削除しました。')
        except OSError:
            print('一時ファイルの削除に失敗しました。')