短時間フーリエ変換（STFT）の理論とPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import stft, istft
from scipy.signal.windows import hann

def compute_stft_manual(x, window, hop_length, n_fft):
    """
    STFTの手動実装。

    Parameters
    ----------
    x : np.ndarray
        入力信号（1D）
    window : np.ndarray
        窓関数（長さ n_fft）
    hop_length : int
        フレームシフト（サンプル数）
    n_fft : int
        FFTの点数（= 窓長）

    Returns
    -------
    stft_matrix : np.ndarray (complex)
        形状 (n_fft//2 + 1, n_frames) の複素スペクトル行列
    """
    n_frames = 1 + (len(x) - n_fft) // hop_length
    stft_matrix = np.zeros((n_fft // 2 + 1, n_frames), dtype=complex)

    for m in range(n_frames):
        # フレームの切り出し
        start = m * hop_length
        frame = x[start : start + n_fft]
        # 窓関数の適用
        windowed_frame = frame * window
        # FFT（正の周波数成分のみ取得）
        stft_matrix[:, m] = np.fft.rfft(windowed_frame)

    return stft_matrix


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

# 時間変化する信号（チャープ + パルス）
# 0〜1秒: 100Hz → 500Hzに線形変化するチャープ
# 1〜2秒: 300Hzの定常正弦波
chirp = np.zeros_like(t)
segment1 = t < 1.0
chirp[segment1] = np.sin(2 * np.pi * (100 + 200 * t[segment1]) * t[segment1])
chirp[~segment1] = np.sin(2 * np.pi * 300 * t[~segment1])
# ホワイトノイズを加える
np.random.seed(42)
signal = chirp + 0.1 * np.random.randn(len(t))

# --- STFTのパラメータ ---
n_fft = 512
hop_length = n_fft // 4   # 75% オーバーラップ
window = hann(n_fft)

# --- 手動実装によるSTFT ---
stft_matrix = compute_stft_manual(signal, window, hop_length, n_fft)

# スペクトログラム（振幅）の計算
spectrogram = np.abs(stft_matrix)

# 時間・周波数軸の計算
n_frames = spectrogram.shape[1]
times = np.arange(n_frames) * hop_length / fs
freqs = np.fft.rfftfreq(n_fft, 1/fs)

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

# 元の信号
axes[0].plot(t, signal, linewidth=0.5)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Input Signal: Chirp (100→500 Hz) + 300 Hz')
axes[0].grid(True, alpha=0.3)

# スペクトログラム
img = axes[1].pcolormesh(
    times, freqs, 20 * np.log10(spectrogram + 1e-10),
    shading='gouraud', cmap='inferno', vmin=-60, vmax=0
)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Frequency [Hz]')
axes[1].set_title(f'Spectrogram (n_fft={n_fft}, hop={hop_length})')
axes[1].set_ylim(0, 800)
plt.colorbar(img, ax=axes[1], label='Magnitude [dB]')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import stft, istft
from scipy.signal.windows import hann

# --- テスト信号（音声系の模擬） ---
fs = 8000
t = np.arange(0, 3.0, 1/fs)

# 時間変化する信号を構成
# 0〜1秒: 200Hz
# 1〜2秒: 200Hz + 500Hz
# 2〜3秒: 800Hz
x = np.zeros_like(t)
mask1 = t < 1.0
mask2 = (t >= 1.0) & (t < 2.0)
mask3 = t >= 2.0
x[mask1] = np.sin(2 * np.pi * 200 * t[mask1])
x[mask2] = (np.sin(2 * np.pi * 200 * t[mask2])
            + 0.7 * np.sin(2 * np.pi * 500 * t[mask2]))
x[mask3] = np.sin(2 * np.pi * 800 * t[mask3])
x += 0.05 * np.random.randn(len(t))

# --- 窓長を変えてSTFT比較 ---
n_fft_list = [128, 512, 2048]
fig, axes = plt.subplots(len(n_fft_list), 1, figsize=(12, 10))

for ax, n_fft in zip(axes, n_fft_list):
    hop = n_fft // 4
    f, t_stft, Zxx = stft(x, fs=fs, window='hann',
                           nperseg=n_fft, noverlap=n_fft - hop)
    magnitude_dB = 20 * np.log10(np.abs(Zxx) + 1e-10)

    img = ax.pcolormesh(t_stft, f, magnitude_dB,
                        shading='gouraud', cmap='inferno',
                        vmin=-60, vmax=0)
    ax.set_ylabel('Frequency [Hz]')
    ax.set_ylim(0, 1200)
    df = fs / n_fft
    dt_ms = hop / fs * 1000
    ax.set_title(
        f'n_fft={n_fft}: Δf={df:.1f} Hz, Δt={dt_ms:.1f} ms '
        f'(time-freq tradeoff)'
    )
    plt.colorbar(img, ax=ax, label='dB')

axes[-1].set_xlabel('Time [s]')
plt.tight_layout()
plt.show()

from scipy.signal import stft, istft

# STFTの計算
n_fft = 512
hop = n_fft // 4
f, t_stft, Zxx = stft(x, fs=fs, window='hann',
                       nperseg=n_fft, noverlap=n_fft - hop)

# 逆STFTで再構成
_, x_reconstructed = istft(Zxx, fs=fs, window='hann',
                            nperseg=n_fft, noverlap=n_fft - hop)

# 再構成誤差の確認
min_len = min(len(x), len(x_reconstructed))
error = np.max(np.abs(x[:min_len] - x_reconstructed[:min_len]))
rms_error = np.sqrt(np.mean((x[:min_len] - x_reconstructed[:min_len]) ** 2))
print(f"最大再構成誤差: {error:.3e}")
print(f"RMS再構成誤差: {rms_error:.3e}")
# 実行結果:
# 最大再構成誤差: 8.882e-16
# RMS再構成誤差: 1.495e-16

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

# --- 音声的な合成信号の作成 ---
fs = 16000
duration = 1.5
t = np.arange(0, duration, 1/fs)

# 母音「あ」を模擬: F0=120Hz, フォルマント F1=800Hz, F2=1200Hz
np.random.seed(0)
glottal = np.zeros_like(t)
T0 = int(fs / 120)  # 基本周期
for i in range(0, len(t), T0):
    if i + 20 < len(t):
        glottal[i:i+20] = np.exp(-np.arange(20) * 0.3)  # 声帯パルス近似

# フォルマントフィルタ（簡易）
from scipy.signal import lfilter, butter

def formant_filter(signal, center_freq, bandwidth, fs):
    """単純なバンドパスフィルタでフォルマントを模擬"""
    b, a = butter(2, [center_freq - bandwidth/2, center_freq + bandwidth/2],
                  btype='band', fs=fs)
    return lfilter(b, a, signal)

voiced = (formant_filter(glottal, 800, 200, fs)
          + 0.7 * formant_filter(glottal, 1200, 150, fs)
          + 0.4 * formant_filter(glottal, 2500, 200, fs))
voiced = voiced / np.max(np.abs(voiced) + 1e-8)

# --- STFTとスペクトログラム ---
n_fft = 512
hop = n_fft // 8  # 87.5% オーバーラップ（音声処理では一般的）

f, t_stft, Zxx = stft(voiced, fs=fs, window='hann',
                       nperseg=n_fft, noverlap=n_fft - hop)

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(t, voiced, linewidth=0.5)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Synthetic Voiced Sound (vowel "a" approximation)')
axes[0].grid(True, alpha=0.3)

img = axes[1].pcolormesh(
    t_stft, f,
    20 * np.log10(np.abs(Zxx) + 1e-10),
    shading='gouraud', cmap='inferno', vmin=-60, vmax=0
)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Frequency [Hz]')
axes[1].set_ylim(0, 4000)
axes[1].set_title('Spectrogram: Formants visible at ~800 Hz, 1200 Hz, 2500 Hz')
plt.colorbar(img, ax=axes[1], label='dB')

plt.tight_layout()
plt.show()