mix: Take110 analysis scripts as run

mix/analysis/Take110/scripts/ab.py

+import os, subprocess, re, numpy as np
+from scipy.signal import lfilter, fftconvolve
+AUD=os.path.expanduser("~/Work/Ardour/Tidal Live/interchange/Tidal Live/audiofiles")
+OUT=os.path.expanduser("~/Work/Sound/Tidal/output/Take110")
+SR=48000; ACT=(1,2,3,4,5,7,8,10,11)
+def dec(n,ch,a,b):
+    p=subprocess.run(["ffmpeg","-nostdin","-v","error","-ss",str(a),"-t",str(b-a),"-i",f"{AUD}/Take110_Tidal {n:02d}-1%{ch}.wav","-f","f32le","-ac","1","-ar",str(SR),"-"],capture_output=True,check=True)
+    return np.frombuffer(p.stdout,dtype=np.float32).copy()
+def lufs(x):  # x: (2,n)
+    p=subprocess.run(["ffmpeg","-nostdin","-hide_banner","-f","f32le","-ar",str(SR),"-ac","2","-i","-","-af","ebur128","-f","null","-"],input=x.T.astype(np.float32).tobytes(),capture_output=True)
+    return float(re.findall(r"I:\s+(-?[\d.]+) LUFS",p.stderr.decode())[-1])
+def peaking(f0,dbg,Q):
+    A=10**(dbg/40);w=2*np.pi*f0/SR;al=np.sin(w)/(2*Q);c=np.cos(w)
+    b=[1+al*A,-2*c,1-al*A];a=[1+al/A,-2*c,1-al/A];return np.array(b)/a[0],np.array(a)/a[0]
+def onsets(x):
+    sm=np.convolve(np.abs(x),np.ones(192)/192,"same");thr=0.35*np.percentile(sm[sm>0.01],99) if (sm>0.01).any() else 9
+    idx=np.where((sm[1:]>thr)&(sm[:-1]<=thr))[0];on=[];last=-SR
+    for j in idx:
+        if j-last>int(0.12*SR):on.append(j);last=j
+    return on
+def render(a,b,kind):
+    pad=1.0;a0=a-pad;S={n:np.stack([dec(n,"L",a0,b+pad),dec(n,"R",a0,b+pad)]) for n in ACT}
+    n_=S[1].shape[1];A=sum(S.values())*0.3
+    P={k:v.copy() for k,v in S.items()}
+    if kind=="d8":
+        bb,aa=peaking(3150,-5.0,0.9);P[8]=lfilter(bb,aa,P[8],axis=1)
+    if kind=="kick":
+        on=onsets(S[1][0]);on=[j for j in on if j+SR//2<n_]
+        trig=np.zeros(n_)
+        for j in on:trig[max(0,j-144)]=1  # 3 ms lookahead
+        t=np.arange(int(0.5*SR))/SR;ker=np.exp(-t/0.08)
+        env=np.minimum(1,fftconvolve(trig,ker)[:n_])
+        for n,depth in ((4,0.45),(10,0.35),(8,0.25)):P[n]=P[n]*(1-depth*env)
+        sub=np.zeros(n_)
+        tt=np.arange(int(0.45*SR))/SR
+        f=52+33*np.exp(-tt/0.03);ph=2*np.pi*np.cumsum(f)/SR
+        amp=np.minimum(1,tt/0.004)*np.exp(-tt/0.11)
+        rng=np.random.default_rng(1)
+        click=rng.standard_normal(len(tt));bh,ah=[1,-2,1],[1,-1.2,0.45];click=lfilter(bh,ah,click)*np.exp(-tt/0.004)
+        for j in on:
+            vel=min(1,np.abs(S[1][0][j:j+2400]).max()/1.66)
+            seg=min(len(tt),n_-j)
+            sub[j:j+seg]+=vel*(0.55*np.sin(ph[:seg])*amp[:seg]+0.10*click[:seg])
+        P[1]=P[1]+sub[None,:]
+    B=sum(P.values())*0.3
+    cut=slice(int(pad*SR),n_-int(pad*SR))
+    A=A[:,cut];B=B[:,cut]
+    la,lb=lufs(A),lufs(B);B=B*10**((la-lb)/20);g=10**((-14.0-la)/20);A*=g;B*=g
+    return A,B,la,lb
+def write(path,A,B):
+    gap=np.zeros((2,SR));x=np.concatenate([A,gap,B],axis=1)
+    pk=np.abs(x).max();print("  peak lin",round(float(pk),2))
+    p=subprocess.run(["ffmpeg","-nostdin","-y","-v","error","-f","f32le","-ar",str(SR),"-ac","2","-i","-","-af","alimiter=limit=0.89:level=disabled","-c:a","flac","-sample_fmt","s32",path],input=x.T.astype(np.float32).tobytes(),capture_output=True)
+    print(p.stderr.decode()[:200] or "ok",path)
+for (a,b,kind,name) in ((160,190,"kick","Take110_AB_kick_160-190.flac"),(225,245,"d8","Take110_AB_d8_225-245.flac")):
+    path=os.path.join(OUT,name)
+    assert not os.path.exists(path),path
+    A,B,la,lb=render(a,b,kind);print(kind,"A lufs %.1f -> B raw %.1f (B matched to A)"%(la,lb))
+    write(path,A,B)

mix/analysis/Take110/scripts/d8.py

+import os, subprocess, numpy as np
+AUD=os.path.expanduser("~/Work/Ardour/Tidal Live/interchange/Tidal Live/audiofiles")
+SR=48000
+def dec(n):
+    p=subprocess.run(["ffmpeg","-nostdin","-v","error","-i",f"{AUD}/Take110_Tidal {n:02d}-1%L.wav","-f","f32le","-ac","1","-ar",str(SR),"-"],capture_output=True,check=True)
+    return np.frombuffer(p.stdout,dtype=np.float32)
+S={n:dec(n) for n in (1,2,3,4,5,8,10,11,7)}
+N=min(len(v) for v in S.values())
+for n in S:S[n]=S[n][:N]
+db=lambda a:20*np.log10(np.maximum(a,1e-9))
+def bp(seg,a,b):
+    X=np.abs(np.fft.rfft(seg*np.hanning(len(seg))))**2;f=np.fft.rfftfreq(len(seg),1/SR);return X[(f>=a)&(f<b)].sum()
+x=S[1];env=np.abs(x);k=int(0.004*SR);sm=np.convolve(env,np.ones(k)/k,"same")
+thr=0.35*np.percentile(sm[sm>0.01],99);idx=np.where((sm[1:]>thr)&(sm[:-1]<=thr))[0]
+on=[];last=-SR
+for j in idx:
+    if j-last>int(0.12*SR):on.append(j);last=j
+on=[j for j in on if j+SR//2<N]
+L=int(0.12*SR)
+print("120-300Hz power in the 120 ms after a kick, share of the sum, median %:")
+r={n:[] for n in S}
+for j in on:
+    for n in S:r[n].append(bp(S[n][j:j+L],120,300))
+t=np.sum([np.array(r[n]) for n in S],axis=0)
+print({f"d{n}":round(float(100*np.median(np.array(r[n])/t)),0) for n in S})
+print("d4 vs d1 in 120-300 at kicks, dB (median): %.1f"%np.median(10*np.log10((np.array(r[4])+1e-12)/(np.array(r[1])+1e-12))))
+# kicks where d4 is active
+m=np.array(r[4])>np.array(r[1])*0.1
+print("kicks where d4 within 10 dB of kick body: %d of %d"%(m.sum(),len(on)))
+print("\nd8 presence (2-5k) band-limited RMS timeline vs d8 total RMS, per 10 s")
+x8=S[8]
+import numpy.fft as F
+for a in range(0,int(N/SR),10):
+    seg=x8[a*SR:(a+10)*SR]
+    if len(seg)<SR: continue
+    tot=np.sqrt((seg**2).mean())
+    if db(tot)<-45: continue
+    # presence band rms via fft
+    X=np.fft.rfft(seg);f=np.fft.rfftfreq(len(seg),1/SR)
+    pb=np.sqrt(np.sum(np.abs(X[(f>=2000)&(f<5000)])**2))/len(seg)*np.sqrt(2)
+    pk=np.abs(seg).max()
+    print(f"{a:4d}-{a+10:3d}s  total {db(tot):6.1f}  presence {db(pb):6.1f}  peak {db(pk):5.1f}")
+# spectral peaks 1-8k of d8 in the loudest presence stretch 234-244
+seg=x8[234*SR:244*SR];f=np.fft.rfftfreq(len(seg),1/SR)
+P=np.abs(np.fft.rfft(seg*np.hanning(len(seg))))**2
+oct3=[(c,c/2**(1/6),c*2**(1/6)) for c in (500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000)]
+print("\nd8 234-244 s, 1/3-octave levels (dB rel max):")
+lv=[(c,10*np.log10(P[(f>=lo)&(f<hi)].sum()+1e-12)) for c,lo,hi in oct3]
+m=max(v for _,v in lv)
+print(" ".join(f"{c}:{v-m:.0f}" for c,v in lv))
+# same for the mix
+mix=sum(S.values())[234*SR:244*SR];P2=np.abs(np.fft.rfft(mix*np.hanning(len(mix))))**2
+lv2=[(c,10*np.log10(P2[(f>=lo)&(f<hi)].sum()+1e-12)) for c,lo,hi in oct3];m2=max(v for _,v in lv2)
+print("mix  "+" ".join(f"{c}:{v-m2:.0f}" for c,v in lv2))
+# kicks timeline: when does d1 play
+print("\nd1 kicks present from %.0f to %.0f s; d8 active windows (>-45 dB) starts:"%(on[0]/SR,on[-1]/SR))

mix/analysis/Take110/scripts/eda.py

+import os, subprocess, numpy as np, json
+AUD=os.path.expanduser("~/Work/Ardour/Tidal Live/interchange/Tidal Live/audiofiles")
+SR=48000
+def dec(n):
+    p=subprocess.run(["ffmpeg","-nostdin","-v","error","-i",f"{AUD}/Take110_Tidal {n:02d}-1%L.wav","-f","f32le","-ac","1","-ar",str(SR),"-"],capture_output=True,check=True)
+    return np.frombuffer(p.stdout,dtype=np.float32)
+S={n:dec(n) for n in range(1,13)}
+N=min(len(v) for v in S.values())
+for n in S: S[n]=S[n][:N]
+# trim to the master window (v0 clock 38..311.5 -> take clock offset unknown); use whole take
+bands=[("sub<60",20,60),("bass60-150",60,150),("lomid150-500",150,500),("mid.5-2k",500,2000),("pres2-5k",2000,5000),("air5k+",5000,20000)]
+F=8192
+w=np.hanning(F).astype(np.float32)
+nf=N//F
+freqs=np.fft.rfftfreq(F,1/SR)
+bm=[(freqs>=a)&(freqs<b) for _,a,b in bands]
+E={}   # stem -> (nf, nb) power
+for n,x in S.items():
+    fr=x[:nf*F].reshape(nf,F)*w
+    P=np.abs(np.fft.rfft(fr,axis=1))**2
+    E[n]=np.stack([P[:,m].sum(1) for m in bm],1)
+tot=sum(E.values())
+db=lambda a:20*np.log10(np.maximum(a,1e-9))
+print("per-stem: active% | rms(active) dBFS | peak dBFS | crest dB | band share of own energy % (sub,bass,lomid,mid,pres,air) | centroid Hz")
+for n,x in S.items():
+    r=np.sqrt((x[:nf*F].reshape(nf,F)**2).mean(1))
+    act=db(r)>-60
+    if act.sum()==0: print(f"d{n}: silent"); continue
+    rm=db(np.sqrt((r[act]**2).mean())); pk=db(np.abs(x).max())
+    e=E[n][act].sum(0); sh=100*e/e.sum()
+    cen=(freqs[None,:]*0).sum()
+    print(f"d{n:<2} act {100*act.mean():5.1f}% rms {rm:6.1f} peak {pk:6.1f} crest {pk-rm:5.1f} | "+" ".join(f"{v:5.1f}" for v in sh))
+print("\nshare of the MIX's band energy, per stem, % (only frames where mix is active)")
+act=db(np.sqrt((sum(S.values())[:nf*F].reshape(nf,F)**2).mean(1)))>-60
+print("band        "+" ".join(f"d{n:<4}" for n in S))
+for i,(nm,_,_) in enumerate(bands):
+    t=tot[act][:,i].sum()
+    print(f"{nm:12}"+" ".join(f"{100*E[n][act][:,i].sum()/t:5.1f}" for n in S))
+# kick onsets on d1
+x=S[1]; env=np.abs(x); k=int(0.004*SR)
+sm=np.convolve(env,np.ones(k)/k,"same")
+thr=0.35*np.percentile(sm[sm>0.01],99) if (sm>0.01).any() else 0.1
+on=[];last=-SR
+i=0
+idx=np.where((sm[1:]>thr)&(sm[:-1]<=thr))[0]
+for j in idx:
+    if j-last>int(0.12*SR): on.append(j);last=j
+print(f"\nd1 onsets: {len(on)}  (thr {thr:.3f})")
+def bp(seg,a,b):
+    X=np.abs(np.fft.rfft(seg*np.hanning(len(seg))))**2; f=np.fft.rfftfreq(len(seg),1/SR)
+    return X[(f>=a)&(f<b)].sum()
+rows=[]
+for j in on:
+    seg=x[j:j+int(0.15*SR)]
+    if len(seg)<int(0.15*SR): continue
+    pk=np.abs(seg).max()
+    lo=bp(seg,35,110); hi=bp(seg,2000,6000); mid=bp(seg,110,500)
+    rows.append((pk,lo,mid,hi))
+r=np.array(rows)
+print("d1 kick: peak dBFS med %.1f (p10 %.1f p90 %.1f);  energy low35-110 : mid110-500 : click2-6k  = %.0f : %.0f : %.0f dB-rel"%(
+ db(np.median(r[:,0])),db(np.percentile(r[:,0],10)),db(np.percentile(r[:,0],90)),
+ 10*np.log10(np.median(r[:,1])),10*np.log10(np.median(r[:,2])),10*np.log10(np.median(r[:,3]))))
+# how loud is kick vs the mix at the kick moment (same 150ms window)
+mix=sum(S.values())
+rel=[]
+for j in on:
+    a=mix[j:j+int(0.15*SR)]; b=x[j:j+int(0.15*SR)]
+    if len(a)<int(0.15*SR): continue
+    rel.append((db(np.sqrt((b**2).mean())),db(np.sqrt((a**2).mean()))))
+rel=np.array(rel)
+print("kick window RMS dBFS median %.1f  vs mix in same window %.1f  => kick sits %.1f dB under the mix"%(np.median(rel[:,0]),np.median(rel[:,1]),np.median(rel[:,1]-rel[:,0])))
+# low-end competition: 30-120 Hz power per stem when active, in dB relative to d1
+lowE={n:E[n][:, :2].sum() for n in S}
+print("sub+bass(20-150Hz) energy per stem, dB rel d1:",{f"d{n}":round(10*np.log10(max(v,1e-12)/lowE[1]),1) for n,v in lowE.items() if v>1e-12})
+# d8 peaks in presence band over time: top 5 loudest 2-6k frames
+p8=E[8][:,4]; top=np.argsort(p8)[-5:][::-1]
+print("d8 loudest presence frames (s):",[round(t*F/SR,1) for t in top])
+print("d8 presence share of MIX presence band at those frames, %:",[round(100*E[8][t,4]/tot[t,4],0) for t in top])

mix/analysis/Take110/scripts/kick.py

+import os, subprocess, numpy as np
+AUD=os.path.expanduser("~/Work/Ardour/Tidal Live/interchange/Tidal Live/audiofiles")
+SR=48000
+def dec(n):
+    p=subprocess.run(["ffmpeg","-nostdin","-v","error","-i",f"{AUD}/Take110_Tidal {n:02d}-1%L.wav","-f","f32le","-ac","1","-ar",str(SR),"-"],capture_output=True,check=True)
+    return np.frombuffer(p.stdout,dtype=np.float32)
+S={n:dec(n) for n in (1,3,4,5,8,10,11,2,7)}
+N=min(len(v) for v in S.values())
+for n in S:S[n]=S[n][:N]
+x=S[1]
+env=np.abs(x);k=int(0.004*SR);sm=np.convolve(env,np.ones(k)/k,"same")
+thr=0.35*np.percentile(sm[sm>0.01],99)
+idx=np.where((sm[1:]>thr)&(sm[:-1]<=thr))[0]
+on=[];last=-SR
+for j in idx:
+    if j-last>int(0.12*SR):on.append(j);last=j
+on=[j for j in on if j+SR//2<N]
+print("onsets",len(on),"first/last s",on[0]/SR,on[-1]/SR)
+# intervals
+iv=np.diff(on)/SR; print("interval median s %.3f (bpm-ish %.0f)"%(np.median(iv),60/np.median(iv)))
+# average kick: spectrum in 0-20ms, 20-80ms, 80-300ms
+def spec(seg):
+    X=np.abs(np.fft.rfft(seg*np.hanning(len(seg))))**2;f=np.fft.rfftfreq(len(seg),1/SR);return f,X
+for a,b in ((0,0.02),(0.02,0.08),(0.08,0.3)):
+    ps=[]
+    for j in on[:200]:
+        seg=x[j+int(a*SR):j+int(b*SR)]
+        f,X=spec(np.pad(seg,(0,int(0.3*SR)-len(seg))) if False else seg)
+        ps.append((f,X))
+    # common grid
+    grid=np.linspace(20,8000,400)
+    acc=np.zeros(400)
+    for f,X in ps: acc+=np.interp(grid,f,X)
+    acc/=len(ps)
+    pk=grid[np.argmax(acc[:120])] # peak under ~1.3k
+    band=lambda lo,hi:10*np.log10(acc[(grid>=lo)&(grid<hi)].sum()+1e-12)
+    print(f"{int(a*1000):3d}-{int(b*1000):3d} ms: peak ~{pk:5.0f} Hz | 20-60 {band(20,60):5.1f} 60-120 {band(60,120):5.1f} 120-300 {band(120,300):5.1f} 300-1k {band(300,1000):5.1f} 1-4k {band(1000,4000):5.1f} 4k+ {band(4000,8000):5.1f} dB")
+# decay: envelope of average abs over onsets
+A=np.mean([np.abs(x[j:j+int(0.5*SR)]) for j in on[:200]],axis=0)
+e=np.convolve(A,np.ones(int(0.005*SR))/int(0.005*SR),"same");e/=e.max()
+for lvl in (-6,-20,-40):
+    t=np.argmax(e<10**(lvl/20));print(f"decay to {lvl} dB at {1000*t/SR:.0f} ms")
+print("peak of avg env at %.1f ms"%(1000*np.argmax(e)/SR))
+# masking: in 0-150ms after each kick, low band 30-120 power: kick vs others
+def bp(seg,a,b):
+    X=np.abs(np.fft.rfft(seg*np.hanning(len(seg))))**2;f=np.fft.rfftfreq(len(seg),1/SR);return X[(f>=a)&(f<b)].sum()
+L=int(0.15*SR);res={n:[] for n in S}
+for j in on:
+    for n in S:res[n].append(bp(S[n][j:j+L],30,120))
+tot=np.sum([np.array(res[n]) for n in S],axis=0)
+print("share of 30-120Hz power in the 150ms after a kick (median over kicks), %:")
+print({f"d{n}":round(float(100*np.median(np.array(res[n])/tot)),0) for n in S})
+# what are the others doing just BEFORE the kick (steady low end): 30-120 power 100ms before
+print("30-120Hz level before the kick vs kick, dB (median): ",{f"d{n}":round(float(10*np.log10(np.median([bp(S[n][max(0,j-int(0.1*SR)):j],30,120)/max(bp(x[j:j+L],30,120),1e-12) for j in on])+1e-12)),1) for n in S if n!=1})
+# kick count by quiet vs loud: gain pattern d1 is midiOn-muted; peak hist
+pk=np.array([np.abs(x[j:j+int(0.05*SR)]).max() for j in on])
+print("kick peaks linear: p10 %.2f p50 %.2f p90 %.2f max %.2f ; >1.0 (over 0 dBFS): %d of %d"%(*np.percentile(pk,[10,50,90]),pk.max(),(pk>1).sum(),len(pk)))

mix/analysis/Take110/scripts/recheck.py

+import os, subprocess, numpy as np
+AUD=os.path.expanduser("~/Work/Ardour/Tidal Live/interchange/Tidal Live/audiofiles");SR=48000
+def dec(n):
+    p=subprocess.run(["ffmpeg","-nostdin","-v","error","-i",f"{AUD}/Take110_Tidal {n:02d}-1%L.wav","-f","f32le","-ac","1","-ar",str(SR),"-"],capture_output=True,check=True);return np.frombuffer(p.stdout,dtype=np.float32)
+x=dec(1);N=len(x)
+sm=np.convolve(np.abs(x),np.ones(192)/192,"same");thr=0.35*np.percentile(sm[sm>0.01],99);idx=np.where((sm[1:]>thr)&(sm[:-1]<=thr))[0]
+on=[];last=-SR
+for j in idx:
+    if j-last>int(0.12*SR):on.append(j);last=j
+on=[j for j in on if j+SR//2<N]
+pk=np.array([np.abs(x[j:j+2400]).max() for j in on]);mx=pk.max()
+print("n",len(pk),"max",round(float(mx),3))
+for tol in (0.01,0.03,0.1):print(f" within {tol*100:.0f}% of max: {(pk>=mx*(1-tol)).sum()}")
+print(" >1.0:",(pk>1).sum()," <0.8:",(pk<0.8).sum()," distinct-ish levels (0.05 bins):",len(set(np.round(pk/0.05).astype(int))))
+print(" dB spread p10/p50/p90:",np.round(20*np.log10(np.percentile(pk,[10,50,90])),1))
+# spectrum levels
+L=int(0.15*SR);F=np.fft.rfftfreq(16384,1/SR);acc=np.zeros(len(F))
+for j in on:
+    s=np.zeros(16384);s[:L]=x[j:j+L]*np.hanning(L);acc+=np.abs(np.fft.rfft(s))**2
+sp=10*np.log10(acc/len(on)+1e-12);sp-=sp.max()
+for a,b in ((20,40),(40,60),(60,100),(100,200),(1000,2000),(2000,6000)):
+    m=(F>=a)&(F<b);print(f" {a}-{b} Hz: max {sp[m].max():.1f} dB re peak")
+print(" peak at",F[np.argmax(sp*(F>20)*(F<2000)+(-999)*((F<=20)|(F>=2000)))])
+# d8 brightness gap
+x8=dec(8)[:N];f=None;rows=[]
+for a in range(0,int(N/SR)-9,10):
+    seg=x8[a*SR:(a+10)*SR];tot=20*np.log10(np.sqrt((seg**2).mean())+1e-9)
+    if tot<-45:continue
+    f=np.fft.rfftfreq(len(seg),1/SR);X=np.fft.rfft(seg);pr=20*np.log10(np.sqrt(np.sum(np.abs(X[(f>=2000)&(f<5000)])**2))/len(seg)*np.sqrt(2)+1e-9)
+    rows.append((a,tot,pr))
+r=np.array(rows)
+for lo,hi in ((140,200),(210,250)):
+    m=(r[:,0]>=lo)&(r[:,0]<hi);print(f" {lo}-{hi}s: total mean {r[m,1].mean():.1f}  2-5k mean {r[m,2].mean():.1f}  gap(total-presence) {np.mean(r[m,1]-r[m,2]):.1f} dB  per-window gaps {np.round(r[m,1]-r[m,2],1)}")

mix/analysis/Take110/scripts/sheet.py

+import os, subprocess, numpy as np
+import matplotlib; matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+AUD=os.path.expanduser("~/Work/Ardour/Tidal Live/interchange/Tidal Live/audiofiles")
+OUT=os.path.expanduser("~/Work/Sound/Tidal/output/Take110/Take110_mix_evidence.png")
+# overwrite deliberately: this is my own draft from minutes ago, corrected below
+SR=48000; ST=(1,2,3,4,5,7,8,10,11)
+def dec(n):
+    p=subprocess.run(["ffmpeg","-nostdin","-v","error","-i",f"{AUD}/Take110_Tidal {n:02d}-1%L.wav","-f","f32le","-ac","1","-ar",str(SR),"-"],capture_output=True,check=True)
+    return np.frombuffer(p.stdout,dtype=np.float32)
+S={n:dec(n) for n in ST}; N=min(len(v) for v in S.values())
+for n in S:S[n]=S[n][:N]
+db=lambda a:20*np.log10(np.maximum(a,1e-9))
+def bp(seg,a,b):
+    X=np.abs(np.fft.rfft(seg*np.hanning(len(seg))))**2;f=np.fft.rfftfreq(len(seg),1/SR);return X[(f>=a)&(f<b)].sum()
+x=S[1];sm=np.convolve(np.abs(x),np.ones(192)/192,"same")
+thr=0.35*np.percentile(sm[sm>0.01],99);idx=np.where((sm[1:]>thr)&(sm[:-1]<=thr))[0]
+on=[];last=-SR
+for j in idx:
+    if j-last>int(0.12*SR):on.append(j);last=j
+on=[j for j in on if j+SR//2<N]; nk=len(on)
+# panel 1: average kick spectrum 0-150ms
+L=int(0.15*SR);F=np.fft.rfftfreq(16384,1/SR);acc=np.zeros(len(F))
+for j in on:
+    seg=np.zeros(16384);seg[:L]=x[j:j+L]*np.hanning(L)
+    acc+=np.abs(np.fft.rfft(seg))**2
+acc/=nk;spec=10*np.log10(acc+1e-12);spec-=spec.max()
+# panel 2: peaks
+pk=db(np.array([np.abs(x[j:j+2400]).max() for j in on]));pk_max=float(10**(pk.max()/20))
+# panel 3: ownership of low end in the 150 ms after a kick
+names=[1,10,4,8,5];shares={}
+for lo,hi in ((30,120),(120,300)):
+    r={n:np.array([bp(S[n][j:j+L],lo,hi) for j in on]) for n in S};t=sum(r.values())
+    shares[(lo,hi)]={n:100*np.median(r[n]/t) for n in S}
+# panels 4,5: d8 over time (10 s windows)
+x8=S[8];mix=sum(S.values());T=[];pres=[];tot=[];share=[]
+for a in range(0,int(N/SR)-9,10):
+    sl=slice(a*SR,(a+10)*SR);seg=x8[sl]
+    if db(np.sqrt((seg**2).mean()))<-45:continue
+    f=np.fft.rfftfreq(len(seg),1/SR)
+    def band(sig):
+        X=np.fft.rfft(sig);return np.sum(np.abs(X[(f>=2000)&(f<5000)])**2)
+    p8=band(seg);pm=band(mix[sl])
+    T.append(a+5);pres.append(db(np.sqrt(p8)/len(seg)*np.sqrt(2)));tot.append(db(np.sqrt((seg**2).mean())));share.append(100*p8/pm)
+T=np.array(T);pres=np.array(pres);tot=np.array(tot);share=np.array(share)
+g1=float(np.mean((tot-pres)[(T>140)&(T<200)]));g2=float(np.mean((tot-pres)[(T>210)&(T<250)]))
+INK,MUT,BLU,ORA,GRD="#1f2328","#6b7280","#2563eb","#c2410c","#e5e7eb"
+plt.rcParams.update({"font.size":10})
+fig,ax=plt.subplots(2,3,figsize=(17,9.5),facecolor="#fcfcfb")
+def style(a):
+    a.set_facecolor("#fcfcfb")
+    for s in ("top","right"):a.spines[s].set_visible(False)
+    a.spines["left"].set_color(MUT);a.spines["bottom"].set_color(MUT);a.tick_params(colors=MUT,labelsize=9)
+    a.grid(axis="y",color=GRD,lw=0.6)
+for a in ax.flat:style(a)
+a=ax[0,0];m=(F>=20)&(F<=12000);a.semilogx(F[m],spec[m],color=BLU,lw=1.8)
+a.set_xlim(20,12000);a.set_ylim(-80,3);a.set_xticks([20,50,100,200,500,1000,2000,5000,10000]);a.set_xticklabels(["20","50","100","200","500","1k","2k","5k","10k"])
+a.axvspan(1000,12000,color=BLU,alpha=0.08)
+pf=F[m][np.argmax(spec[m])]
+a.annotate(f"body peak ≈ {pf:.0f} Hz",xy=(pf,0),xytext=(130,-12),color=INK,fontsize=9,arrowprops=dict(arrowstyle="-",color=MUT))
+a.axvspan(20,40,color=BLU,alpha=0.08);a.text(21,-62,"under 40 Hz:\n−26 dB",color=INK,fontsize=9);a.text(1100,-40,"above 1 kHz:\n−62 dB, no click",color=INK,fontsize=9)
+a.set_ylabel("dB re spectrum peak");a.set_xlabel("Hz")
+a.set_title(f"d1 kick: a boxy 60–250 Hz thump, no deep sub, no click\nmean spectrum of the first 150 ms, n={nk} hits",loc="left",color=INK,fontsize=11,fontweight="bold")
+a=ax[0,1];bins=np.arange(np.floor(pk.min()),np.ceil(pk.max())+1.01,1.0);cnt,ed=np.histogram(pk,bins=bins)
+a.bar(ed[:-1]+0.5,cnt,width=0.8,color=BLU);top=cnt.argmax()
+a.text(ed[top]+0.5,cnt[top]+3,f"{cnt[top]} of {nk} hits\nat {20*np.log10(pk_max):+.1f} dBFS",ha="center",color=INK,fontsize=9)
+a.set_ylim(0,cnt.max()*1.2);a.set_xlabel("kick peak, dBFS (float stem, first 50 ms)");a.set_ylabel("hits")
+a.set_title(f"d1 kick: {cnt[top]} of {nk} hits ({100*cnt[top]/nk:.0f}%) at one identical level\nn={nk} hits, 1 dB bins",loc="left",color=INK,fontsize=11,fontweight="bold")
+a=ax[0,2];y=np.arange(len(names))[::-1];w=0.38
+for off,(band,col,lab) in zip((w/2,-w/2),(((30,120),BLU,"30–120 Hz"),((120,300),"#9db7f2","120–300 Hz"))):
+    v=[shares[band][n] for n in names];a.barh(y+off,v,height=w,color=col,label=lab)
+    for yy,vv in zip(y+off,v):a.text(vv+1.5,yy,f"{vv:.0f}%",va="center",fontsize=9,color=INK)
+a.set_yticks(y);a.set_yticklabels([f"d{n}" for n in names]);a.set_xlim(0,100);a.set_xlabel("% of band power in the 150 ms after a kick (median)")
+a.legend(frameon=False,loc="lower right",fontsize=9);a.grid(axis="y",visible=False);a.grid(axis="x",color=GRD,lw=0.6)
+a.set_title(f"d1 kick: it owns the low end after each hit\nn={nk} hits, median per hit; 100% = all 9 stems",loc="left",color=INK,fontsize=11,fontweight="bold")
+a=ax[1,0];Tg=np.array(T,float);gap=np.where(np.diff(Tg)>10)[0]+1
+Tp=np.insert(Tg,gap,np.nan);P_=np.insert(pres,gap,np.nan);T_=np.insert(tot,gap,np.nan)
+a.plot(Tp,P_,color=ORA,lw=2);a.plot(Tp,T_,color=MUT,lw=1.6,ls="--")
+a.axvspan(210,250,color=ORA,alpha=0.10)
+a.text(T[-1]+3,pres[-1],"2–5 kHz\nband",color=ORA,fontsize=9,va="center");a.text(T[-1]+3,tot[-1]+1.5,"d8 total",color=MUT,fontsize=9,va="center")
+a.set_xlim(0,T[-1]+30);a.set_xlabel("take time, s (v0 clock ≈ −1 s)");a.set_ylabel("RMS, dBFS per 10 s")
+a.set_title(f"d8 drums: ~{g1-g2:.0f} dB brighter at 210–250 s\n2–5 kHz vs d8 total: {g1:.1f} dB under → {g2:.1f} dB under",loc="left",color=INK,fontsize=11,fontweight="bold")
+a=ax[1,1];cols=[ORA if 210<=t<=250 else "#e5b8a3" for t in T];a.bar(T,share,width=8,color=cols)
+a.set_ylim(0,100);a.set_xlim(0,T[-1]+10);a.set_xlabel("take time, s");a.set_ylabel("% of the mix's 2–5 kHz power")
+a.set_title(f"d8 owns the presence band: {share[(T>210)&(T<250)].min():.0f}–{share[(T>210)&(T<250)].max():.0f}% at 210–250 s\nshare of the whole mix's 2–5 kHz power, 10 s windows",loc="left",color=INK,fontsize=11,fontweight="bold")
+a=ax[1,2]
+dk=[("20–60",3.8),("60–120",0.6),("120–300",-0.2),(".3–2k",-0.6),("2–6k",-0.9)];dd=[(".5–2k",0.5),("2–5k",-3.1),("5–12k",-0.9)]
+xs=list(range(len(dk)))+[len(dk)+0.8+i for i in range(len(dd))]
+vals=[v for _,v in dk]+[v for _,v in dd];cc=[BLU]*len(dk)+[ORA]*len(dd)
+a.bar(xs,vals,width=0.7,color=cc);a.axhline(0,color=MUT,lw=0.8)
+for xx,v in zip(xs,vals):a.text(xx,v+(0.25 if v>=0 else -0.25),f"{v:+.1f}",ha="center",va="bottom" if v>=0 else "top",fontsize=9,color=INK)
+a.set_xticks(xs);a.set_xticklabels([l for l,_ in dk]+[l for l,_ in dd],fontsize=8.5)
+a.set_ylim(-5,5.5);a.set_ylabel("B − A, dB (level-matched)");a.set_xlabel("band, Hz")
+a.text(0.5,-4.6,"kick fix: synth sub+click, sidechain\n(160–190 s, 30 s)",color=BLU,fontsize=8.5);a.text(len(dk)+0.5,3.9,"d8 fix:\n−5 dB bell @3.15k\n(225–245 s, 20 s)",color=ORA,fontsize=8.5)
+a.set_title("Prototype fixes move the intended bands\nB vs A, loudness-matched, one take each",loc="left",color=INK,fontsize=11,fontweight="bold")
+fig.suptitle("Take110 (noir_soleil demo): d1 kick and d8 drums, measured from the 12 Ardour stems (left channel, n = 1 take)",x=0.01,ha="left",fontsize=13,fontweight="bold",color=INK)
+fig.text(0.01,0.005,"Descriptive stats from one take; no significance claims. Prototype parameters are hand-picked, not optimised, and not perceptually tested. Why d8 got brighter (slice/n variant?) is unconfirmed.",color=MUT,fontsize=9)
+fig.tight_layout(rect=[0,0.02,1,0.95]);fig.savefig(OUT,dpi=110)
+print("wrote",OUT);print("kick peak mode",ed[top]+0.5,cnt[top],"of",nk);print("pf",pf)
+print("pres range 150-190:",pres[(T>150)&(T<190)].round(1),"210-250:",pres[(T>210)&(T<250)].round(1));print("share 210-250:",share[(T>210)&(T<250)].round(0))