picadae calibration programs
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

plot_seq.py 17KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559
  1. ##| Copyright: (C) 2019-2020 Kevin Larke <contact AT larke DOT org>
  2. ##| License: GNU GPL version 3.0 or above. See the accompanying LICENSE file.
  3. import os, sys
  4. import matplotlib.pyplot as plt
  5. import numpy as np
  6. from common import parse_yaml_cfg
  7. from rms_analysis import rms_analysis_main
  8. from rms_analysis import select_first_stable_note_by_delta_db
  9. from rms_analysis import select_first_stable_note_by_dur
  10. from rms_analysis import samples_to_linear_residual
  11. def is_nanV( xV ):
  12. for i in range(xV.shape[0]):
  13. if np.isnan( xV[i] ):
  14. return True
  15. return False
  16. def _find_max_take_id( inDir ):
  17. id = 0
  18. while os.path.isdir( os.path.join(inDir, "%i" % id) ):
  19. id += 1
  20. if id > 0:
  21. id -= 1
  22. return id
  23. def form_final_pulse_list( inDir, midi_pitch, analysisArgsD, take_id=None ):
  24. # append the midi pitch to the input directory
  25. #inDir = os.path.join( inDir, "%i" % (midi_pitch))
  26. dirL = os.listdir(inDir)
  27. pkL = []
  28. # for each take in this directory
  29. for idir in dirL:
  30. take_number = int(idir)
  31. if not os.path.isfile(os.path.join( inDir,idir, "seq.json")):
  32. continue
  33. # analyze this takes audio and locate the note peaks
  34. r = rms_analysis_main( os.path.join(inDir,idir), midi_pitch, **analysisArgsD['rmsAnalysisArgs'] )
  35. # store the peaks in pkL[ (db,us) ]
  36. for db,us in zip(r.pkDbL,r.pkUsL):
  37. pkL.append( (db,us) )
  38. # sort the peaks on increasing attack pulse microseconds
  39. pkL = sorted( pkL, key= lambda x: x[1] )
  40. # merge sample points that separated by less than 'minSampleDistUs' milliseconds
  41. pkL = merge_close_sample_points( pkL, analysisArgsD['minSampleDistUs'] )
  42. # split pkL
  43. pkDbL,pkUsL = tuple(zip(*pkL))
  44. #-------------------------------------------
  45. # locate the first and last note
  46. min_pk_idx, max_pk_idx = find_min_max_peak_index( pkDbL, analysisArgsD['minAttkDb'], analysisArgsD['maxDbOffset'] )
  47. db1 = pkDbL[ max_pk_idx ]
  48. db0 = pkDbL[ min_pk_idx ]
  49. pulseUsL = []
  50. pulseDbL = []
  51. multValL = []
  52. for out_idx in range(128):
  53. # calc the target volume
  54. db = db0 + (out_idx * (db1-db0)/127.0)
  55. multi_value_count = 0
  56. # look for the target between each of the sampled points
  57. for i in range(1,len(pkDbL)):
  58. # if the target volume is between these two sample points
  59. if pkDbL[i-1] <= db and db < pkDbL[i]:
  60. # if the target has not already been located
  61. if len(pulseUsL) == out_idx:
  62. # interpolate the pulse time from between the sampled points
  63. frac = (db - pkDbL[i-1]) / (pkDbL[i] - pkDbL[i-1])
  64. us = pkUsL[i-1] + frac * (pkUsL[i] - pkUsL[i-1])
  65. db = pkDbL[i-1] + frac * (pkDbL[i] - pkDbL[i-1])
  66. pulseUsL.append(us)
  67. pulseDbL.append(db)
  68. else:
  69. # this target db value was found between multiple sampled points
  70. # therefore the sampled volume function is not monotonic
  71. multi_value_count += 1
  72. if multi_value_count > 0:
  73. multValL.append((out_idx,multi_value_count))
  74. if len(multValL) > 0:
  75. # print("Multi-value pulse locations were found during velocity table formation: ",multValL)
  76. pass
  77. return pulseUsL,pulseDbL,r.holdDutyPctL
  78. def merge_close_sample_points( pkDbUsL, minSampleDistanceUs ):
  79. avg0Us = np.mean(np.diff([ x[1] for x in pkDbUsL ]))
  80. n0 = len(pkDbUsL)
  81. while True and n0>0:
  82. us0 = None
  83. db0 = None
  84. for i,(db,us) in enumerate(pkDbUsL):
  85. if i > 0 and us - us0 < minSampleDistanceUs:
  86. us1 = (us0 + us)/2
  87. db1 = (db0 + db)/2
  88. pkDbUsL[i-1] = (db1,us1)
  89. del pkDbUsL[i]
  90. break
  91. else:
  92. us0 = us
  93. db0 = db
  94. if i+1 == len(pkDbUsL):
  95. break
  96. avg1Us = np.mean(np.diff([ x[1] for x in pkDbUsL ]))
  97. print("%i sample points deleted by merging close points." % (n0 - len(pkDbUsL)))
  98. print("Mean time between samples - before:%f after:%f " % (avg0Us,avg1Us))
  99. print("Min time between samples: %i " % (np.min(np.diff([x[1] for x in pkDbUsL]))))
  100. return pkDbUsL
  101. def _calc_resample_points( dPkDb, pkUs0, pkUs1, samplePerDb, minSampleDistUs ):
  102. dPkUs = pkUs1 - pkUs0
  103. sampleCnt = max(int(round(abs(dPkDb) * samplePerDb)),samplePerDb)
  104. dUs = max(int(round(dPkUs/sampleCnt)),minSampleDistUs)
  105. sampleCnt = int(round(dPkUs/dUs))
  106. dUs = int(round(dPkUs/sampleCnt))
  107. usL = [ pkUs0 + dUs*j for j in range(sampleCnt+1)]
  108. return usL
  109. def calc_resample_ranges( pkDbL, pkUsL, min_pk_idx, max_pk_idx, maxDeltaDb, samplePerDb, minSampleDistUs ):
  110. if min_pk_idx == 0:
  111. print("No silent notes were generated. Decrease the minimum peak level or the hold voltage.")
  112. return None
  113. resampleUsSet = set()
  114. refPkDb = pkDbL[min_pk_idx]
  115. #pkDbL = pkDbL[ pkIdxL ]
  116. for i in range( min_pk_idx, max_pk_idx+1 ):
  117. d = pkDbL[i] - pkDbL[i-1]
  118. usL = []
  119. # if this peak is less than maxDeltaDb above the previous pk or
  120. # it is below the previous max peak
  121. if d > maxDeltaDb or d <= 0 or pkDbL[i] < refPkDb:
  122. usL = _calc_resample_points( d, pkUsL[i-1], pkUsL[i], samplePerDb, minSampleDistUs )
  123. if d <= 0 and i + 1 < len(pkDbL):
  124. d = pkDbL[i+1] - pkDbL[i]
  125. usL += _calc_resample_points( d, pkUsL[i-1], pkUsL[i], samplePerDb, minSampleDistUs )
  126. if pkDbL[i] > refPkDb:
  127. refPkDb = pkDbL[i]
  128. if usL:
  129. resampleUsSet = resampleUsSet.union( usL )
  130. return resampleUsSet
  131. def form_resample_pulse_time_list( inDir, analysisArgsD ):
  132. """" This function merges all available data from previous takes to form
  133. a new list of pulse times to sample.
  134. """
  135. # the last folder is always the midi pitch of the note under analysis
  136. midi_pitch = int( inDir.split("/")[-1] )
  137. dirL = os.listdir(inDir)
  138. pkL = []
  139. # for each take in this directory
  140. for idir in dirL:
  141. take_number = int(idir)
  142. # analyze this takes audio and locate the note peaks
  143. r = rms_analysis_main( os.path.join(inDir,idir), midi_pitch, **analysisArgsD['rmsAnalysisArgs'] )
  144. # store the peaks in pkL[ (db,us) ]
  145. for db,us in zip(r.pkDbL,r.pkUsL):
  146. pkL.append( (db,us) )
  147. # sort the peaks on increasing attack pulse microseconds
  148. pkL = sorted( pkL, key= lambda x: x[1] )
  149. # merge sample points that separated by less than 'minSampleDistUs' milliseconds
  150. pkL = merge_close_sample_points( pkL, analysisArgsD['minSampleDistUs'] )
  151. # split pkL
  152. pkDbL,pkUsL = tuple(zip(*pkL))
  153. # locate the first and last note
  154. min_pk_idx, max_pk_idx = find_min_max_peak_index( pkDbL, analysisArgsD['minAttkDb'], analysisArgsD['maxDbOffset'] )
  155. # estimate the microsecond locations to resample
  156. resampleUsSet = calc_resample_ranges( pkDbL, pkUsL, min_pk_idx, max_pk_idx, analysisArgsD['maxDeltaDb'], analysisArgsD['samplesPerDb'], analysisArgsD['minSampleDistUs'] )
  157. resampleUsL = sorted( list(resampleUsSet) )
  158. #print(resampleUsL)
  159. return resampleUsL, pkDbL, pkUsL
  160. def plot_curve( ax, pulseUsL, rmsDbV ):
  161. coeff = np.polyfit(pulseUsL,rmsDbV,5)
  162. func = np.poly1d(coeff)
  163. ax.plot( pulseUsL, func(pulseUsL), color='red')
  164. def plot_resample_pulse_times_0( inDir, analysisArgsD ):
  165. newPulseUsL, rmsDbV, pulseUsL = form_resample_pulse_time_list( inDir, analysisArgsD )
  166. midi_pitch = int( inDir.split("/")[-1] )
  167. velTblUsL,velTblDbL,_ = form_final_pulse_list( inDir, midi_pitch, analysisArgsD, take_id=None )
  168. fig,ax = plt.subplots()
  169. ax.plot(pulseUsL,rmsDbV,marker='.' )
  170. for us in newPulseUsL:
  171. ax.axvline( x = us )
  172. ax.plot(velTblUsL,velTblDbL,marker='.',linestyle='None')
  173. plt.show()
  174. def plot_resample_pulse_times( inDir, analysisArgsD ):
  175. newPulseUsL, rmsDbV, pulseUsL = form_resample_pulse_time_list( inDir, analysisArgsD )
  176. midi_pitch = int( inDir.split("/")[-1] )
  177. velTblUsL,velTblDbL,_ = form_final_pulse_list( inDir, midi_pitch, analysisArgsD, take_id=None )
  178. fig,axL = plt.subplots(2,1,gridspec_kw={'height_ratios': [2, 1]})
  179. axL[0].plot(pulseUsL,rmsDbV,marker='.' )
  180. #plot_curve( ax, velTblUsL,velTblDbL)
  181. scoreV = samples_to_linear_residual( pulseUsL, rmsDbV)
  182. axL[0].plot(pulseUsL,rmsDbV + scoreV)
  183. axL[0].plot(pulseUsL,rmsDbV + np.power(scoreV,2.0))
  184. axL[0].plot(pulseUsL,rmsDbV - np.power(scoreV,2.0))
  185. axL[1].axhline(0.0,color='black')
  186. axL[1].axhline(1.0,color='black')
  187. axL[1].plot(pulseUsL,np.abs(scoreV * 100.0 / rmsDbV))
  188. axL[1].set_ylim((0.0,50))
  189. plt.show()
  190. def find_min_max_peak_index( pkDbL, minDb, maxDbOffs ):
  191. """
  192. Find the min db and max db peak.
  193. """
  194. # select only the peaks from rmsV[] to work with
  195. yV = pkDbL
  196. # get the max volume note
  197. max_i = np.argmax( yV )
  198. maxDb = yV[ max_i ]
  199. min_i = max_i
  200. # starting from the max volume peak go backwards
  201. for i in range( max_i, 0, -1 ):
  202. # if this peak is within maxDbOffs of the loudest then choose this one instead
  203. #if maxDb - yV[i] < maxDbOffs:
  204. # max_i = i
  205. # if this peak is less than minDb then the previous note is the min note
  206. if yV[i] < minDb:
  207. break
  208. min_i = i
  209. if min_i >= max_i:
  210. min_i = 0
  211. max_i = len(pkDbL)-1
  212. if min_i == 0:
  213. print("No silent notes were generated. Decrease the minimum peak level or the hold voltage.")
  214. return min_i, max_i
  215. def find_skip_peaks( rmsV, pkIdxL, min_pk_idx, max_pk_idx ):
  216. """ Fine peaks associated with longer attacks pulses that are lower than peaks with a shorter attack pulse.
  217. These peaks indicate degenerate portions of the pulse/db curve which must be skipped during velocity table formation
  218. """
  219. skipPkIdxL = []
  220. yV = rmsV[pkIdxL]
  221. refPkDb = yV[min_pk_idx]
  222. for i in range( min_pk_idx+1, max_pk_idx+1 ):
  223. if yV[i] > refPkDb:
  224. refPkDb = yV[i]
  225. else:
  226. skipPkIdxL.append(i)
  227. return skipPkIdxL
  228. def find_out_of_range_peaks( rmsV, pkIdxL, min_pk_idx, max_pk_idx, maxDeltaDb ):
  229. """ Locate peaks which are more than maxDeltaDb from the previous peak.
  230. If two peaks are separated by more than maxDeltaDb then the range must be resampled
  231. """
  232. oorPkIdxL = []
  233. yV = rmsV[pkIdxL]
  234. for i in range( min_pk_idx, max_pk_idx+1 ):
  235. if i > 0:
  236. d = yV[i] - yV[i-1]
  237. if d > maxDeltaDb or d < 0:
  238. oorPkIdxL.append(i)
  239. return oorPkIdxL
  240. def plot_spectrum( ax, srate, binHz, specV, midiPitch, harmN ):
  241. """ Plot a single spectrum, 'specV' and the harmonic peak location boundaries."""
  242. binN = specV.shape[0]
  243. harmLBinL,harmMBinL,harmUBinL = calc_harm_bins( srate, binHz, midiPitch, harmN )
  244. fundHz = harmMBinL[0] * binHz
  245. maxPlotHz = fundHz * (harmN+1)
  246. maxPlotBinN = int(round(maxPlotHz/binHz))
  247. hzV = np.arange(binN) * (srate/(binN*2))
  248. specV = 20.0 * np.log10(specV)
  249. ax.plot(hzV[0:maxPlotBinN], specV[0:maxPlotBinN] )
  250. for h0,h1,h2 in zip(harmLBinL,harmMBinL,harmUBinL):
  251. ax.axvline( x=h0 * binHz, color="blue")
  252. ax.axvline( x=h1 * binHz, color="black")
  253. ax.axvline( x=h2 * binHz, color="blue")
  254. ax.set_ylabel(str(midiPitch))
  255. def plot_spectral_ranges( inDir, pitchL, rmsWndMs=300, rmsHopMs=30, harmN=5, dbLinRef=0.001 ):
  256. """ Plot the spectrum from one note (7th from last) in each attack pulse length sequence referred to by pitchL."""
  257. plotN = len(pitchL)
  258. fig,axL = plt.subplots(plotN,1)
  259. for plot_idx,midiPitch in enumerate(pitchL):
  260. # get the audio and meta-data file names
  261. seqFn = os.path.join( inDir, str(midiPitch), "seq.json")
  262. audioFn = os.path.join( inDir, str(midiPitch), "audio.wav")
  263. # read the meta data object
  264. with open( seqFn, "rb") as f:
  265. r = json.load(f)
  266. # read the audio file
  267. srate, signalM = wavfile.read(audioFn)
  268. sigV = signalM / float(0x7fff)
  269. # calc. the RMS envelope in the time domain
  270. rms0DbV, rms0_srate = audio_rms( srate, sigV, rmsWndMs, rmsHopMs, dbLinRef )
  271. # locate the sample index of the peak of each note attack
  272. pkIdx0L = locate_peak_indexes( rms0DbV, rms0_srate, r['eventTimeL'] )
  273. # select the 7th to last note for spectrum measurement
  274. #
  275. # TODO: come up with a better way to select the note to measure
  276. #
  277. spectrumSmpIdx = pkIdx0L[ len(pkIdx0L) - 7 ]
  278. # calc. the RMS envelope by taking the max spectral peak in each STFT window
  279. rmsDbV, rms_srate, specV, specHopIdx, binHz = audio_stft_rms( srate, sigV, rmsWndMs, rmsHopMs, dbLinRef, spectrumSmpIdx)
  280. # specV[] is the spectrum of the note at spectrumSmpIdx
  281. # plot the spectrum and the harmonic selection ranges
  282. plot_spectrum( axL[plot_idx], srate, binHz, specV, midiPitch, harmN )
  283. plt.show()
  284. def td_plot( ax, inDir, midi_pitch, id, analysisArgsD ):
  285. r = rms_analysis_main( inDir, midi_pitch, **analysisArgsD['rmsAnalysisArgs'] )
  286. min_pk_idx, max_pk_idx = find_min_max_peak_index( r.pkDbL, analysisArgsD['minAttkDb'], analysisArgsD['maxDbOffset'] )
  287. skipPkIdxL = find_skip_peaks( r.rmsDbV, r.pkIdxL, min_pk_idx, max_pk_idx )
  288. jmpPkIdxL = find_out_of_range_peaks( r.rmsDbV, r.pkIdxL, min_pk_idx, max_pk_idx, analysisArgsD['maxDeltaDb'] )
  289. secV = np.arange(0,len(r.rmsDbV)) / r.rms_srate
  290. ax.plot( secV, r.rmsDbV )
  291. ax.plot( np.arange(0,len(r.tdRmsDbV)) / r.rms_srate, r.tdRmsDbV, color="black" )
  292. # print beg/end boundaries
  293. for i,(begMs, endMs) in enumerate(r.eventTimeL):
  294. pkSec = r.pkIdxL[i] / r.rms_srate
  295. endSec = pkSec + r.statsL[i].durMs / 1000.0
  296. ax.axvline( x=begMs/1000.0, color="green")
  297. ax.axvline( x=endMs/1000.0, color="red")
  298. ax.axvline( x=pkSec, color="black")
  299. ax.axvline( x=endSec, color="black")
  300. ax.text(begMs/1000.0, 20.0, str(i) )
  301. # plot peak markers
  302. for i,pki in enumerate(r.pkIdxL):
  303. marker = 4 if i==min_pk_idx or i==max_pk_idx else 5
  304. color = "red" if i in skipPkIdxL else "black"
  305. ax.plot( [pki / r.rms_srate], [ r.rmsDbV[pki] ], marker=marker, color=color)
  306. if i in jmpPkIdxL:
  307. ax.plot( [pki / r.rms_srate], [ r.rmsDbV[pki] ], marker=6, color="blue")
  308. return r
  309. def do_td_plot( inDir, analysisArgs ):
  310. fig,axL = plt.subplots(3,1)
  311. fig.set_size_inches(18.5, 10.5, forward=True)
  312. id = int(inDir.split("/")[-1])
  313. midi_pitch = int(inDir.split("/")[-2])
  314. r = td_plot(axL[0],inDir,midi_pitch,id,analysisArgs)
  315. qualityV = np.array([ x.quality for x in r.statsL ]) * np.max(r.pkDbL)
  316. durMsV = np.array([ x.durMs for x in r.statsL ])
  317. avgV = np.array([ x.durAvgDb for x in r.statsL ])
  318. #durMsV[ durMsV < 400 ] = 0
  319. #durMsV = durMsV * np.max(r.pkDbL)/np.max(durMsV)
  320. #durMsV = durMsV / 100.0
  321. dV = np.diff(r.pkDbL) / r.pkDbL[1:]
  322. axL[1].plot( r.pkUsL, r.pkDbL, marker='.',label="pkDb" )
  323. axL[1].plot( r.pkUsL, qualityV, marker='.',label="quality" )
  324. axL[1].plot( r.pkUsL, avgV, marker='.',label="avgDb" )
  325. #axL[2].plot( r.pkUsL, durMsV, marker='.' )
  326. axL[2].plot( r.pkUsL[1:], dV, marker='.',label='d')
  327. axL[2].set_ylim([-1,1])
  328. axL[1].legend()
  329. sni = select_first_stable_note_by_dur( durMsV )
  330. if sni is not None:
  331. axL[1].plot( r.pkUsL[sni], r.pkDbL[sni], marker='*', color='red')
  332. sni = select_first_stable_note_by_delta_db( r.pkDbL )
  333. if sni is not None:
  334. axL[2].plot( r.pkUsL[sni], dV[sni-1], marker='*', color='red')
  335. for i,s in enumerate(r.statsL):
  336. axL[1].text( r.pkUsL[i], r.pkDbL[i] + 1, "%i" % (i))
  337. for i in range(1,len(r.pkUsL)):
  338. axL[2].text( r.pkUsL[i], dV[i-1], "%i" % (i))
  339. plt.show()
  340. def do_td_multi_plot( inDir, analysisArgs ):
  341. midi_pitch = int(inDir.split("/")[-1])
  342. dirL = os.listdir(inDir)
  343. fig,axL = plt.subplots(len(dirL),1)
  344. for id,(idir,ax) in enumerate(zip(dirL,axL)):
  345. td_plot(ax, os.path.join(inDir,str(id)), midi_pitch, id, analysisArgs )
  346. plt.show()
  347. if __name__ == "__main__":
  348. inDir = sys.argv[1]
  349. cfgFn = sys.argv[2]
  350. take_id = None if len(sys.argv)<4 else sys.argv[3]
  351. cfg = parse_yaml_cfg( cfgFn )
  352. if take_id is not None:
  353. inDir = os.path.join(inDir,take_id)
  354. do_td_plot(inDir,cfg.analysisArgs)
  355. else:
  356. #do_td_multi_plot(inDir,cfg.analysisArgs)
  357. #plot_spectral_ranges( inDir, [ 24, 36, 48, 60, 72, 84, 96, 104] )
  358. plot_resample_pulse_times( inDir, cfg.analysisArgs )