piccal/plot_seq.py

##| Copyright: (C) 2019-2020 Kevin Larke <contact AT larke DOT org> 
##| License: GNU GPL version 3.0 or above. See the accompanying LICENSE file.
import os, sys
import matplotlib.pyplot as plt
import numpy as np
from common import parse_yaml_cfg
from rms_analysis import rms_analysis_main
from rms_analysis import select_first_stable_note_by_delta_db
from rms_analysis import select_first_stable_note_by_dur
from rms_analysis import samples_to_linear_residual

def is_nanV( xV ):
    
    for i in range(xV.shape[0]):
        if np.isnan( xV[i] ):
            return True
        
    return False

def _find_max_take_id( inDir ):

    id = 0
    while os.path.isdir( os.path.join(inDir, "%i" % id) ):
        id += 1

    if id > 0:
        id -= 1
        
    return id


def form_final_pulse_list( inDir, midi_pitch, analysisArgsD, take_id=None ):

    # append the midi pitch to the input directory
    #inDir = os.path.join( inDir, "%i" % (midi_pitch))

    dirL =  os.listdir(inDir)

    pkL = []

    # for each take in this directory
    for idir in dirL:

        take_number = int(idir)


        if not os.path.isfile(os.path.join( inDir,idir, "seq.json")):
            continue

        # analyze this takes audio and locate the note peaks
        r = rms_analysis_main( os.path.join(inDir,idir), midi_pitch, **analysisArgsD['rmsAnalysisArgs'] )

        # store the peaks in pkL[ (db,us) ]
        for db,us in zip(r.pkDbL,r.pkUsL):
            pkL.append( (db,us) )
            
    # sort the peaks on increasing attack pulse microseconds
    pkL = sorted( pkL, key= lambda x: x[1] )

    # merge sample points that separated by less than 'minSampleDistUs' milliseconds
    pkL = merge_close_sample_points( pkL, analysisArgsD['minSampleDistUs'] )
    
    # split pkL 
    pkDbL,pkUsL = tuple(zip(*pkL))

    #-------------------------------------------
        
    # locate the first and last note 
    min_pk_idx, max_pk_idx = find_min_max_peak_index( pkDbL, analysisArgsD['minAttkDb'], analysisArgsD['maxDbOffset'] )

    db1 = pkDbL[ max_pk_idx ]
    db0 = pkDbL[ min_pk_idx ]

    pulseUsL = []
    pulseDbL = []
    multValL = []
    for out_idx in range(128):

        # calc the target volume
        db = db0 + (out_idx * (db1-db0)/127.0)

        multi_value_count = 0

        # look for the target between each of the sampled points
        for i in range(1,len(pkDbL)):

            # if the target volume is between these two sample points
            if pkDbL[i-1] <= db and db < pkDbL[i]:

                # if the target has not already been located
                if len(pulseUsL) == out_idx:

                    # interpolate the pulse time from between the sampled points
                    frac = (db - pkDbL[i-1]) / (pkDbL[i] - pkDbL[i-1])
                    us = pkUsL[i-1] + frac * (pkUsL[i] - pkUsL[i-1])
                    db = pkDbL[i-1] + frac * (pkDbL[i] - pkDbL[i-1])
                    pulseUsL.append(us)
                    pulseDbL.append(db)

                else:
                    # this target db value was found between multiple sampled points
                    # therefore the sampled volume function is not monotonic
                    multi_value_count += 1

        if multi_value_count > 0:
            multValL.append((out_idx,multi_value_count))

    if len(multValL) > 0:
        # print("Multi-value pulse locations were found during velocity table formation: ",multValL)
        pass

    return pulseUsL,pulseDbL,r.holdDutyPctL
    


