sound_plt_comp.py

#!/usr/bin/env python
import sys
import numpy as np
import matplotlib.pyplot as plt
import wave
import scipy.signal as signal
import math
from scipy.interpolate import interp1d
LIMIT_dB = 20

def downsampling(conversion_rate, data, fr) :
    decimation_sampleNum = conversion_rate-1
    nyqF = (fr/conversion_rate)/2.0
    cF = (fr/conversion_rate/2.0-500.)/nyqF
    taps = 511
    b = signal.firwin(taps, cF)
    data = signal.lfilter(b, 1, data)
    down_data = []
    for i in range(0, len(data), decimation_sampleNum+1):
        down_data.append(data[i])

    return (down_data, int(fr/conversion_rate))

def smoothing(input, window):
    output = []
    for i in range(input.shape[0]):
        if i < window:
            output.append(np.mean(input[:i+window+1]))
        elif i > input.shape[0] - 1 - window:
            output.append(np.mean(input[i:]))
        else:
            output.append(np.mean(input[i-window:i+window+1]))
    return np.array(output)

def to_db(data, N):
    pad = np.zeros(N//2)
    pad_data = np.concatenate([pad, data, pad])
    rms = np.array([np.sqrt((1/N) * (np.sum(pad_data[i:i+N]))**2) \
                    for i in range(len(data))])
    with np.errstate(divide='ignore'):
            db = 20 * np.log10(rms)
    return db

def valid_convolve(data, size):
    b = np.ones(size)/size
    data_mean = np.convolve(data, b, mode="same")
    n_conv = math.ceil(size/2)
    data_mean[0] *= size/n_conv
    for i in range(1, n_conv):
        data_mean[i] *= size/(i+n_conv)
        data_mean[-i] *= size/(i + n_conv - (size % 2)) 

    return data_mean

def sound_plt(wav_fname, wav_fname1, title):
        plt.rcParams['font.family'] = 'sans-serif'
        plt.rcParams['font.sans-serif'] = ['IPAPGothic', 'VL PGothic']
        fig = plt.figure(figsize=(16, 8), dpi=100, facecolor='tan', tight_layout=True)
        ax = fig.add_subplot(111, fc='w', xlabel='分(min)', ylabel='音量(dB)')
        ax.set_title(title + " 音量推移")

        wave_file = wave.open(wav_fname,"rb") 
        fr = wave_file.getframerate()
        nframes = wave_file.getnframes()
        data = wave_file.readframes(nframes) 
        data = np.frombuffer(data, dtype= "int16")
        down_rate = 2
        down_fr = int(fr / (down_rate * 1000))
        down_data, down_fr = downsampling(down_fr , data, fr)

        N = int(fr / 42)
        db = to_db(down_data, N)

        time = np.arange(0, db.shape[0] / down_fr, 1 / down_fr) / 60 / down_rate
        sm_db = smoothing(db, 1000)

        sm_db_x = [i for i in sm_db if i >= LIMIT_dB]
        db_mean = np.mean(sm_db_x)
        db_max = (np.max(sm_db_x))
        db_max_time = time[np.argmax(sm_db)]

        sm_db_s = [i for i in sm_db if i < LIMIT_dB]
        
        ax.plot(time, sm_db, 'c', label='No.1 signal')

        cnt = 0
        while True:
            sm_db1 = valid_convolve(sm_db, fr)
            time1 = np.linspace(np.min(time), np.max(time), np.size(sm_db1))
            sm_db = sm_db1
            time = time1
            cnt += 1
            if cnt >= 5:
                break
            
        ax.set_ylim(LIMIT_dB, 70)
        ax.plot(time, sm_db , 'b', label='No.1 moving average')
        ax.axhline(y=db_mean , color='b',linestyle='dashed', linewidth=1)
        ax.text(-0.1, db_mean-2, "No.1 average:" + str(round(db_mean,1)) + "dB", size=10)
        ax.axvline(x=db_max_time , color='b',linestyle='dashed', linewidth=1)
        ax.text(db_max_time, db_max+2, \
                "No.1 max:" + str(round(db_max,1)) + \
                'dB(' + str(round(db_max_time,1)) +'min)', \
                size=10)

        wave_file1 = wave.open(wav_fname1,"rb") 
        fr1 = wave_file1.getframerate()
        nframes1 = wave_file1.getnframes()
        data1 = wave_file1.readframes(nframes1) 
        data1 = np.frombuffer(data1, dtype= "int16")
        down_rate1 = 2
        down_fr1 = int(fr1 / (down_rate * 1000))
        down_data1, down_fr1 = downsampling(down_fr1 , data1, fr1)

        N = int(fr1 / 42)
        db1 = to_db(down_data1, N)

        time1 = np.arange(0, db1.shape[0] / down_fr1, 1 / down_fr1) / 60 / down_rate1
        sm_db1 = smoothing(db1, 1000)

        sm_db_x1 = [i for i in sm_db1 if i >= LIMIT_dB]
        db_mean1 = np.mean(sm_db_x1)
        db_max1 = (np.max(sm_db_x1))
        db_max_time1 = time1[np.argmax(sm_db1)]

        sm_db_s1 = [i for i in sm_db1 if i < LIMIT_dB]
        
        ax.plot(time1, sm_db1, 'm', label='No.2 Signal')

        cnt = 0
        while True:
            sm_db11 = valid_convolve(sm_db1, fr)
            time11 = np.linspace(np.min(time1), np.max(time1), np.size(sm_db11))
            sm_db1 = sm_db11
            time1 = time11
            cnt += 1
            if cnt >= 5:
                break
        ax.plot(time1, sm_db1 , 'r', label='No.2 moving average')
        ax.axhline(y=db_mean1 , color='r',linestyle='dashed', linewidth=1)
        ax.text(-0.1, db_mean1+1, "No.2 average:" + str(round(db_mean1,1)) + "dB", size=10)
        ax.axvline(x=db_max_time1 , color='r',linestyle='dashed', linewidth=1)
        ax.text(db_max_time1, db_max1+1, \
                "No.2 max:" + str(round(db_max1,1)) + \
                'dB(' + str(round(db_max_time1,1)) +'min)', \
                size=10)
        ax.grid()
        plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=0)
        fig.savefig(title + str(LIMIT_dB) + 'ret.png', facecolor=fig.get_facecolor())

        plt.show()

if __name__ == '__main__':
    wav_fname = sys.argv[1]
    wav_fname1 = sys.argv[2]
    title = sys.argv[3]
    sound_plt(wav_fname, wav_fname1, title)