picadae calibration programs
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

plot_seq_1.py 35KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088
  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,json,math
  4. import matplotlib.pyplot as plt
  5. import numpy as np
  6. from common import parse_yaml_cfg
  7. import rms_analysis
  8. import elbow
  9. def fit_to_reference( pkL, refTakeId ):
  10. us_outL = []
  11. db_outL = []
  12. dur_outL = []
  13. tid_outL = []
  14. dbL,usL,durMsL,takeIdL = tuple(zip(*pkL))
  15. us_refL,db_refL,dur_refL = zip(*[(usL[i],dbL[i],durMsL[i]) for i in range(len(usL)) if takeIdL[i]==refTakeId])
  16. for takeId in set(takeIdL):
  17. us0L,db0L,dur0L = zip(*[(usL[i],dbL[i],durMsL[i]) for i in range(len(usL)) if takeIdL[i]==takeId ])
  18. if takeId == refTakeId:
  19. db_outL += db0L
  20. else:
  21. db1V = elbow.fit_points_to_reference(us0L,db0L,us_refL,db_refL)
  22. if db1V is not None:
  23. db_outL += db1V.tolist()
  24. us_outL += us0L
  25. dur_outL+= dur0L
  26. tid_outL+= [takeId] * len(us0L)
  27. return zip(db_outL,us_outL,dur_outL,tid_outL)
  28. def get_merged_pulse_db_measurements( inDir, midi_pitch, analysisArgsD, takeId=None ):
  29. inDir = os.path.join(inDir,"%i" % (midi_pitch))
  30. takeIdL = []
  31. if takeId is not None:
  32. takeIdL.append(takeId)
  33. else:
  34. takeDirL = os.listdir(inDir)
  35. # for each take in this directory
  36. for take_folder in takeDirL:
  37. takeIdL.append(int(take_folder))
  38. pkL = []
  39. refTakeId = None
  40. usRefL = None
  41. dbRefL = None
  42. # for each take in this directory
  43. for take_number in takeIdL:
  44. if refTakeId is None:
  45. refTakeId = take_number
  46. # analyze this takes audio and locate the note peaks
  47. r = rms_analysis.rms_analysis_main( os.path.join(inDir,str(take_number)), midi_pitch, **analysisArgsD )
  48. # store the peaks in pkL[ (db,us) ]
  49. for db,us,stats in zip(r.pkDbL,r.pkUsL,r.statsL):
  50. pkL.append( (db,us,stats.durMs,take_number) )
  51. pkUsL = []
  52. pkDbL = []
  53. durMsL = []
  54. takeIdL = []
  55. holdDutyPctL = []
  56. if refTakeId is None:
  57. print("No valid data files at %s pitch:%i" % (inDir,midi_pitch))
  58. else:
  59. pkL = fit_to_reference( pkL, refTakeId )
  60. # sort the peaks on increasing attack pulse microseconds
  61. pkL = sorted( pkL, key= lambda x: x[1] )
  62. # merge sample points that separated by less than 'minSampleDistUs' milliseconds
  63. #pkL = merge_close_sample_points( pkL, analysisArgsD['minSampleDistUs'] )
  64. # split pkL
  65. pkDbL,pkUsL,durMsL,takeIdL = tuple(zip(*pkL))
  66. return pkUsL,pkDbL,durMsL,takeIdL,r.holdDutyPctL
  67. def select_resample_reference_indexes( noiseIdxL ):
  68. resampleIdxS = set()
  69. # for each noisy sample index store that index and the index
  70. # before and after it
  71. for i in noiseIdxL:
  72. resampleIdxS.add( i )
  73. if i+1 < len(noiseIdxL):
  74. resampleIdxS.add( i+1 )
  75. if i-1 >= 0:
  76. resampleIdxS.add( i-1 )
  77. resampleIdxL = list(resampleIdxS)
  78. # if a single sample point is left out of a region of
  79. # contiguous sample points then include this as a resample point also
  80. for i in resampleIdxL:
  81. if i + 1 not in resampleIdxL and i + 2 in resampleIdxL: # BUG BUG BUG: Hardcoded constant
  82. if i+1 < len(noiseIdxL):
  83. resampleIdxL.append(i+1)
  84. return resampleIdxL
  85. def locate_resample_regions( usL, dbL, resampleIdxL ):
  86. # locate regions of points to resample
  87. regionL = [] # (bi,ei)
  88. inRegionFl = False
  89. bi = None
  90. for i in range(len(usL)):
  91. if inRegionFl:
  92. if i not in resampleIdxL:
  93. regionL.append((bi,i-1))
  94. inRegionFl = False
  95. bi = None
  96. else:
  97. if i in resampleIdxL:
  98. inRegionFl = True
  99. bi = i
  100. if bi is not None:
  101. regionL.append((bi,len(usL)-1))
  102. # select points around and within the resample regions
  103. # to resample
  104. reUsL = []
  105. reDbL = []
  106. for bi,ei in regionL:
  107. for i in range(bi,ei+2):
  108. if i == 0:
  109. us = usL[i]
  110. db = dbL[i]
  111. elif i >= len(usL):
  112. us = usL[i-1]
  113. db = dbL[i-1]
  114. else:
  115. us = usL[i-1] + (usL[i]-usL[i-1])/2
  116. db = dbL[i-1] + (dbL[i]-dbL[i-1])/2
  117. reUsL.append(us)
  118. reDbL.append(db)
  119. return reUsL,reDbL
  120. def get_dur_skip_indexes( durMsL, dbL, takeIdL, scoreL, minDurMs, minDb, noiseLimitPct ):
  121. firstAudibleIdx = None
  122. firstNonSkipIdx = None
  123. # get the indexes of samples which do not meet the duration, db level, or noise criteria
  124. skipIdxL = [ i for i,(ms,db,score) in enumerate(zip(durMsL,dbL,scoreL)) if ms < minDurMs or db < minDb or score > noiseLimitPct ]
  125. # if a single sample point is left out of a region of
  126. # contiguous skipped points then skip this point also
  127. for i in range(len(durMsL)):
  128. if i not in skipIdxL and i-1 in skipIdxL and i+1 in skipIdxL:
  129. skipIdxL.append(i)
  130. # find the first set of 3 contiguous samples that
  131. # are greater than minDurMs - all samples prior
  132. # to these will be skipped
  133. xL = []
  134. for i in range(len(durMsL)):
  135. if i in skipIdxL:
  136. xL = []
  137. else:
  138. xL.append(i)
  139. if len(xL) == 3: # BUG BUG BUG: Hardcoded constant
  140. firstAudibleIdx = xL[0]
  141. break
  142. # decrease by one decibel to locate the first non-skip
  143. # TODO: what if no note exists that is one decibel less
  144. # The recordings of very quiet notes do not give reliabel decibel measures
  145. # so this may not be the best backup criteria
  146. if firstAudibleIdx is not None:
  147. i = firstAudibleIdx-1
  148. while abs(dbL[i] - dbL[firstAudibleIdx]) < 1.0: # BUG BUG BUG: Hardcoded constant
  149. i -= 1
  150. firstNonSkipIdx = i
  151. return skipIdxL, firstAudibleIdx, firstNonSkipIdx
  152. def get_resample_points( usL, dbL, durMsL, takeIdL, minDurMs, minDb, noiseLimitPct ):
  153. scoreV = np.abs( rms_analysis.samples_to_linear_residual( usL, dbL) * 100.0 / dbL )
  154. skipIdxL, firstAudibleIdx, firstNonSkipIdx = get_dur_skip_indexes( durMsL, dbL, takeIdL, scoreV.tolist(), minDurMs, minDb, noiseLimitPct )
  155. skipL = [ (usL[i],dbL[i]) for i in skipIdxL ]
  156. noiseIdxL = [ i for i in range(scoreV.shape[0]) if scoreV[i] > noiseLimitPct ]
  157. noiseL = [ (usL[i],dbL[i]) for i in noiseIdxL ]
  158. resampleIdxL = select_resample_reference_indexes( noiseIdxL )
  159. if firstNonSkipIdx is not None:
  160. resampleIdxL = [ i for i in resampleIdxL if i >= firstNonSkipIdx ]
  161. resampleL = [ (usL[i],dbL[i]) for i in resampleIdxL ]
  162. reUsL,reDbL = locate_resample_regions( usL, dbL, resampleIdxL )
  163. return reUsL, reDbL, noiseL, resampleL, skipL, firstAudibleIdx, firstNonSkipIdx
  164. def get_resample_points_wrap( inDir, midi_pitch, analysisArgsD ):
  165. usL, dbL, durMsL,_,_ = get_merged_pulse_db_measurements( inDir, midi_pitch, analysisArgsD['rmsAnalysisArgs'] )
  166. reUsL,_,_,_,_,_,_ = get_resample_points( usL, dbL, durMsL, analysisArgsD['resampleMinDurMs'], analysisArgsD['resampleMinDb'], analysisArgsD['resampleNoiseLimitPct'] )
  167. return reUsL
  168. def plot_us_db_curves( ax, inDir, keyMapD, midi_pitch, analysisArgsD, plotResamplePointsFl=False, plotTakesFl=True, usMax=None ):
  169. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, analysisArgsD['rmsAnalysisArgs'] )
  170. reUsL, reDbL, noiseL, resampleL, skipL, firstAudibleIdx, firstNonSkipIdx = get_resample_points( usL, dbL, durMsL, takeIdL, analysisArgsD['resampleMinDurMs'], analysisArgsD['resampleMinDb'], analysisArgsD['resampleNoiseLimitPct'] )
  171. # plot first audible and non-skip position
  172. if False:
  173. if firstNonSkipIdx is not None:
  174. ax.plot( usL[firstNonSkipIdx], dbL[firstNonSkipIdx], markersize=15, marker='+', linestyle='None', color='red')
  175. if firstAudibleIdx is not None:
  176. ax.plot( usL[firstAudibleIdx], dbL[firstAudibleIdx], markersize=15, marker='*', linestyle='None', color='red')
  177. # plot the resample points
  178. if plotResamplePointsFl:
  179. ax.plot( reUsL, reDbL, markersize=13, marker='x', linestyle='None', color='green')
  180. # plot the noisy sample positions
  181. if noiseL:
  182. nUsL,nDbL = zip(*noiseL)
  183. ax.plot( nUsL, nDbL, marker='o', markersize=9, linestyle='None', color='black')
  184. # plot the noisy sample positions and the neighbors included in the noisy region
  185. if resampleL:
  186. nUsL,nDbL = zip(*resampleL)
  187. ax.plot( nUsL, nDbL, marker='+', markersize=8, linestyle='None', color='red')
  188. # plot actual sample points
  189. elbow_us = None
  190. elbow_db = None
  191. elbow_len = None
  192. usL,dbL,takeIdL = zip(*[(us,dbL[i],takeIdL[i]) for i,us in enumerate(usL) if usMax is None or us <= usMax])
  193. if plotTakesFl:
  194. for takeId in list(set(takeIdL)):
  195. # get the us,db points included in this take
  196. xL,yL = zip(*[(usL[i],dbL[i]) for i in range(len(usL)) if takeIdL[i]==takeId ])
  197. ax.plot(xL,yL, marker='.',label=takeId)
  198. for i,(x,y) in enumerate(zip(xL,yL)):
  199. ax.text(x,y,str(i))
  200. #if elbow_len is None or len(xL) > elbow_len:
  201. if takeId+1 == len(set(takeIdL)):
  202. elbow_us,elbow_db = elbow.find_elbow(xL,yL)
  203. elbow_len = len(xL)
  204. else:
  205. ax.plot(usL, dbL, marker='.')
  206. ax.plot([elbow_us],[elbow_db],marker='*',markersize=12,color='red',linestyle='None')
  207. # plot the skip points in yellow
  208. if False:
  209. if skipL:
  210. nUsL,nDbL = zip(*skipL)
  211. ax.plot( nUsL, nDbL, marker='.', linestyle='None', color='yellow')
  212. # plot the locations where the hold duty cycle changes with vertical black lines
  213. for us_duty in holdDutyPctL:
  214. us,duty = tuple(us_duty)
  215. if us > 0:
  216. ax.axvline(us,color='black')
  217. # plot the 'minDb' reference line
  218. ax.axhline(analysisArgsD['resampleMinDb'] ,color='black')
  219. if os.path.isfile("minInterpDb.json"):
  220. with open("minInterpDb.json","r") as f:
  221. r = json.load(f)
  222. if midi_pitch in r['pitchL']:
  223. ax.axhline( r['minDbL'][ r['pitchL'].index(midi_pitch) ], color='blue' )
  224. ax.axhline( r['maxDbL'][ r['pitchL'].index(midi_pitch) ], color='blue' )
  225. ax.set_ylabel( "%i %s %s" % (midi_pitch, keyMapD[midi_pitch]['type'],keyMapD[midi_pitch]['class']))
  226. def plot_us_db_take_curve( ax, inDir, keyMapD, midi_pitch, takeId, analysisArgsD ):
  227. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, analysisArgsD['rmsAnalysisArgs'] )
  228. reUsL, reDbL, noiseL, resampleL, skipL, firstAudibleIdx, firstNonSkipIdx = get_resample_points( usL, dbL, durMsL, takeIdL, analysisArgsD['resampleMinDurMs'], analysisArgsD['resampleMinDb'], analysisArgsD['resampleNoiseLimitPct'] )
  229. # plot first audible and non-skip position
  230. if False:
  231. if firstNonSkipIdx is not None:
  232. ax.plot( usL[firstNonSkipIdx], dbL[firstNonSkipIdx], markersize=15, marker='+', linestyle='None', color='red')
  233. if firstAudibleIdx is not None:
  234. ax.plot( usL[firstAudibleIdx], dbL[firstAudibleIdx], markersize=15, marker='*', linestyle='None', color='red')
  235. # plot the resample points
  236. if plotResamplePointsFl:
  237. ax.plot( reUsL, reDbL, markersize=13, marker='x', linestyle='None', color='green')
  238. # plot the noisy sample positions
  239. if noiseL:
  240. nUsL,nDbL = zip(*noiseL)
  241. ax.plot( nUsL, nDbL, marker='o', markersize=9, linestyle='None', color='black')
  242. # plot the noisy sample positions and the neighbors included in the noisy region
  243. if resampleL:
  244. nUsL,nDbL = zip(*resampleL)
  245. ax.plot( nUsL, nDbL, marker='+', markersize=8, linestyle='None', color='red')
  246. # plot actual sample points
  247. elbow_us = None
  248. elbow_db = None
  249. elbow_len = None
  250. usL,dbL,takeIdL = zip(*[(us,dbL[i],takeIdL[i]) for i,us in enumerate(usL) if usMax is None or us <= usMax])
  251. if plotTakesFl:
  252. for takeId in list(set(takeIdL)):
  253. # get the us,db points included in this take
  254. xL,yL = zip(*[(usL[i],dbL[i]) for i in range(len(usL)) if takeIdL[i]==takeId ])
  255. ax.plot(xL,yL, marker='.',label=takeId)
  256. for i,(x,y) in enumerate(zip(xL,yL)):
  257. ax.text(x,y,str(i))
  258. #if elbow_len is None or len(xL) > elbow_len:
  259. if takeId+1 == len(set(takeIdL)):
  260. elbow_us,elbow_db = elbow.find_elbow(xL,yL)
  261. elbow_len = len(xL)
  262. else:
  263. ax.plot(usL, dbL, marker='.')
  264. ax.plot([elbow_us],[elbow_db],marker='*',markersize=12,color='red',linestyle='None')
  265. # plot the skip points in yellow
  266. if False:
  267. if skipL:
  268. nUsL,nDbL = zip(*skipL)
  269. ax.plot( nUsL, nDbL, marker='.', linestyle='None', color='yellow')
  270. # plot the locations where the hold duty cycle changes with vertical black lines
  271. for us_duty in holdDutyPctL:
  272. us,duty = tuple(us_duty)
  273. if us > 0:
  274. ax.axvline(us,color='black')
  275. # plot the 'minDb' reference line
  276. ax.axhline(analysisArgsD['resampleMinDb'] ,color='black')
  277. if os.path.isfile("minInterpDb.json"):
  278. with open("minInterpDb.json","r") as f:
  279. r = json.load(f)
  280. if midi_pitch in r['pitchL']:
  281. ax.axhline( r['minDbL'][ r['pitchL'].index(midi_pitch) ], color='blue' )
  282. ax.axhline( r['maxDbL'][ r['pitchL'].index(midi_pitch) ], color='blue' )
  283. ax.set_ylabel( "%i %s %s" % (midi_pitch, keyMapD[midi_pitch]['type'],keyMapD[midi_pitch]['class']))
  284. def plot_us_db_curves_main( inDir, cfg, pitchL, plotTakesFl=True, usMax=None, printDir="" ):
  285. analysisArgsD = cfg.analysisArgs
  286. keyMapD = { d['midi']:d for d in cfg.key_mapL }
  287. axN = len(pitchL)
  288. fig,axL = plt.subplots(axN,1,sharex=True)
  289. if axN == 1:
  290. axL = [axL]
  291. fig.set_size_inches(18.5, 10.5*axN)
  292. for ax,midi_pitch in zip(axL,pitchL):
  293. plot_us_db_curves( ax,inDir, keyMapD, midi_pitch, analysisArgsD, plotTakesFl=plotTakesFl, usMax=usMax )
  294. if plotTakesFl:
  295. plt.legend()
  296. if printDir:
  297. plt.savefig(os.path.join(printDir,"us_db.png"),format="png")
  298. plt.show()
  299. def _plot_us_db_takes( inDir, cfg, pitchL, takeIdL, printDir="", printFn="" ):
  300. assert( len(pitchL) == len(takeIdL) )
  301. analysisArgsD = cfg.analysisArgs
  302. keyMapD = { d['midi']:d for d in cfg.key_mapL }
  303. fig,ax = plt.subplots(1,1)
  304. fig.set_size_inches(18.5, 10.5)
  305. for midi_pitch,takeId in zip(pitchL,takeIdL):
  306. usL, dbL, durMsL, _, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, analysisArgsD['rmsAnalysisArgs'], takeId=takeId )
  307. ax.plot(usL,dbL, marker='.',label="%i:%i %s %s" % (midi_pitch,takeId,keyMapD[midi_pitch]['class'],keyMapD[midi_pitch]['type']))
  308. for i,(x,y) in enumerate(zip(usL,dbL)):
  309. ax.text(x,y,str(i))
  310. f_usL,f_dbL = filter_us_db(usL,dbL)
  311. ax.plot(f_usL,f_dbL, marker='.')
  312. elbow_us,elbow_db = elbow.find_elbow(usL,dbL)
  313. ax.plot([elbow_us],[elbow_db],marker='*',markersize=12,color='red',linestyle='None')
  314. elb_idx = nearest_sample_point( dbL, usL, elbow_db, elbow_us )
  315. if printDir:
  316. plt.savefig(os.path.join(printDir,printFn),format="png")
  317. plt.legend()
  318. plt.show()
  319. def get_pitches_and_takes( inDir ):
  320. pitchD = {}
  321. inDirL = os.listdir( inDir )
  322. for pitch in inDirL:
  323. path = os.path.join( inDir, pitch )
  324. takeIdL = os.listdir( path )
  325. takeIdL = sorted([ int(takeId) for takeId in takeIdL ])
  326. takeIdL = [ str(x) for x in takeIdL ]
  327. pitchD[int(pitch)] = takeIdL
  328. return pitchD
  329. def plot_us_db_takes( inDir, cfg, pitchL, printDir=""):
  330. takeIdL = None
  331. takeIdL = [ pitchL[i] for i in range(1,len(pitchL),2) ]
  332. pitchL = [ pitchL[i] for i in range(0,len(pitchL),2) ]
  333. return _plot_us_db_takes( inDir, cfg, pitchL, takeIdL, printDir, "us_db_takes.png")
  334. def plot_us_db_takes_last( inDir, cfg, pitchL, printDir ):
  335. pitchD = get_pitches_and_takes(inDir)
  336. takeIdL = []
  337. for pitch in pitchL:
  338. takeIdL.append( int(pitchD[pitch][-1]) )
  339. return _plot_us_db_takes( inDir, cfg, pitchL, takeIdL, printDir, "us_db_takes_last.png")
  340. def plot_all_noise_curves( inDir, cfg, pitchL=None ):
  341. pitchFolderL = os.listdir(inDir)
  342. if pitchL is None:
  343. pitchL = [ int( int(pitchFolder) ) for pitchFolder in pitchFolderL ]
  344. fig,ax = plt.subplots()
  345. for midi_pitch in pitchL:
  346. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, cfg.analysisArgs['rmsAnalysisArgs'] )
  347. scoreV = np.abs( rms_analysis.samples_to_linear_residual( usL, dbL) * 100.0 / dbL )
  348. minDurMs = cfg.analysisArgs['resampleMinDurMs']
  349. minDb = cfg.analysisArgs['resampleMinDb']
  350. noiseLimitPct = cfg.analysisArgs['resampleNoiseLimitPct']
  351. skipIdxL, firstAudibleIdx, firstNonSkipIdx = get_dur_skip_indexes( durMsL, dbL, scoreV.tolist(), takeIdL, minDurMs, minDb, noiseLimitPct )
  352. if False:
  353. ax.plot( usL[firstAudibleIdx], scoreV[firstAudibleIdx], markersize=10, marker='*', linestyle='None', color='red')
  354. ax.plot( usL, scoreV, label="%i"%(midi_pitch) )
  355. ax.set_xlabel('us')
  356. else:
  357. xL = [ (score,db,i) for i,(score,db) in enumerate(zip(scoreV,dbL)) ]
  358. xL = sorted(xL, key=lambda x: x[1] )
  359. scoreV,dbL,idxL = zip(*xL)
  360. ax.plot( dbL[idxL[firstAudibleIdx]], scoreV[idxL[firstAudibleIdx]], markersize=10, marker='*', linestyle='None', color='red')
  361. ax.plot( dbL, scoreV, label="%i"%(midi_pitch) )
  362. ax.set_xlabel('db')
  363. ax.set_ylabel("noise db %")
  364. plt.legend()
  365. plt.show()
  366. def nearest_sample_point( dbL, usL, db0, us0 ):
  367. xL = np.array([ abs(us-us0) for db,us in zip(dbL,usL) ])
  368. return np.argmin(xL)
  369. def plot_min_max_2_db( inDir, cfg, pitchL=None, takeId=2, printDir=None ):
  370. takeIdArg = takeId
  371. pitchTakeD = get_pitches_and_takes(inDir)
  372. pitchFolderL = os.listdir(inDir)
  373. if pitchL is None:
  374. pitchL = [ int( int(pitchFolder) ) for pitchFolder in pitchFolderL ]
  375. okL = []
  376. outPitchL = []
  377. minDbL = []
  378. maxDbL = []
  379. for midi_pitch in pitchL:
  380. takeId = None
  381. if takeIdArg == -1:
  382. takeId = pitchTakeD[midi_pitch][-1]
  383. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, cfg.analysisArgs['rmsAnalysisArgs'], takeId )
  384. okL.append(False)
  385. db_maxL = sorted(dbL)
  386. maxDbL.append( np.mean(db_maxL[-5:]) )
  387. usL,dbL = zip(*[(usL[i],dbL[i]) for i in range(len(usL)) if takeIdL[i]==takeId ])
  388. if len(set(takeIdL)) == 3:
  389. okL[-1] = True
  390. elbow_us,elbow_db = elbow.find_elbow(usL,dbL)
  391. minDbL.append(elbow_db)
  392. outPitchL.append(midi_pitch)
  393. smp_idx = nearest_sample_point( dbL, usL, elbow_db, elbow_us )
  394. print(" %i:[-1,%i], " % (midi_pitch,smp_idx))
  395. p_dL = sorted( zip(outPitchL,minDbL,maxDbL,okL), key=lambda x: x[0] )
  396. outPitchL,minDbL,maxDbL,okL = zip(*p_dL)
  397. fig,ax = plt.subplots()
  398. ax.plot(outPitchL,minDbL)
  399. ax.plot(outPitchL,maxDbL)
  400. keyMapD = { d['midi']:d for d in cfg.key_mapL }
  401. for pitch,min_db,max_db,okFl in zip(outPitchL,minDbL,maxDbL,okL):
  402. c = 'black' if okFl else 'red'
  403. ax.text( pitch, min_db, "%i %s %s" % (pitch, keyMapD[pitch]['type'],keyMapD[pitch]['class']), color=c)
  404. ax.text( pitch, max_db, "%i %s %s" % (pitch, keyMapD[pitch]['type'],keyMapD[pitch]['class']), color=c)
  405. if printDir:
  406. plt.savefig(os.path.join(printDir,"min_max_db_2.png"),format="png")
  407. plt.show()
  408. def plot_min_db_manual( inDir, cfg, printDir=None, absMaxDb=27, absMinDb=3 ):
  409. pitchTakeD = get_pitches_and_takes(inDir)
  410. pitchL = list(cfg.manualMinD.keys())
  411. outPitchL = []
  412. maxDbL = []
  413. minDbL = []
  414. okL = []
  415. anchorMinDbL = []
  416. anchorMaxDbL = []
  417. for midi_pitch in pitchL:
  418. if cfg.manualLastFl:
  419. manual_take_id = pitchTakeD[midi_pitch][-1]
  420. takeId = manual_take_id
  421. else:
  422. manual_take_id = cfg.manualMinD[midi_pitch][0]
  423. takeId = None
  424. manual_sample_idx = cfg.manualMinD[midi_pitch][1]
  425. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, cfg.analysisArgs['rmsAnalysisArgs'], takeId )
  426. okL.append(False)
  427. if takeId is None:
  428. takeId = len(set(takeIdL))-1
  429. # most pitches have 3 sample takes that do not
  430. if len(set(takeIdL)) == 3 and manual_take_id == takeId:
  431. okL[-1] = True
  432. else:
  433. okL[-1] = True
  434. # maxDb is computed on all takes (not just the specified take)
  435. db_maxL = sorted(dbL)
  436. max_db = min(absMaxDb,np.mean(db_maxL[-4:]))
  437. maxDbL.append( max_db )
  438. # get the us,db values for the specified take
  439. usL,dbL = zip(*[(usL[i],dbL[i]) for i in range(len(usL)) if takeIdL[i]==manual_take_id ])
  440. # min db from the sample index manually specified in cfg
  441. manualMinDb = max(absMinDb,dbL[ manual_sample_idx ])
  442. minDbL.append( manualMinDb )
  443. outPitchL.append(midi_pitch)
  444. if midi_pitch in cfg.manualAnchorPitchMinDbL:
  445. anchorMinDbL.append( manualMinDb )
  446. if midi_pitch in cfg.manualAnchorPitchMaxDbL:
  447. anchorMaxDbL.append( max_db )
  448. # Form the complete set of min/max db levels for each pitch by interpolating the
  449. # db values between the manually selected anchor points.
  450. interpMinDbL = np.interp( pitchL, cfg.manualAnchorPitchMinDbL, anchorMinDbL )
  451. interpMaxDbL = np.interp( pitchL, cfg.manualAnchorPitchMaxDbL, anchorMaxDbL )
  452. fig,ax = plt.subplots()
  453. ax.plot(outPitchL,minDbL) # plot the manually selected minDb values
  454. ax.plot(outPitchL,maxDbL) # plot the max db values
  455. # plot the interpolated minDb/maxDb values
  456. ax.plot(pitchL,interpMinDbL)
  457. ax.plot(pitchL,interpMaxDbL)
  458. keyMapD = { d['midi']:d for d in cfg.key_mapL }
  459. for pitch,min_db,max_db,okFl in zip(outPitchL,minDbL,maxDbL,okL):
  460. c = 'black' if okFl else 'red'
  461. ax.text( pitch, min_db, "%i %s %s" % (pitch, keyMapD[pitch]['type'],keyMapD[pitch]['class']), color=c)
  462. ax.text( pitch, max_db, "%i %s %s" % (pitch, keyMapD[pitch]['type'],keyMapD[pitch]['class']), color=c)
  463. with open("minInterpDb.json",'w') as f:
  464. json.dump( { "pitchL":pitchL, "minDbL":list(interpMinDbL), "maxDbL":list(interpMaxDbL) }, f )
  465. if printDir:
  466. plt.savefig(os.path.join(printDir,"manual_db.png"),format="png")
  467. plt.show()
  468. def plot_min_max_db( inDir, cfg, pitchL=None, printDir=None ):
  469. pitchFolderL = os.listdir(inDir)
  470. if pitchL is None:
  471. pitchL = [ int( int(pitchFolder) ) for pitchFolder in pitchFolderL ]
  472. maxDbL = []
  473. minDbL = []
  474. for midi_pitch in pitchL:
  475. print(midi_pitch)
  476. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, midi_pitch, cfg.analysisArgs['rmsAnalysisArgs'] )
  477. scoreV = np.abs( rms_analysis.samples_to_linear_residual( usL, dbL) * 100.0 / dbL )
  478. minDurMs = cfg.analysisArgs['resampleMinDurMs']
  479. minDb = cfg.analysisArgs['resampleMinDb']
  480. noiseLimitPct = cfg.analysisArgs['resampleNoiseLimitPct']
  481. skipIdxL, firstAudibleIdx, firstNonSkipIdx = get_dur_skip_indexes( durMsL, dbL, takeIdL, scoreV.tolist(), minDurMs, minDb, noiseLimitPct )
  482. minDbL.append( dbL[firstAudibleIdx] )
  483. dbL = sorted(dbL)
  484. x = np.mean(dbL[-3:])
  485. x = np.max(dbL)
  486. maxDbL.append( x )
  487. fig,ax = plt.subplots()
  488. fig.set_size_inches(18.5, 10.5)
  489. p_dL = sorted( zip(pitchL,maxDbL), key=lambda x: x[0] )
  490. pitchL,maxDbL = zip(*p_dL)
  491. ax.plot(pitchL,maxDbL)
  492. ax.plot(pitchL,minDbL)
  493. for pitch,db in zip(pitchL,maxDbL):
  494. keyMapD = { d['midi']:d for d in cfg.key_mapL }
  495. ax.text( pitch, db, "%i %s %s" % (pitch, keyMapD[pitch]['type'],keyMapD[pitch]['class']))
  496. if printDir:
  497. plt.savefig(os.path.join(printDir,"min_max_db.png"),format="png")
  498. plt.show()
  499. def estimate_us_to_db_map( inDir, cfg, minMapDb=16.0, maxMapDb=26.0, incrMapDb=0.5, pitchL=None ):
  500. pitchFolderL = os.listdir(inDir)
  501. if pitchL is None:
  502. pitchL = [ int( int(pitchFolder) ) for pitchFolder in pitchFolderL ]
  503. mapD = {} # pitch:{ loDb: { hiDb, us_avg, us_cls, us_std, us_min, us_max, db_avg, db_std, cnt }}
  504. # where: cnt=count of valid sample points in this db range
  505. # us_cls=us of closest point to center of db range
  506. dbS = set() # { (loDb,hiDb) } track the set of db ranges
  507. for pitch in pitchL:
  508. print(pitch)
  509. # get the sample measurements for pitch
  510. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, pitch, cfg.analysisArgs['rmsAnalysisArgs'] )
  511. # calc the fit to local straight line curve fit at each point
  512. scoreV = np.abs( rms_analysis.samples_to_linear_residual( usL, dbL) * 100.0 / dbL )
  513. minDurMs = cfg.analysisArgs['resampleMinDurMs']
  514. minDb = cfg.analysisArgs['resampleMinDb']
  515. noiseLimitPct = cfg.analysisArgs['resampleNoiseLimitPct']
  516. # get the set of samples that are not valid (too short, too quiet, too noisy)
  517. skipIdxL, firstAudibleIdx, firstNonSkipIdx = get_dur_skip_indexes( durMsL, dbL, takeIdL, scoreV.tolist(), minDurMs, minDb, noiseLimitPct )
  518. mapD[ pitch ] = {}
  519. # get the count of db ranges
  520. N = int(round((maxMapDb - minMapDb) / incrMapDb)) + 1
  521. # for each db range
  522. for i in range(N):
  523. loDb = minMapDb + (i*incrMapDb)
  524. hiDb = loDb + incrMapDb
  525. dbS.add((loDb,hiDb))
  526. # get the valid (pulse,db) pairs for this range
  527. u_dL = [(us,db) for i,(us,db) in enumerate(zip(usL,dbL)) if i not in skipIdxL and loDb<=db and db<hiDb ]
  528. us_avg = 0
  529. us_cls = 0
  530. us_std = 0
  531. us_min = 0
  532. us_max = 0
  533. db_avg = 0
  534. db_std = 0
  535. if len(u_dL) == 0:
  536. print("No valid samples for pitch:",pitch," db range:",loDb,hiDb)
  537. else:
  538. us0L,db0L = zip(*u_dL)
  539. if len(us0L) == 1:
  540. us_avg = us0L[0]
  541. us_cls = us_avg
  542. us_min = us_avg
  543. us_max = us_avg
  544. db_avg = db0L[0]
  545. elif len(us0L) > 1:
  546. us_avg = np.mean(us0L)
  547. us_cls = us0L[ np.argmin(np.abs(np.array(db0L)-(loDb - (hiDb-loDb)/2.0 ))) ]
  548. us_min = np.min(us0L)
  549. us_max = np.max(us0L)
  550. us_std = np.std(us0L)
  551. db_avg = np.mean(db0L)
  552. db_std = np.std(db0L)
  553. us_avg = int(round(us_avg))
  554. mapD[pitch][loDb] = { 'hiDb':hiDb, 'us_avg':us_avg, 'us_cls':us_cls, 'us_std':us_std,'us_min':us_min,'us_max':us_max, 'db_avg':db_avg, 'db_std':db_std, 'cnt':len(u_dL) }
  555. return mapD, list(dbS)
  556. def plot_us_to_db_map( inDir, cfg, minMapDb=16.0, maxMapDb=26.0, incrMapDb=1.0, pitchL=None, printDir=None ):
  557. fig,ax = plt.subplots()
  558. mapD, dbRefL = estimate_us_to_db_map( inDir, cfg, minMapDb, maxMapDb, incrMapDb, pitchL )
  559. # for each pitch
  560. for pitch, dbD in mapD.items():
  561. u_dL = [ (d['us_avg'],d['us_cls'],d['db_avg'],d['us_std'],d['us_min'],d['us_max'],d['db_std']) for loDb, d in dbD.items() if d['us_avg'] != 0 ]
  562. if u_dL:
  563. # get the us/db lists for this pitch
  564. usL,uscL,dbL,ussL,usnL,usxL,dbsL = zip(*u_dL)
  565. # plot central curve and std dev's
  566. p = ax.plot(usL,dbL, marker='.', label=str(pitch))
  567. ax.plot(uscL,dbL, marker='x', label=str(pitch), color=p[0].get_color(), linestyle='None')
  568. ax.plot(usL,np.array(dbL)+dbsL, color=p[0].get_color(), alpha=0.3)
  569. ax.plot(usL,np.array(dbL)-dbsL, color=p[0].get_color(), alpha=0.3)
  570. # plot us error bars
  571. for db,us,uss,us_min,us_max in zip(dbL,usL,ussL,usnL,usxL):
  572. ax.plot([us_min,us_max],[db,db], color=p[0].get_color(), alpha=0.3 )
  573. ax.plot([us-uss,us+uss],[db,db], color=p[0].get_color(), alpha=0.3, marker='.', linestyle='None' )
  574. plt.legend()
  575. if printDir:
  576. plt.savefig(os.path.join(printDir,"us_db_map.png"),format="png")
  577. plt.show()
  578. def report_take_ids( inDir ):
  579. pitchDirL = os.listdir(inDir)
  580. for pitch in pitchDirL:
  581. pitchDir = os.path.join(inDir,pitch)
  582. takeDirL = os.listdir(pitchDir)
  583. if len(takeDirL) == 0:
  584. print(pitch," directory empty")
  585. else:
  586. fn = os.path.join(pitchDir,'0','seq.json')
  587. if not os.path.isfile(fn):
  588. print("Missing sequence file:",fn)
  589. else:
  590. with open( fn, "rb") as f:
  591. r = json.load(f)
  592. if len(r['eventTimeL']) != 81:
  593. print(pitch," ",len(r['eventTimeL']))
  594. if len(takeDirL) != 3:
  595. print("***",pitch,len(takeDirL))
  596. def filter_us_db( us0L, db0L ):
  597. us1L = [us0L[-1]]
  598. db1L = [db0L[-1]]
  599. dDb = 0
  600. lastIdx = 0
  601. for i,(us,db) in enumerate(zip( us0L[::-1],db0L[::-1])):
  602. db1 = db1L[-1]
  603. if db < db1 and db1-db >= dDb/2:
  604. dDb = db1 - db
  605. us1L.append(us)
  606. db1L.append(db)
  607. lastIdx = i
  608. lastIdx = len(us0L) - lastIdx - 1
  609. usL = [ us0L[lastIdx] ]
  610. dbL = [ db0L[lastIdx] ]
  611. dDb = 0
  612. for us,db in zip(us0L[lastIdx::],db0L[lastIdx::]):
  613. db1 = dbL[-1]
  614. if db > db1:
  615. dDb = db-db1
  616. usL.append(us)
  617. dbL.append(db)
  618. return usL,dbL
  619. def cache_us_db( inDir, cfg, outFn ):
  620. pitchTakeD = get_pitches_and_takes(inDir)
  621. pitch_usDbD = {}
  622. pitchDirL = os.listdir(inDir)
  623. for pitch,takeIdL in pitchTakeD.items():
  624. pitch = int(pitch)
  625. takeId = takeIdL[-1]
  626. print(pitch)
  627. usL, dbL, durMsL, takeIdL, holdDutyPctL = get_merged_pulse_db_measurements( inDir, pitch, cfg.analysisArgs['rmsAnalysisArgs'], takeId )
  628. pitch_usDbD[pitch] = { 'usL':usL, 'dbL':dbL, 'durMsL':durMsL, 'takeIdL':takeIdL, 'holdDutyPctL': holdDutyPctL }
  629. with open(outFn,"w") as f:
  630. json.dump(pitch_usDbD,f)
  631. def gen_vel_map( inDir, cfg, minMaxDbFn, dynLevelN, cacheFn ):
  632. velMapD = {} # { pitch:[ us ] }
  633. pitchDirL = os.listdir(inDir)
  634. # pitchUsDbD = { pitch:
  635. with open(cacheFn,"r") as f:
  636. pitchUsDbD = json.load(f)
  637. # form minMaxDb = { pitch:(minDb,maxDb) }
  638. with open("minInterpDb.json","r") as f:
  639. r = json.load(f)
  640. minMaxDbD = { pitch:(minDb,maxDb) for pitch,minDb,maxDb in zip(r['pitchL'],r['minDbL'],r['maxDbL']) }
  641. pitchL = sorted( [ int(pitch) for pitch in pitchUsDbD.keys()] )
  642. # for each pitch
  643. for pitch in pitchL:
  644. # get the us/db map for this
  645. d = pitchUsDbD[str(pitch)]
  646. usL, dbL = filter_us_db( d['usL'], d['dbL'] )
  647. #usL = d['usL']
  648. #dbL = np.array(d['dbL'])
  649. dbL = np.array(dbL)
  650. velMapD[pitch] = []
  651. maxDb = minMaxDbD[pitch][1]
  652. minDb = minMaxDbD[pitch][0]
  653. dynLevelN = len(cfg.velTableDbL)
  654. # for each dynamic level
  655. for i in range(dynLevelN):
  656. #db = minDb + (i * (maxDb - minDb)/ dynLevelN)
  657. db = cfg.velTableDbL[i]
  658. usIdx = np.argmin( np.abs(dbL - db) )
  659. velMapD[pitch].append( (usL[ usIdx ],db) )
  660. with open("velMapD.json","w") as f:
  661. json.dump(velMapD,f)
  662. mtx = np.zeros((len(velMapD),dynLevelN+1))
  663. print(mtx.shape)
  664. for i,(pitch,usDbL) in enumerate(velMapD.items()):
  665. for j in range(len(usDbL)):
  666. mtx[i,j] = usDbL[j][1]
  667. fig,ax = plt.subplots()
  668. ax.plot(pitchL,mtx)
  669. plt.show()
  670. if __name__ == "__main__":
  671. printDir = os.path.expanduser("~/temp") # os.path.expanduser( "~/src/picadae_ac_3/doc")
  672. cfgFn = sys.argv[1]
  673. inDir = sys.argv[2]
  674. mode = sys.argv[3]
  675. if len(sys.argv) <= 4:
  676. pitchL = None
  677. else:
  678. pitchL = [ int(sys.argv[i]) for i in range(4,len(sys.argv)) ]
  679. cfg = parse_yaml_cfg( cfgFn )
  680. if mode == 'us_db':
  681. plot_us_db_curves_main( inDir, cfg, pitchL, plotTakesFl=True,usMax=None, printDir=printDir )
  682. elif mode == 'us_db_pitch_take':
  683. plot_us_db_takes( inDir, cfg, pitchL, printDir=printDir)
  684. elif mode == 'us_db_pitch_last':
  685. plot_us_db_takes_last( inDir, cfg, pitchL, printDir=printDir)
  686. elif mode == 'noise':
  687. plot_all_noise_curves( inDir, cfg, pitchL )
  688. elif mode == 'min_max':
  689. plot_min_max_db( inDir, cfg, pitchL, printDir=printDir )
  690. elif mode == 'min_max_2':
  691. takeId = pitchL[-1]
  692. del pitchL[-1]
  693. plot_min_max_2_db( inDir, cfg, pitchL, takeId=takeId, printDir=printDir )
  694. elif mode == 'us_db_map':
  695. plot_us_to_db_map( inDir, cfg, pitchL=pitchL, printDir=printDir )
  696. elif mode == 'audacity':
  697. rms_analysis.write_audacity_label_files( inDir, cfg.analysisArgs['rmsAnalysisArgs'] )
  698. elif mode == 'rpt_take_ids':
  699. report_take_ids( inDir )
  700. elif mode == 'manual_db':
  701. plot_min_db_manual( inDir, cfg, printDir=printDir )
  702. elif mode == 'gen_vel_map':
  703. gen_vel_map( inDir, cfg, "minInterpDb.json", 12, "cache_us_db.json" )
  704. elif mode == 'cache_us_db':
  705. cache_us_db( inDir, cfg, "cache_us_db.json")
  706. else:
  707. print("Unknown mode:",mode)