def merge_close_sample_points( pkDbUsL, minSampleDistanceUs ):

    avg0Us = np.mean(np.diff([ x[1] for x in pkDbUsL ]))
    n0 = len(pkDbUsL)
    
    while True and n0>0:
        us0 = None
        db0 = None
        
        for i,(db,us) in enumerate(pkDbUsL):
            if i > 0 and us - us0 < minSampleDistanceUs:
                us1 = (us0 + us)/2
                db1  = (db0 + db)/2
                pkDbUsL[i-1] = (db1,us1)
                del pkDbUsL[i]
                break
            else:
                us0 = us
                db0 = db

        if i+1 == len(pkDbUsL):
            break

    avg1Us = np.mean(np.diff([ x[1] for x in pkDbUsL ]))
    
    print("%i sample points deleted by merging close points." % (n0 - len(pkDbUsL)))
    print("Mean time between samples - before:%f after:%f " % (avg0Us,avg1Us))
    print("Min time between samples: %i " % (np.min(np.diff([x[1] for x in pkDbUsL]))))
    
    return pkDbUsL


def _calc_resample_points( dPkDb, pkUs0, pkUs1, samplePerDb, minSampleDistUs ):

    dPkUs = pkUs1 - pkUs0
    sampleCnt = max(int(round(abs(dPkDb) * samplePerDb)),samplePerDb)
    dUs       = max(int(round(dPkUs/sampleCnt)),minSampleDistUs)
    sampleCnt = int(round(dPkUs/dUs))
    dUs = int(round(dPkUs/sampleCnt))
    usL       = [ pkUs0 + dUs*j  for j in range(sampleCnt+1)]

    return usL

def calc_resample_ranges( pkDbL, pkUsL, min_pk_idx, max_pk_idx, maxDeltaDb, samplePerDb, minSampleDistUs ):

    if min_pk_idx == 0:
        print("No silent notes were generated.  Decrease the minimum peak level or the hold voltage.")
        return None

    resampleUsSet = set()
    refPkDb       = pkDbL[min_pk_idx]

    #pkDbL = pkDbL[ pkIdxL ]

    for i in range( min_pk_idx, max_pk_idx+1 ):

        d = pkDbL[i] - pkDbL[i-1]

        usL = []

        # if this peak is less than maxDeltaDb above the previous pk or
        # it is below the previous max peak
        if d > maxDeltaDb or d <= 0 or pkDbL[i] < refPkDb:

            usL = _calc_resample_points( d, pkUsL[i-1], pkUsL[i], samplePerDb, minSampleDistUs )
            
            if d <= 0 and i + 1 < len(pkDbL):
                d = pkDbL[i+1] - pkDbL[i]

                usL += _calc_resample_points( d, pkUsL[i-1], pkUsL[i], samplePerDb, minSampleDistUs )


        if pkDbL[i] > refPkDb:
            refPkDb = pkDbL[i]
                
        if usL:
            resampleUsSet = resampleUsSet.union( usL )

    return resampleUsSet
            


def form_resample_pulse_time_list( inDir, analysisArgsD ):
    """" This function merges all available data from previous takes to form
    a new list of pulse times to sample.
    """

    # the last folder is always the midi pitch of the note under analysis
    midi_pitch = int( inDir.split("/")[-1] )

    dirL =  os.listdir(inDir)

    pkL = []

    # for each take in this directory
    for idir in dirL:

        take_number = int(idir)

        # analyze this takes audio and locate the note peaks
        r = rms_analysis_main( os.path.join(inDir,idir), midi_pitch, **analysisArgsD['rmsAnalysisArgs'] )

        # store the peaks in pkL[ (db,us) ]
        for db,us in zip(r.pkDbL,r.pkUsL):
            pkL.append( (db,us) )
            
    # sort the peaks on increasing attack pulse microseconds
    pkL = sorted( pkL, key= lambda x: x[1] )

    # merge sample points that separated by less than 'minSampleDistUs' milliseconds
    pkL = merge_close_sample_points( pkL, analysisArgsD['minSampleDistUs'] )
    
    # split pkL 
    pkDbL,pkUsL = tuple(zip(*pkL))


    # locate the first and last note 
    min_pk_idx, max_pk_idx = find_min_max_peak_index( pkDbL, analysisArgsD['minAttkDb'], analysisArgsD['maxDbOffset'] )    

    # estimate the microsecond locations to resample
    resampleUsSet = calc_resample_ranges( pkDbL, pkUsL, min_pk_idx, max_pk_idx, analysisArgsD['maxDeltaDb'], analysisArgsD['samplesPerDb'], analysisArgsD['minSampleDistUs'] )

    resampleUsL = sorted( list(resampleUsSet) )

    #print(resampleUsL)

    return resampleUsL, pkDbL, pkUsL

