74 lines
2.8 KiB
Python
74 lines
2.8 KiB
Python
#!/usr/bin/env python3
|
|
"""unify.py — global unified-law fit: cut = alpha*ln1p(am*s/res^p / beta)+c
|
|
across multi6 (6 peaks) + tt-series (8 points, both peaks)."""
|
|
import numpy as np, glob, wave
|
|
from scipy.optimize import least_squares
|
|
|
|
def deepest(dirname,b):
|
|
best=None
|
|
for fn in sorted(glob.glob(dirname+'/ph*.npz')):
|
|
d=np.load(fn)
|
|
if '0x540628' not in d.files: continue
|
|
s=d['0x540628'].astype(np.float64)
|
|
if best is None or s[b]<best[0]: best=(s[b],s)
|
|
return best[1] if best else None
|
|
|
|
def loadwav(p):
|
|
w=wave.open(p,'rb'); n=w.getnframes(); ch=w.getnchannels(); sw=w.getsampwidth()
|
|
d=w.readframes(n); w.close()
|
|
if sw==2: return np.frombuffer(d,dtype=np.int16).astype(np.float64).reshape(-1,ch).mean(1)/32768
|
|
raw=np.frombuffer(d,dtype=np.uint8).reshape(-1,ch,3)
|
|
s=raw[:,:,0].astype(np.int64)|(raw[:,:,1].astype(np.int64)<<8)|(raw[:,:,2].astype(np.int64)<<16)
|
|
return np.where(s>=0x800000,s-0x1000000,s).astype(np.float64).reshape(-1,ch).mean(1)/8388608
|
|
|
|
def ta(x,f,sr=44100,L=None):
|
|
x=x[-L:]; t=np.arange(len(x))/sr; c=np.cos(2*np.pi*f*t); sn=np.sin(2*np.pi*f*t)
|
|
return np.hypot(2*(x*c).sum(),2*(x*sn).sum())/len(x)
|
|
|
|
pts=[]
|
|
# multi6 six peaks
|
|
S=deepest('/tmp/opencode/sc_multi6',85)
|
|
cut=np.maximum(-S*8.685889638,0)
|
|
am=np.zeros(len(S)); rs=np.zeros(len(S))
|
|
for ln in open('/tmp/opencode/tract_multi6.txt'):
|
|
if ln.startswith('#'): continue
|
|
p=ln.split(); kk=int(p[0])
|
|
if kk<len(S): am[kk]=float(p[1]); rs[kk]=float(p[2])
|
|
for kk in range(1,len(S)):
|
|
if cut[kk]>1.0 and am[kk]>0.01:
|
|
pts.append(('m6',am[kk],rs[kk],cut[kk]))
|
|
# tt-series 4 drives x 2 peaks (ref-domain cuts via Goertzel)
|
|
inp0=loadwav('/tmp/opencode/tt0_in.wav')
|
|
rc={0:(10.34,11.84),6:(8.12,9.57),12:(6.12,7.45),18:(4.40,5.55)}
|
|
for att,(c1v,c2v) in rc.items():
|
|
tf='/tmp/opencode/tract_tt%d.txt'%att
|
|
d43=d171=None; r43=r171=None
|
|
for ln in open(tf):
|
|
if ln.startswith('#'): continue
|
|
p=ln.split(); kk=int(p[0])
|
|
if kk==43: r43=float(p[2]); d43=am0=None; d43=float(p[1])
|
|
if kk==171: r171=float(p[2]); d171=float(p[1])
|
|
pts.append(('tt%d-pk1'%att,d43,r43,c1v))
|
|
pts.append(('tt%d-pk2'%att,d171,r171,c2v))
|
|
print('points:',len(pts))
|
|
|
|
def resid(p):
|
|
al,be,cc,pp,ss=p
|
|
out=[]
|
|
for nm,a,r,c in pts:
|
|
X=a*ss/max(r**pp,1e-12)
|
|
out.append(al*np.log1p(X/be)+cc-c)
|
|
return np.array(out)
|
|
|
|
p0=[3.28,0.56,0.77,1.25,1.0]
|
|
lb=[0.5,0.01,-5,0.05,0.05]
|
|
ub=[8,5,5,4,50]
|
|
r=least_squares(resid,p0,bounds=(lb,ub),max_nfev=8000)
|
|
mm=-resid(r.x)+np.array([c for _,_,_,c in pts])
|
|
tgt=np.array([c for _,_,_,c in pts])
|
|
rms=np.sqrt(((mm-tgt)**2).mean())
|
|
print('UNIFIED: alpha=%.3f beta=%.4f c=%+.3f p=%.3f s=%.3f ; rms=%.4f max=%.3f' % (
|
|
*r.x,rms,np.abs(mm-tgt).max()))
|
|
for (nm,a,rr,c),pv in zip(pts,mm):
|
|
print(' %-10s am=%.4f res=%.4f cut=%6.2f pred=%6.2f err=%+.2f' % (nm,a,rr,c,pv,c-pv))
|