音圧の移動平均微分を調べる

spl_sma_diff.py

import datetime
from os.path import exists
import pickle
from scipy.io.wavfile import read as read_wave
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import IndexLocator, FuncFormatter

mpl.rcParams["agg.path.chunksize"] = 100_000


class PCM:
    def __init__(self, title: str, file=None):
        self.title = title
        self.rate, self.data = read_wave(file or title + ".wav")
        self.samples, self.channels = self.data.shape
        self.duration = self.samples / self.rate
        self.quantize = self.data.dtype

def unit_convolve(array: np.ndarray, span: int) -> np.ndarray:
    length = array.shape[0]
    iteration = length - span + 1

    conv = np.empty(iteration)
    conv[0] = array[0:span].sum()
    for i in range(iteration - 1):
        conv[i+1] = conv[i] - array[i] + array[i+span]

    return conv

def sma(pcm: PCM, span: int, normalize=True) -> np.ndarray:
    width = pcm.rate * span // 1000
    monaural = pcm.data.sum(axis=1) / pcm.channels
    conv = unit_convolve(monaural**2, width)

    if normalize:
        average = np.sqrt(conv) / (width * np.iinfo(pcm.quantize).max)
    else:
        average = np.sqrt(conv) / width
    return average


# Milli-Seconds
SMA_SPAN = 200

# Names of source file
songs = [
    "未来の僕らは知ってるよ",
    "虹色Passions!",
]

# Load WAV Files and Calculate
data = {}
for title in songs:
    pcm = PCM(title)

    saved = f"{title}_sma{SMA_SPAN}"
    if exists(saved):
        with open(saved, 'rb') as f:
            dsma = pickle.load(f)
    else:
        print(f"Calculating Dsma: title={title}, SMA span={SMA_SPAN}")
        dsma = np.diff(sma(pcm, SMA_SPAN), n=1)
        with open(saved, 'wb') as f:
            pickle.dump(dsma, f)

    data[title] = {"PCM": pcm, "Dsma": dsma}


# Plot
def time_formatter(x, pos):
    h = int(x / 3600)
    m = int((x - h*3600) / 60)
    s = int(x - h*3600 - m*60)
    us = int((x - h*3600 - m*60 - s) * 1_000_000)
    time = datetime.time(h, m, s, us)
    return time.strftime("%M:%S")

def log_formatter(x, pos):
    dB = 10 * np.log10(abs(x))
    return r"$-\infty$" if np.isinf(dB) else "{:.0f}".format(dB)

fig, subplots = plt.subplots(len(songs), 1, figsize=(20,5*len(songs)), constrained_layout=True)

limit = [0, 0]
y_margin = 0.1
for d in data.values():
    limit[0] = max(limit[0], d["PCM"].duration)
    limit[1] = max(limit[1], d["Dsma"].max() * (1+y_margin), abs(d["Dsma"].min()) * (1+y_margin))

for i, title in enumerate(songs):
    print(f"Plotting: title={title}, SMA span={SMA_SPAN}")
    pcm = data[title]["PCM"]
    y = data[title]["Dsma"]
    x = np.linspace(0, pcm.duration, len(y))
    subplots[i].plot(x, y, linewidth=0.8)
    subplots[i].set_title(title, fontname="Yu Gothic")
    subplots[i].set_xlim(0, limit[0])
    subplots[i].set_ylim(-limit[1], limit[1])
    subplots[i].xaxis.set_major_locator(IndexLocator(15,0))
    subplots[i].xaxis.set_major_formatter(FuncFormatter(time_formatter))
    subplots[i].yaxis.set_major_formatter(FuncFormatter(log_formatter))
    subplots[i].set_ylabel("dB")
    subplots[i].grid(alpha=0.2)

fig.savefig(f"figure_{SMA_SPAN}.png", dpi=300, bbox_inches="tight")