CDイメージからPCMデータを抽出して離散図・MS波形・スペクトルを描画

ms_plot.py

import numpy as np
from scipy import signal
from scipy import fftpack
from scipy import stats
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

samplerate = 44100
fft_size = 8192

start_time = [15, 30, 7]  # [min, sec, 0.01sec]
stop_time = [18, 54, 67]
start_sec = (start_time[0] * 60) + start_time[1] + (start_time[2] * 0.01)
stop_sec = (stop_time[0] * 60) + stop_time[1] + (stop_time[2] * 0.01)
start_sample = int(start_sec * samplerate)
stop_sample = int(stop_sec * samplerate)

num_sample = stop_sample - start_sample
num_seek = start_sample * 2 * 2  # 16bit, 2ch

time_min = np.linspace(0., (stop_sec-start_sec)/60, num_sample, endpoint=False)

file = open("Image.bin", mode="rb")
file.seek(num_seek)

data = np.fromfile(file, "<i2", num_sample * 2)  # 16bit, 2ch
data = data.reshape((-1,2))

data_L = data[:,0].astype(np.int32)
data_R = data[:,1].astype(np.int32)
del(data)

data_mid = ((data_L + data_R) * 0.5).astype(np.int32)
data_side = ((data_L - data_R) * 0.5).astype(np.int32)

r, _ = stats.pearsonr(data_L, data_R)
print(r)
fig = plt.figure(figsize=(8, 8))
plt.scatter(data_L/32768., data_R/32768., alpha=0.01, color="black", s=1)
plt.xlim([-1., 1])
plt.ylim([-1., 1])
plt.xlabel("L channel amplitude")
plt.ylabel("R channel amplitude")
plt.text(-0.8, 0.8, "r={:.3f}".format(r))
plt.savefig("graph_phase.png")
del(fig)
del(data_L, data_R)

data_mid_nor = data_mid / 32768.
data_side_nor = data_side / 32768.

fig = plt.figure(figsize=(8, 5))
plt.plot(time_min, data_mid_nor, label="Mid", linewidth=0.5)
plt.plot(time_min, data_side_nor, label="Side", linewidth=0.5)
plt.xlabel("Time [min]")
plt.ylabel("Amplitude")
plt.legend(loc="upper left")
plt.savefig("graph_ms.png")
del(fig)

remain = num_sample % fft_size
count_fft = 2 * int((num_sample + remain) / fft_size) - 1
hann_window = signal.hann(fft_size)
freq_space = fftpack.fftfreq(n=fft_size, d=1./samplerate)[0:int(fft_size/2)]

data_mid_argnd = np.concatenate([data_mid_nor, np.zeros(remain)])
data_side_argnd = np.concatenate([data_side_nor, np.zeros(remain)])

def calc_fft(indata):
    f_average = np.zeros(fft_size)
    f_peak = np.zeros(fft_size)
    for c in range(count_fft):
        s_smpl = c * int(fft_size/2)
        e_smpl = s_smpl + fft_size
        d = indata[s_smpl:e_smpl] * hann_window
        f = np.abs(fftpack.fft(d))
        f_average += f
        f_peak = np.max([f, f_peak], axis=0)

    f_average /= count_fft
    return f_peak[0:int(fft_size/2)], f_average[0:int(fft_size/2)]

fft_mid_peak, fft_mid_average = calc_fft(data_mid_argnd)
fft_side_peak, fft_side_average = calc_fft(data_side_argnd)

fig = plt.figure(figsize=(8, 5))
plt.plot(freq_space, fft_mid_peak, label="Mid peak", color="C0", linewidth=0.75)
plt.plot(freq_space, fft_mid_average, label="Mid average", color="C0", linestyle="--", linewidth=0.75)
plt.plot(freq_space, fft_side_peak, label="Side peak", color="C1", linewidth=0.75)
plt.plot(freq_space, fft_side_average, label="Side average", color="C1", linestyle="--", linewidth=0.75)
plt.xscale('log')
# plt.yscale('log')
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude")
plt.legend(loc="upper right")
plt.savefig("graph_ms_fft.png")