5200b59Write raw float pressure from the render command; add post-processing scripts for matching a recordingJeremy Magland 1import sys, wave, numpy as np
2from scipy.signal import stft
3def load(path):
4 w=wave.open(path); fs=w.getframerate(); n=w.getnframes()
5 x=np.frombuffer(w.readframes(n),dtype='<i2').astype(float)/32768
6 if w.getnchannels()>1: x=x.reshape(-1,w.getnchannels()).mean(1)
7 return fs,x
8def onset(x,fs):
9 e=np.abs(x); i=np.argmax(e>0.3*e.max()); return max(0,i-int(0.005*fs))
10bands=[(30,60),(60,120),(120,250),(250,500),(500,1000),(1000,2000),(2000,4000),(4000,8000),(8000,16000)]
11def analyze(path):
12 fs,x=load(path); o=onset(x,fs); x=x[o:]
13 print(f'== {path}: fs {fs}, {len(x)/fs:.2f} s after onset, peak {np.abs(x).max():.3f} FS')
14 f,t,Z=stft(x,fs,nperseg=4096,noverlap=4096-512)
15 P=np.abs(Z)**2
16 print(' band lvl@0.1s T60 fit(0.2-2s) T60 fit(0.2-end) level@2s')
17 for lo,hi in bands:
18 m=(f>=lo)&(f<hi); e=10*np.log10(P[m].sum(0)+1e-20)
19 i0=np.searchsorted(t,0.1); e0=e[i0]
20 def fit(a,b):
21 ia,ib=np.searchsorted(t,a),np.searchsorted(t,b); ib=min(ib,len(t)-1)
22 if ib-ia<4: return float('nan')
23 k=np.polyfit(t[ia:ib],e[ia:ib],1)[0]; return -60/k if k<0 else float('inf')
24 i2=min(np.searchsorted(t,2.0),len(t)-1)
25 print(f' {lo:5d}-{hi:5d} {e0:7.1f} {fit(0.2,2.0):8.2f} {fit(0.2,t[-1]):8.2f} {e[i2]-e0:7.1f}')
26 # strongest partials in 0.3-1.5 s
27 s=x[int(0.3*fs):int(1.5*fs)]; S=np.abs(np.fft.rfft(s*np.hanning(len(s))))**2; fr=np.fft.rfftfreq(len(s),1/fs)
28 idx=[i for i in range(2,len(S)-2) if S[i]>S[i-1] and S[i]>S[i+1] and fr[i]>25]
29 idx=sorted(idx,key=lambda i:-S[i])[:14]
30 print(' strongest partials (Hz, dB rel max):', ', '.join(f'{fr[i]:.0f}({10*np.log10(S[i]/S[idx[0]]):.0f})' for i in sorted(idx)))
31 return fs,x
32for p in sys.argv[1:]: analyze(p)