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 )