def plot_curve( ax, pulseUsL, rmsDbV ):

    coeff = np.polyfit(pulseUsL,rmsDbV,5)
    func  = np.poly1d(coeff)

    ax.plot( pulseUsL, func(pulseUsL), color='red')

def plot_resample_pulse_times_0( inDir, analysisArgsD ):

    newPulseUsL, rmsDbV, pulseUsL = form_resample_pulse_time_list( inDir, analysisArgsD )

    midi_pitch = int( inDir.split("/")[-1] )
    velTblUsL,velTblDbL,_ = form_final_pulse_list( inDir, midi_pitch, analysisArgsD, take_id=None )
    
    fig,ax = plt.subplots()

    ax.plot(pulseUsL,rmsDbV,marker='.' )

    for us in newPulseUsL:
        ax.axvline( x = us )

    ax.plot(velTblUsL,velTblDbL,marker='.',linestyle='None')

    
    plt.show()

def plot_resample_pulse_times( inDir, analysisArgsD ):

    newPulseUsL, rmsDbV, pulseUsL = form_resample_pulse_time_list( inDir, analysisArgsD )

    midi_pitch = int( inDir.split("/")[-1] )
    velTblUsL,velTblDbL,_ = form_final_pulse_list( inDir, midi_pitch, analysisArgsD, take_id=None )
    
    fig,axL = plt.subplots(2,1,gridspec_kw={'height_ratios': [2, 1]})

    axL[0].plot(pulseUsL,rmsDbV,marker='.' )

    #plot_curve( ax, velTblUsL,velTblDbL)
    
    scoreV = samples_to_linear_residual( pulseUsL, rmsDbV)

    axL[0].plot(pulseUsL,rmsDbV + scoreV)
    axL[0].plot(pulseUsL,rmsDbV + np.power(scoreV,2.0))
    axL[0].plot(pulseUsL,rmsDbV - np.power(scoreV,2.0))

    axL[1].axhline(0.0,color='black')
    axL[1].axhline(1.0,color='black')
    axL[1].plot(pulseUsL,np.abs(scoreV * 100.0 / rmsDbV))
    axL[1].set_ylim((0.0,50))
    plt.show()



def find_min_max_peak_index( pkDbL, minDb, maxDbOffs ):
    """
    Find the min db and max db peak.
    """
    # select only the peaks from rmsV[] to work with
    yV = pkDbL

    # get the max volume note
    max_i = np.argmax( yV )
    maxDb = yV[ max_i ]

    min_i = max_i

    # starting from the max volume peak go backwards
    for i in range( max_i, 0, -1 ):

        # if this peak is within maxDbOffs of the loudest then choose this one instead
        #if maxDb - yV[i] < maxDbOffs:
        #    max_i = i

        # if this peak is less than minDb then the previous note is the min note
        if yV[i] < minDb:
            break
        
        min_i = i

    if min_i >= max_i:
        min_i = 0
        max_i = len(pkDbL)-1
    

    if min_i == 0:
        print("No silent notes were generated.  Decrease the minimum peak level or the hold voltage.")
    
    return min_i, max_i

def find_skip_peaks( rmsV, pkIdxL, min_pk_idx, max_pk_idx ):
    """ Fine peaks associated with longer attacks pulses that are lower than peaks with a shorter attack pulse.
    These peaks indicate degenerate portions of the pulse/db curve which must be skipped during velocity table formation
    """
    skipPkIdxL = []
    yV         = rmsV[pkIdxL]
    refPkDb    = yV[min_pk_idx]

    for i in range( min_pk_idx+1, max_pk_idx+1 ):
        if yV[i] > refPkDb:
            refPkDb = yV[i]
        else:
            skipPkIdxL.append(i)
            
            
    return skipPkIdxL

def find_out_of_range_peaks( rmsV, pkIdxL, min_pk_idx, max_pk_idx, maxDeltaDb ):
    """ Locate peaks which are more than maxDeltaDb from the previous peak. 
    If two peaks are separated by more than maxDeltaDb then the range must be resampled
    """

    oorPkIdxL = []
    yV         = rmsV[pkIdxL]

    for i in range( min_pk_idx, max_pk_idx+1 ):
        if i > 0:
            d = yV[i] - yV[i-1]
            if d > maxDeltaDb or d < 0:
                oorPkIdxL.append(i)

    return oorPkIdxL



def plot_spectrum( ax, srate, binHz, specV, midiPitch, harmN ):
    """ Plot a single spectrum, 'specV' and the harmonic peak location boundaries."""
    
    binN      = specV.shape[0]
    harmLBinL,harmMBinL,harmUBinL = calc_harm_bins( srate, binHz, midiPitch, harmN )

    fundHz      = harmMBinL[0] * binHz
    maxPlotHz   = fundHz * (harmN+1)
    maxPlotBinN = int(round(maxPlotHz/binHz))
    
    hzV = np.arange(binN) * (srate/(binN*2))
    
    specV = 20.0 * np.log10(specV) 
    
    ax.plot(hzV[0:maxPlotBinN], specV[0:maxPlotBinN] )

    for h0,h1,h2 in zip(harmLBinL,harmMBinL,harmUBinL):
        ax.axvline( x=h0 * binHz, color="blue")
        ax.axvline( x=h1 * binHz, color="black")
        ax.axvline( x=h2 * binHz, color="blue")

    ax.set_ylabel(str(midiPitch))
        
def plot_spectral_ranges( inDir, pitchL, rmsWndMs=300, rmsHopMs=30, harmN=5, dbLinRef=0.001 ):
    """ Plot the spectrum from one note (7th from last) in each attack pulse length sequence referred to by pitchL."""
    
    plotN = len(pitchL)
    fig,axL = plt.subplots(plotN,1)

    for plot_idx,midiPitch in enumerate(pitchL):

        # get the audio and meta-data file names
        seqFn   = os.path.join( inDir, str(midiPitch), "seq.json")
        audioFn = os.path.join( inDir, str(midiPitch), "audio.wav")

        # read the meta data object
        with open( seqFn, "rb") as f:
            r = json.load(f)

        # read the audio file
        srate, signalM  = wavfile.read(audioFn)
        sigV  = signalM / float(0x7fff)

        # calc. the RMS envelope in the time domain
        rms0DbV, rms0_srate = audio_rms( srate, sigV, rmsWndMs, rmsHopMs, dbLinRef )

        # locate the sample index of the peak of each note attack 
        pkIdx0L = locate_peak_indexes( rms0DbV, rms0_srate, r['eventTimeL'] )

        # select the 7th to last note for spectrum measurement

        #
        # TODO: come up with a better way to select the note to measure
        #
        spectrumSmpIdx = pkIdx0L[ len(pkIdx0L) - 7 ]


        # calc. the RMS envelope by taking the max spectral peak in each STFT window 
        rmsDbV, rms_srate, specV, specHopIdx, binHz = audio_stft_rms( srate, sigV, rmsWndMs, rmsHopMs, dbLinRef, spectrumSmpIdx)

        # specV[] is the spectrum of the note at spectrumSmpIdx

        # plot the spectrum and the harmonic selection ranges
        plot_spectrum( axL[plot_idx], srate, binHz, specV, midiPitch, harmN )
    plt.show()

        
        
def td_plot( ax, inDir, midi_pitch, id, analysisArgsD ):

    r = rms_analysis_main( inDir, midi_pitch, **analysisArgsD['rmsAnalysisArgs'] )

    min_pk_idx, max_pk_idx = find_min_max_peak_index( r.pkDbL, analysisArgsD['minAttkDb'], analysisArgsD['maxDbOffset'] )    
    
    skipPkIdxL = find_skip_peaks( r.rmsDbV, r.pkIdxL, min_pk_idx, max_pk_idx )

    jmpPkIdxL = find_out_of_range_peaks( r.rmsDbV, r.pkIdxL, min_pk_idx, max_pk_idx, analysisArgsD['maxDeltaDb'] )

    secV = np.arange(0,len(r.rmsDbV)) / r.rms_srate
    
    ax.plot( secV, r.rmsDbV )
    ax.plot( np.arange(0,len(r.tdRmsDbV)) / r.rms_srate, r.tdRmsDbV, color="black" )

    
    # print beg/end boundaries
    for i,(begMs, endMs) in enumerate(r.eventTimeL):

        pkSec = r.pkIdxL[i] / r.rms_srate
        endSec = pkSec + r.statsL[i].durMs / 1000.0
        
        ax.axvline( x=begMs/1000.0, color="green")
        ax.axvline( x=endMs/1000.0, color="red")
        ax.axvline( x=pkSec,        color="black")
        ax.axvline( x=endSec,       color="black")
        ax.text(begMs/1000.0, 20.0, str(i) )


    # plot peak markers
    for i,pki in enumerate(r.pkIdxL):
        marker = 4 if i==min_pk_idx or i==max_pk_idx else 5
        color  = "red" if i in skipPkIdxL else "black"
        ax.plot( [pki / r.rms_srate], [ r.rmsDbV[pki] ], marker=marker, color=color)

        if i in jmpPkIdxL:
            ax.plot( [pki / r.rms_srate], [ r.rmsDbV[pki] ], marker=6, color="blue")


    return r



def do_td_plot( inDir, analysisArgs ):
    
    fig,axL = plt.subplots(3,1)
    fig.set_size_inches(18.5, 10.5, forward=True)

    id         = int(inDir.split("/")[-1])
    midi_pitch = int(inDir.split("/")[-2])

    r = td_plot(axL[0],inDir,midi_pitch,id,analysisArgs)

    qualityV = np.array([ x.quality  for x in r.statsL ]) * np.max(r.pkDbL)
    durMsV   = np.array([ x.durMs    for x in r.statsL ])
    avgV     = np.array([ x.durAvgDb for x in r.statsL ])
    
    #durMsV[ durMsV < 400 ] = 0
    #durMsV = durMsV * np.max(r.pkDbL)/np.max(durMsV)
    #durMsV  = durMsV / 100.0

    dV = np.diff(r.pkDbL) / r.pkDbL[1:] 
    
    axL[1].plot( r.pkUsL, r.pkDbL, marker='.',label="pkDb" )
    axL[1].plot( r.pkUsL, qualityV, marker='.',label="quality" )
    axL[1].plot( r.pkUsL, avgV,     marker='.',label="avgDb" )
    #axL[2].plot( r.pkUsL, durMsV,   marker='.' )
    axL[2].plot( r.pkUsL[1:], dV,       marker='.',label='d')
    axL[2].set_ylim([-1,1])
    axL[1].legend()


    sni = select_first_stable_note_by_dur( durMsV )
    if sni is not None:
        axL[1].plot( r.pkUsL[sni], r.pkDbL[sni], marker='*', color='red')
        
    sni = select_first_stable_note_by_delta_db( r.pkDbL )
    if sni is not None:
        axL[2].plot( r.pkUsL[sni], dV[sni-1], marker='*', color='red')
    

    for i,s in enumerate(r.statsL):
        axL[1].text( r.pkUsL[i], r.pkDbL[i] + 1, "%i" % (i))

    for i in range(1,len(r.pkUsL)):
        axL[2].text( r.pkUsL[i],  dV[i-1], "%i" % (i))
        
    
    plt.show()

def do_td_multi_plot( inDir, analysisArgs ):

    midi_pitch = int(inDir.split("/")[-1])

    dirL =  os.listdir(inDir)

    fig,axL = plt.subplots(len(dirL),1)

    
    for id,(idir,ax) in enumerate(zip(dirL,axL)):

        td_plot(ax, os.path.join(inDir,str(id)), midi_pitch, id, analysisArgs )
    
    plt.show()

    
if __name__ == "__main__":

    inDir = sys.argv[1]
    cfgFn = sys.argv[2]
    take_id = None if len(sys.argv)<4 else sys.argv[3]

    cfg = parse_yaml_cfg( cfgFn )

    if take_id is not None:
        inDir = os.path.join(inDir,take_id)
        do_td_plot(inDir,cfg.analysisArgs)
    else:    
        #do_td_multi_plot(inDir,cfg.analysisArgs)

        #plot_spectral_ranges( inDir, [ 24, 36, 48, 60, 72, 84, 96, 104] )

        plot_resample_pulse_times( inDir, cfg.analysisArgs )