From f344abd67cbee006b11813a0ac7f2fca7dd80101 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 00:49:35 +0100 Subject: [PATCH] [FEATURE] add envelope, start_ir, time_shift and other editing utilities Adds to pytta.functions nine editing utilities that the ITA-Toolbox has: envelope (ita_envelope, the magnitude of the analytic signal), start_ir (ita_start_IR, the start of an impulse response by ISO 3382-1, A.3), time_shift (ita_time_shift, cyclic, fractional shifts as a linear phase, or the start of the response to the first sample), append (ita_append), normalize (ita_normalize_dat and ita_normalize_spk), amplify_to (ita_amplify_to, to a mean RMS level), zero_phase, randomize_phase and quantize (ita_zerophase, ita_randomize_phase, ita_quantize). Tests check start_ir and quantize against the ITA-Toolbox on golden data from tests/golden/make_dsp_edit_ita.m (identical), the envelope of an AM tone against its formula and against ITA (1.2e-7), integer and fractional shifts, the peak and level rules, and that the phase edits keep the magnitude spectrum. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/functions.rst | 10 +- pytta/__init__.py | 13 ++- pytta/functions.py | 193 +++++++++++++++++++++++++++++++ tests/golden/dsp_edit_ita.mat | Bin 0 -> 32755 bytes tests/golden/make_dsp_edit_ita.m | 22 ++++ tests/test_dsp_edit.py | 135 +++++++++++++++++++++ 6 files changed, 371 insertions(+), 2 deletions(-) create mode 100644 tests/golden/dsp_edit_ita.mat create mode 100644 tests/golden/make_dsp_edit_ita.m create mode 100644 tests/test_dsp_edit.py diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..d61eb01 100644 --- a/docs/pyttadoc/functions.rst +++ b/docs/pyttadoc/functions.rst @@ -26,4 +26,12 @@ Functions .. autofunction:: pytta.plot_bars .. autofunction:: pytta.plot_spectrogram .. autofunction:: pytta.plot_waterfall - +.. autofunction:: pytta.envelope +.. autofunction:: pytta.start_ir +.. autofunction:: pytta.time_shift +.. autofunction:: pytta.append +.. autofunction:: pytta.normalize +.. autofunction:: pytta.amplify_to +.. autofunction:: pytta.zero_phase +.. autofunction:: pytta.randomize_phase +.. autofunction:: pytta.quantize diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..d2927b8 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -113,7 +113,9 @@ list_devices, print_devices, get_device_from_user,\ fft_convolve, find_delay, corr_coef, resample,\ peak_time, save, load, plot_time, plot_time_dB, plot_freq,\ - plot_bars, plot_spectrogram, plot_waterfall, SPL + plot_bars, plot_spectrogram, plot_waterfall, SPL,\ + envelope, start_ir, time_shift, append, normalize, amplify_to,\ + zero_phase, randomize_phase, quantize # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +148,15 @@ 'resample', 'corr_coef', 'peak_time', + 'envelope', + 'start_ir', + 'time_shift', + 'append', + 'normalize', + 'amplify_to', + 'zero_phase', + 'randomize_phase', + 'quantize', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..f318da7 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -23,6 +23,15 @@ >>> pytta.plot_bars(signalObj1, signalObj2, ..., signalObjN) >>> pytta.save(fileName, obj1, ..., objN) >>> pytta.load(fileName) + >>> pytta.envelope(signalObj) + >>> pytta.start_ir(signalObj) + >>> pytta.time_shift(signalObj, shift) + >>> pytta.append(signalObj1, signalObj2) + >>> pytta.normalize(signalObj) + >>> pytta.amplify_to(signalObj, level) + >>> pytta.zero_phase(signalObj) + >>> pytta.randomize_phase(signalObj) + >>> pytta.quantize(signalObj, bits) For further information, check the function specific documentation. @@ -274,6 +283,190 @@ def peak_time(signal): else: return peaks_time[0] +def _new_signal(signal, data, domain='time'): + """SignalObj like signal (sampling rate, channels) holding data.""" + newSignal = SignalObj(data, domain, samplingRate=signal.samplingRate, + signalType=signal.signalType) + newSignal.channels = cp.deepcopy(signal.channels) + return newSignal + + +def envelope(signal): + """ + Envelope of each channel, the magnitude of the analytic signal + (Hilbert transform), computed on twice the signal length to keep the + circular wrap out of it, as ita_envelope. + + >>> env = pytta.envelope(signalObj) + """ + x = np.asarray(signal.timeSignal, dtype=float) + n = x.shape[0] + return _new_signal(signal, np.abs(ss.hilbert(x, N=2 * n, axis=0)[:n])) + + +def start_ir(signal, threshold=20.): + """ + Start of each channel of an impulse response, ISO 3382-1:2009, A.3: the + last sample before the peak that lies more than threshold dB below the + peak of the squared response, as ita_start_IR. + + If a sample 6 dB above that level comes earlier, the start moves to the + last sample below the threshold before it (an oscillating onset). A + channel whose peak is less than 20 dB above its noise (the mean of the + last 10 %) or lies in that last 10 % gets 0. + + Returns + ------- + numpy.ndarray + 0-based start sample of each channel. + + >>> pytta.start_ir(ir) + """ + x2 = np.asarray(signal.timeSignal, dtype=float) ** 2 + n = x2.shape[0] + noise = np.mean(x2[int(round(0.9 * n)) - 1:], axis=0) + peak, at = np.max(x2, axis=0), np.argmax(x2, axis=0) + thr = -abs(threshold) + start = np.zeros(x2.shape[1], dtype=int) + for c in range(x2.shape[1]): + if peak[c] < 100 * noise[c] or at[c] + 1 > round(0.9 * n) or at[c] == 0: + continue + with np.errstate(divide='ignore'): + level = 10 * np.log10(x2[:at[c] + 1, c]) - 10 * np.log10(peak[c]) + below = np.flatnonzero(level < thr) + start[c] = below[-1] if below.size else 0 + above = np.flatnonzero(level[:start[c] + 1] > thr + 6) + if above.size: + before = np.flatnonzero(level[:above[0] + 1] < thr) + start[c] = before[-1] if before.size else 0 + return start + + +def time_shift(signal, shift=None, unit='time'): + """ + Cyclic shift of each channel, as ita_time_shift. + + Parameters + ---------- + signal : SignalObj + shift : float, array_like or None + Shift in seconds (unit='time') or samples (unit='samples'), one for + all channels or one per channel; a fraction of a sample is applied + as a linear phase. None moves the earliest start of the impulse + response (start_ir) to the first sample. + unit : str + 'time' (default) or 'samples'. + + >>> shifted = pytta.time_shift(ir, 1e-3) + """ + x = np.asarray(signal.timeSignal, dtype=float) + fs = signal.samplingRate + if shift is None: + samples = np.full(x.shape[1], -np.min(start_ir(signal)), dtype=float) + elif unit == 'samples': + samples = np.broadcast_to(np.asarray(shift, dtype=float), (x.shape[1],)) + elif unit == 'time': + samples = np.broadcast_to(np.round(np.asarray(shift, dtype=float) * fs), (x.shape[1],)) + else: + raise ValueError("unit must be 'time' or 'samples'.") + out = np.empty_like(x) + n = x.shape[0] + for c in range(x.shape[1]): + k = samples[c] + if k == round(k): + out[:, c] = np.roll(x[:, c], int(k)) + else: + f = np.fft.rfftfreq(n) + out[:, c] = np.fft.irfft(np.fft.rfft(x[:, c]) * np.exp(-2j * np.pi * f * k), n) + return _new_signal(signal, out) + + +def append(signal1, signal2): + """ + signal2 after signal1 in time, as ita_append; both need the same + sampling rate and number of channels. + + >>> longer = pytta.append(signalObj1, signalObj2) + """ + if signal1.samplingRate != signal2.samplingRate: + raise ValueError('the sampling rates differ.') + if signal1.numChannels != signal2.numChannels: + raise ValueError('the numbers of channels differ.') + return _new_signal(signal1, np.vstack([np.asarray(signal1.timeSignal, dtype=float), + np.asarray(signal2.timeSignal, dtype=float)])) + + +def normalize(signal, domain='time', per_channel=False): + """ + Signal scaled so that its largest magnitude, in time or in the + spectrum, is 1, as ita_normalize_dat and ita_normalize_spk. + + Returns + ------- + (SignalObj, numpy.ndarray) + The scaled signal and the gain removed (one value, or one per + channel with per_channel=True). + + >>> normed, peak = pytta.normalize(signalObj) + """ + x = np.asarray(signal.timeSignal, dtype=float) + if domain == 'time': + mag = np.abs(x) + elif domain == 'freq': + mag = np.abs(np.fft.rfft(x, axis=0)) + else: + raise ValueError("domain must be 'time' or 'freq'.") + gain = np.max(mag, axis=0) if per_channel else np.max(mag) + return _new_signal(signal, x / gain), np.atleast_1d(gain) + + +def amplify_to(signal, level): + """ + Signal scaled so that the mean RMS of its channels is level dB, re + 20 uPa for channels in Pa and re 1 otherwise, as ita_amplify_to. + + >>> at_70dB = pytta.amplify_to(signalObj, 70) + """ + x = np.asarray(signal.timeSignal, dtype=float) + ref = 2e-5 if signal.channels[signal.channels.mapping[0]].unit == 'Pa' else 1. + rms = np.mean(np.sqrt(np.mean(x ** 2, axis=0))) + return _new_signal(signal, x * 10 ** ((level - 20 * np.log10(rms / ref)) / 20)) + + +def zero_phase(signal): + """Signal with the magnitude spectrum of signal and zero phase (ita_zerophase).""" + x = np.asarray(signal.timeSignal, dtype=float) + return _new_signal(signal, np.fft.irfft(np.abs(np.fft.rfft(x, axis=0)), x.shape[0], axis=0)) + + +def randomize_phase(signal, seed=None): + """ + Signal with the magnitude spectrum of signal and a uniformly random + phase (ita_randomize_phase); DC and Nyquist keep their sign so that the + result is real. + """ + x = np.asarray(signal.timeSignal, dtype=float) + X = np.fft.rfft(x, axis=0) + phase = np.random.default_rng(seed).uniform(0, 2 * np.pi, X.shape) + phase[0] = np.angle(X[0]) + if x.shape[0] % 2 == 0: + phase[-1] = np.angle(X[-1]) + return _new_signal(signal, np.fft.irfft(np.abs(X) * np.exp(1j * phase), x.shape[0], axis=0)) + + +def quantize(signal, bits=24): + """ + Signal rounded to 2^(bits - 1) + 1 intervals across its peak, halves + away from zero, as ita_quantize. + """ + x = np.asarray(signal.timeSignal, dtype=float) + peak = np.max(np.abs(x)) + steps = 2 ** (bits - 1) + 1 + v = x / peak * steps / 2 + v = np.sign(v) * np.floor(np.abs(v) + 0.5) # halves away from zero + return _new_signal(signal, v / steps * 2 * peak) + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/golden/dsp_edit_ita.mat b/tests/golden/dsp_edit_ita.mat new file mode 100644 index 0000000000000000000000000000000000000000..329c698de6b9a3439ec7a55a1bbe754d1d978194 GIT binary patch literal 32755 zcma&rV~-{b6D8nj+qP}nwr$&XPun)9ZQFL=ZQHip`@Ea%_dQ7^m6Q4fb;${-$_fb+ zve7dV$_c5`nOoVK(Gn`!8o8P~INNdm*NF(Rv(gfZIGY){nwb(h*mD!AxY-lRJGc`P zG7&N|arC z$Cge^hFWJs=$ssr3{;bIy(8byl@gt9Hr=}#*OC~seH;I>nwBYfew6r1# zRhC;hrpzvvfONocI80B4(>pbXA;q;*l{@fCmcvgoblnsZEG&>8Q-YgKY<5jcmEJ0e z$*3?H+hlwLpYqTzBm3wEh3vrexGH^Th3Ik@nK98fqS^^_2`+_MM0uje0rlKsmqU;A zxYZOJB4C%(WOG>Pg=3($q1o4%I!~cEsc1RYGh9Xip5$+v*};R{PoB*fT8p{pP1c`akuv9}IU zEnnC>fYQ)S2gnK^jQ1mNw&wKxsb!{YkUyUw&x;awoTP&t<9et=_9G!!LK}Vz;G&u8 zA0m?-$5_FN)0SgPt~+Zy982r_+_+YkGV;Qg35%B+IoSidQa5S%*$!Jxyr8w0K2u!E zN}g7n_{o_{16R>~mNr>diE19NP|(Aesc@N8LCY0mcVx5?M{!@u9cjOd($JGP9w-%= z^N}Lx=L_uOL+beB&uf%+k2@E6`u@eqO?!BEA5n6ep?P#@%*jL=W}r4tdy4~TidM?G zMH^|(!wYArE0*eno7`fmh1M z9aJepBv*q-7NE-rVi7Xc1WVjPlC1Y4ob3afvqTY0G5@2Y{-$Dd5HHnDu+Q#~Y4N$> zbVP`>jlld#yqwedI_W!M?bUc8ewSN744 zqEv>okBV&+X>w^bqc9Z8Z^sb8I4IRbf8R+zmpBG2@mN1QCzo64x*bR#S7P`zZ) znFQ?-9^BEIenvz^tF*&X9zZz6Cnha3iH%oYF;&?zZqK;dXs3=}Ss)-ecy}wABPAD^ zIXs-f40`5_GIrch*=6%hj?_H#vIT|dMIQFdPE^U31bk2y?~$~R;Rkl^XU+0^ zH-0)N-2@5=hn^p?4Ffk%oN1<5aVnn+1?%Z%B%g}UQt`X3-eBQAh^23%2D0g8dpAS#gS+Mww8n2Pcya#CTxzqC7mnTvRO)V3 zPr_5Ek*`+H^@)c*`>ADY%7At8n!$rpv-~lW?xkX?7Kv0uAoIUf6Tayoe(6JnTBjOj z7c2c%B)rZ79vy(TZBsW4uLZ6odk9$_Uo@th6{qmJ%AaJ?AdNS!!DV}P=8a1?IfUm* zdYvG0pb7^lbCCyNBwA5R>@gvM;}P=@N1b8O4SM88=p2gSHlQfI-^6d5ksygM|E!&R zpe15|jNMd94j)_yxR2xpTs?)nFoPk{pAh4+iP;*E9yInwaOsFX;ZPO!g-5Wlh;HqW zK2bc`bg74I5+T{5Hq@(&7S6#3kRM!1Ax+60X)YMRqj+#B8z0|UCz|9p@?v2>Y;DlR zUx0gUIE+6&l=HHaLswtd^@qYigF>zDYc>G?_NWjZW&`)ypm8-RBaH>1Ci)bPLurt~ zgTkyn6(J9>PmYIEAKs~Dwv1A{rd1^P$bwA(tJznpjC_YudtX><-l&pP%ExN_R@__F zaHld|WIUe>Kflqx{ceXkkJMI=-%7`)>m}}12?>XTYv8hRy~L9zA#I=eSK;JtXiILk zKz>(-dQ`Ie+}ITL74Q&*I+b!XLU?3(Gpd$QkgM@!{#H=DsBuq!REW5{Fl;=#t^3KR zuYmO^1#D#(%zJA_d^mo267*ud91m9Y%H}M^UlRL$)Gm&xZ3nh8f*3c$%~x@VviC=H zD8zDd6$#*hgJ@!~2cFx++p2+XpR4Oehq29-C{U+k?12;W8RROTTZ{Svz>?PxGTfNpNX8jDq(428B05 zQG!G}sQ^m0VrsScbn_w4%^G|ft|+$44f{LgO6zLDsP8}?dd-5An=5uD4q?nbebR7( zE8RD}ukx+S>vu}jH<&nX{|Kd`(_eW3)*gl7JOhAw|v}cWvDE|sv#JP zJ;P-0ixx7np&czVnaIK$Yc1r2oj%k;!ckzj2KC36HJm>M9d(c~eE^6Srja)^6Gi4Y zg-gZr$(pp8?zEG!yy-LR*{ZxHXAe39R2$h#4(1c38LC(ZFBTvbxvdoz&NUMU zj#04QAl+trf@IOlO_gr)N8WX`x^*)E!QHm-xZV<{2~Ygtn=c)l zD@V@P8uZ6JkC^mRCNXOKr3LSzkZ$X*Vqmt&7tKYZp=Z!;*$dgiJ27b1oeRW2-UcBj zQO0pcXe^QypnL^4(e$(hdgUM{rKtysx>0~CSO4^Vwi?`?MrUiDsYr~A z$>|N~uOC%MXFD?U_tYGf8sl6)w8ky!qg#!6{>H2+{AFDG3iAqKIH9F_rh>2~7Pc#v zS}_MJx}po4^{>`6CzZ=LZ0+z|MTlSYivxho<@vno?3l*0Zr+;U%nEpQF`TUDazAdg zxbh0CGg(dIItTVU{cZ^0iQL}7?m*84W$f2#(9@Eo2%W=)uqpGQ3Xo#iy2sZQYQj&a ziD$RJyol}&Yj^w;qYgW$?L=j1wzREY&8_3aZddf$%`R_ib?D=clNf7v<+m&DVa7y5 z9_iH~OB*A=Bo@T7gda(@l7eVSPvjaKQq@JKOOu$eNjY#`5wh@#9CjwjWkg&uuP*&J z85L(hUBF2-JAQVrDt<$UC9?w^hP+YM2h|5fgxpe8 zw=n}T;tu{~Ng+6HC7gFBGpd^EJ_EfefG&LRFrX_OFt%xCnl@ffl`#L95aVUUtPivE z*`zia55~NNF*%n1TtZ2kPV7K2 z0uYqm!+?_dm26j}A|jR%o>Pa4){v5HO(A4Sy1P|ilEIIf=cZbtlaci<&p4ME?#8$o zkWynJbJ#pXh9cxi+s0j@i zZB3xcuvZT^Y)dkC2WNeBDr6M@U~PJ{M)`!!hIF(te(Wbt+jTgyHzz?zXM51ClZZc@f+g!cm3}}zjdJJvYDvMjSu$h;uE^+at!#K8n z*Vi0^Tg8s^K+1!VV!P3)RJzcL-F>IlZiNiP(H1*ZUOvb)Q?j+Ps_H7)?-(X?o2f*WOYlFZ1}AJXh(M=|Ub z)%lBvHXM(*k`h!T8?l!ItcZguLx&>RF&8~#+2AkhG|bWY;ve@~s6ELt0T%eT!^O?~ zLR^8o`1?b}LWadEgol%PAoafcg09Ft&Y|y|j9aA#*kWwYxR-DYiiC}Z;(gaR7brn);6;)A_tGj zoe@>8xQrH~4ODJHG}EjWX-+mu3aRGEtGYHz=DLB)vNsp7 zc?+-G>TVeHO-o&5o}|SHlB}@kG?zbPvuc%cA}e-yOGOEnHTd4~GCiF;&2>*o+$n!x z`Wv}c8&1~iUR|3jIJeu%`imdm0r1?y*t{bPz%7c4AL)0bT{Mp~5{_hZo~=&Ewj60# zo0shph)RMgn!9Tqe%_b7Xcl+2O$RG$6#w@I$gwFsdD>*pQ79Q||23t2Q5Tu;z!4sJ zAFxG(BR4CH@lnaik&{Zyf^^A-l1d!q^k6&bTKx0VyFtL0`!gJl z*>IXfR7D5dfWTxpy!)ZsvPwbUO3=jm>EthF(As?=a@EXFxivrJ6NYcb;o~w?!V!OS zesEC`eUq#&dw(fbXA!jT>!x^iCUw-BH{H1t{07EX&%d4vvl4Vnx>25GRpg=@jf1Xz z)w3DzBa~UqBaJ#r zNwiZr9H}qO8lYOtB$eoagNkHBS?DYVR7MCk8jpfYQr7=X?Ug1rVG_t&6+U(NDIWJ2 zf0SU8QlwCUO93_p_ncT{(uNCxcYLQ$IOZb4RSV-J=9u!qH=q$pY)vAz2DO9DDdd=* z8DU`Xv$p;FGDc5y8^7JQ7;yDa_`3$+xgWE4$ldEI3%GIx8{lU)lPknFE zQ%pMtZK77?sK&6PCpi@p``_>og)}V+?kbs3X3-#WiJo-^8{~TzHc->@IYz1cjf9)Kmb;Oa8HTgJJJCuX4GMuriER;hY0#l^x)t+t zB)JW4m(f=}j>>hjrk_wWA&H&;k+sWZHamNE7t6V>M-0G+k+}?Fq|X}>e<69MjF0QA zZfLrVJ4W^@G&+txx2Cf@Z#R1_z$ZDeSVHrc|EhN*{DjMI;-`+mS?AM9S0rEF!0*Ef z%aMN}&<8tA(;U08HXUiePP`gxk>J$hTkDqX&-MMNZLeL7a3@r!i=31SzT|PoUrJB- z!kw9z;%swW0QNMXKt%J&7xIuV{Ge6`lGpu*#@Z8vBYLBAz zh)kByF%OTRZ*AWT+jNtnan%a>uXXbJqX{$E*cF!^PiS#_h3pPbY!gnY39CQ27m0=P zy%jl`f<kjA|3r1mNEq!0d3Dfh#*fW-cLu7c1PFtIQ?vws>277Ip0 z_pkBMC&Pf23l$4$PD7H5Q9lZXwLYzh_oMSgmipbM4OwlrI+OdJ^~KR@f`>I6UW1L4 z?rl`jIk&r;PW(&ivu=C)-$_;u@Kc9kn>%uRF$Z(6dY&{P9Zbvyw631bHG*cmqGKLO z=VwBTGx_NlJ}A!%sacjT$>?psl$uXDURV0$OwX3@@km0YA8hqK`C5R05)*&-J#U~+ zztEPSZ$}`(1I2f{)n)vHRKOA)$lD^OpyjS$XDMjmCrLs(Kq#gIRWw135u0m3HVMI~ zP6N`(L%99eH{$a;T$7g%b@fD17D}Xy4wh}xF{~~8afCZ)H3_m5KCZZ#H_0YTNfrv_ zzzdZRUBYA?6VA_A5h{Y{O$fFStN~MX#M#FZSEKOotr#yp$bpyW0gBXz!tj9)NTE)Y zy%*hB2bslzv|;p|5yl$+(cX}`vUj4|KQZR#2-O=;%BWM9T3tp3%a<5!mfwmuQtogy zoArsktm$2vTiS&MGR5p4CwVf@a_02vtZ-TuSxh=DG+!eXZz`ZnJHVVd^4Imw#qr=vES6oBO@|W;*=z{b-hkP(oX9JVH1H1^Ql@p4ia~ zsrEG5vqB}N!-ENb!MO7><=>g|;AS^+%1-*7p($0RQHkVC+Jc`Z78lQyV~57*{{Atn z0@c`;4A*8Z4JJ|5-amUX3qKk!_g3VweKpGt4W>s=w8~#Kr`g^toSo$ut4i}hi~@Z; z-Jy$Iww#H&{jGTHX<8el_==H3Z908)x3jlpC6B5&H}q%H+q%;UfH$l+Ykw)Ytul6m z_fo%8i+G;A)zM!ck}{PL&SlBHDsk~l)_vxhNUTJ;;bo3{vEEGgOY;YjEG(DKOyC-H zQ+hgCIs1a$Wb_-$7UQ<&)JGrf#HHp1-_F$P_?m7dtD^1p*j|)P%*KDQ%(v|P;cX2i zdiC2i?8V&1z++M8O$yVK$V{IE(H%t$MAm5`aU#h2CsL4oG)4dzsQ6Bq%te_=1CnA) z9cd~YhLT@ZD%r%wMKS3UWGb=G647_*iJF8)Wopiw8dNzstfA3KMqL`rn}0EZ%vgFj zW}Lg}iL`67*TFHfc+9#6kTacBChdDQaVCQS)S-6PrXi&hR7D+W?~1J zYYS#p(6O2Yo6`&_9t-l8Lpu`MSqBPlTb@($R+X)~T=Z-A1%oM0N2z+_)EEyU=9%ln z{aZp3;Ru6ZQ(7FlqG?6Pf>2nH0SA$4Gg#5Epjg9+^PtaZGuTGt#;h3CMA#WgCd%Z4 zb`fbPWqOa)#Dck2^&hI3OA@yvyqs326G!f zBm#EyhHoz-lS>$>Y`usQE25+>Fxfh4!=C|Q{0;Q_It`#9?8kM~EQM0*UPX@6duvU6 z7=_~`7S4394z3WJQ%1v8B=GN5q&ogGxd-r zGZpN&kr6K}l-XgmBDfDhq2LeVCqzg>7hA^py3lkY_yd43;Yta}w|!qhfGprE?+uA- zHl6|I+fn(^q&$pZ9om07_1*cVhoM2dw_M@rxGx(5GM>MGdEr}gvU^?WXk?b&C$p$8 zyKyGha@c#1*^P(&e$x}TE0ig`?H}`UZtU20gXTfpuZTyl!lv+ECmpNY&SYLOo$Yo# zS{=&`?bo%u+BCU%1c&qIfH(S%vsv>zb|<7ej$Lz3R7?r4Z&L^dMh+#cc8Kzbu<#2V z@zor$Z*)4M;`#K{{?yy+xb(w<26DFv-7m?4Ywxi!5nY-#i_Qe}2mHm_EjxBdm@0i|(M{KrIfJ+0?rJb$M8A~*l%ZxsJCxST^1 zLqoUyq$K5arcSgJmP$Y*@ydqmj86!r!ypE_@D+Ex$(DftB*zN;8*tK_xAUa{2L|-SsD3n!EgSr zMxP%M@tIv4vRq;-1}A(T+BC5R<#NzS;Cr<^}V>bf}xvP zCD^@llb{TP)W4;svJNhPkH~}mL<5x0GlZM@*HK;~J620W%TX3!d zcd^d){=2`9+pyzeba=$NPx8<>*2;HSw;IthpDZJr)0@K%tgSFpUv2j<1TZ5oZ7Yzu z44qN=nCN#+n{*@LDdLg$Rp+9?P;#O@t|$SVRwQWwU24^NMRfW*ITYgT<`GSY+@S&= zn<16R_Xnk4^h~l5sEeD^A?9iY!W75wPVTm_^JTC3Z|gPBZ=)Zt!%Z>lL~r&KH%@%y zz6R-JTPf~pmL8pQ(?@Jhg@3*GMr zRO|w%l)izxhljQ%BN3xuKTmUiPL?yNblr<{2r}A9{og#x>Pz1;JfG?~JajiDEIDh# zekn-CRzILal}}Dbcmf00%#SB-sR4*#H6eqF!HdQ^1j5yBhis=`XEaD1b*TYt+LHmX zywt88xHPEPosSFfJfVe8;y=a#4bLb`xta_n5+kn!c@!SC`g#}{%lqs$K<;NANa z%wx3>8Qj9qI=Gp;&a$sz01EVh@Q{|VR5+_YI-A0o5h^MF{3ihhZ5~a|?NjjB&0C}y z!g!R6dV6F%>;-h00~0cW!74J$&jpkZd2?t$F9X@n5;ru4p&m}dg@Aq`Pru5xOs`6KB zt!;17cKj-b>hZw?vQsVyTsM#^#uOI%{XY%0FPanFz9H(=z1dl!w4Qo2>h8L+;GN^W zk@wArxFDS2E>NIE4~;xYUt7s(&a0ur6~1B7DO=%BGiZT4?J=Gj4#e!M*V9!qsNptX z|A)mBtf=`s;h$?32LFyaSR;X^%<{ha7{&SD$>?RYkhWKJq!;)XS($LY`(I;6Q@{%5W@Z!i; zE))F5`_-8G4i()Qzc-QLXbu$G*^^N5^Uyz*qq4+DLlT31lvx_r8-9_(vdClkL=C}oFF3q-}|1NGZpTiPW`mkfcJ=Oy!UV(Q~c>{O1B3B z{|)4ojvi4JGZF&zd@)LznZbrjkM|}88&A9u>8r~vV}6GD2!RgB2Xu^!h~8mw_RBV` zfRjPDLB!(-=+u$nf_#hjhSPOXIg&}P$`3>wURkKaD}=WXi_3^!KO)haP7hmumQi_v znIj}2zVt52P(>&AlTG%2pe9!`#VOv8*cUXAWrgYOvOi`AJH5U){^W#LM}rZiq3InI z9O@6MTmU+kvjGvh*taN&@CcVr)VHNW^kKtIjGIk&T=8~1ikaCJBO3U$y~NQP%h!(9mf)2bY@=jt!^c+4ac|E%yXb`8*5PztSiqX6k3v#*9|d0F8Of!yLSIs^d0$Qt3dv`>F!=~hhE!dknedwKkJ~O zc!)4OeGn6$jf`fE16gx>(Y*x@etY5Bf1?T<@-&QRPx!=iD)#m`i=KgkVdD0#lctxm z^MkPEOn_tj`i%A!jdt2mH;CNr77PB{stGhk$z^Yp4{@7l8`*~wb-1NR<~i)ehh$PG zy#%*om>llHxjpmzRk9D>onw2#y?P|H+fwUUrHsJ;8cf^=!Li?l@BMdEF1$M?lYAAO z=IC_yv%hBI2jIua6xy$RaUV>v$+_ePUF#(gH^1|5_LI}gB^~_2sm4wE>qCrgE2EHt zgI)GckN?&2Hq(1+7}~;jlgaLCbJ@r9HOP?nHH~_;%^>e~lV@C;R(F)ASM`*yOZ;ww zq3z#F7<^8ETATy_zu<ra5lw+Ex8bOLF*xlbr|oo&^oi5E z`|t#=oMCsjC--pH#xTcZoVol<-_HI6=RIYgtBx(3_5SDCj}TR;0_D=FEX(Tr@LW9(iCwnKSw z`GH5%yawraXwP+7XdnHnkZ}6$P_&r@BeK4V`sFr~_Hh>}-FRbA4yprU|9&T27O&i} zIx5-$cz(yDtDvg(#iq-Da>iIc*0CsgrAZThSt?WJn?F(1I50Jg@Gt)IIeV!0XK%JHYsyUz__j8UuH`6Ki|6+m{w$KL@${*#bqEA=X67L- zDTpLy(_9A6o@zMjJ8X5_cqlVAtsrn%UYhT&W;yo23ENl5{o#205wHG7m$qG69g`jZ zZ606U30#z9+ff!y-x$Y*-adr}%JL{M@Ky+Vb3*`NPolIyNds_>*%&5HkvRs3@7TUV zSgn07AtEfjj$VaxrLh}CN%RJQZSPDSa^Q?L(DhiN zK3lBy=4DT-DvlPGU>>eK2T3^)Tl#sx+)SaX$uH`m3yk<1#oJ38U&ksp3`mjN#{9Ul zYt3e$eASxSvo;UicRQEbWWHZ8W7hTC5_te;j947l^`kH5TW8>)1m1ci1#0%jz&jTg@tt zXorm~;y0~+CKq)0LZkQoHK(QoE1&ZP((%l2vii82twgb&(9_jm^(!9Zu&i6{A_m#S3Th(rg;H@*N_*1ohQFEL1(sb& z7+UhfKg56A4K(&u)gWJqmdmLLpbm+yzvine7#UF1-u$Lk=jV@C$~l=Khi$ro-hjbL z;+){Ne;pQu;;KIA=5cILM!O8YG>JB?ZPNCe$k;t|qO5iq&0^0TBY*Z1VzCix|1kz9 ze-#Jc^s;~}^?gAthpOrSIEXcJSoc?dF@pDb7XE`?=$lmo4;9iHmM8@6UM$KM133c#ms_oY7KSZ0w#QY+r0J)l^p-(5g<;mfzyB_iv_}xt9e>wS1*2z zoB&zc>`|fjiLr zX=}=muRS4MVrO1ubF^OZz37|n!t zbrPq+Xk0|r#4d|${dE@T!_oV%>j-ykgYVSPdr4Q;qd5jjJz*G^ci55i3$r`FJ5BG& zoX>}n15%dEqhPZaV1j|OVMMB(;2=oi2u>vs>UZnAFs=C4D=Vw(*Uk2fvO1D!Wi2$8I5 zm?jv#${WI62;Oc|4u)O!ZA7l-ri{HkJKq3B)JdWg`OW2u#J(dF(u<7CWucAZqLXf+ zP9|%j_xe_dRRVN3pLl$S`;8Kk+btx3^f53@(0f+l14miNIou8*W#nE}vgY?fK4&NO z$2R$_FhAaN=wva!U?MA+(MuhE%ab>k>h}|s?JaL|8c!v}WOA1OcXYnR)Wqp`1+L2h;hQYO?d4y68?l@_sjc4Dv^& zj3B}~1E(bF4nqz6S0@FIMQ2MCPr73o|8hw24m+sorcAj|dV8#L*#@v!dm7=9N7{eIy*}he) z;+-!3H_I9Dw(CWrdle!g%MdgtN7jJj>C}<~w9x^EdqMmT_|Q`=oNR9`^W&q$pYQU~ zznOE2;g?kwf@sW1+Q2VD-J=-TwwlmY57u*!NECY&4lH{m0UIdawYl9%c)%Fq7yT=- z*yE2&%X8D4QC80_MGU@a;u;q8f$ZV2M(aPC+!wPZcf;;M%p;>`?+=tD-^?h@K1)Rq zQy+7gtIGuoZ;;YxD=I(trLsAvV)aBY3sG(ZG2FT zKvsw#w^H1nKr@nGMzVkd{=4+oT)a>W`8CnQ;JRWCfG|J&(Pi0m$m9)=BC?*p8`kUT z1|AcVFO<#kRVu(e`Y$T}@Vx3|yScO#%Wu~YDQmSI|4lB{%o#zUOlOg6g8VK z(oL%f6deSjv|Evg15k}|Mtz3sk*ve`6?4b;kL&cF5$3=R5RB*XnxAy)C}fX6{7C)2oW5rhTVN;I@qZ;q#D5Fyax^GLw4`M=J0!SNN0qBvSr)b7=Yc zL-!P6%kgfmj8?gqFE2nBRe1cBY#Di^&0VgV9fzN*00%o)K8IV9tZyJVwHg;C*aNbS zv?PeulNoMKcU7XB(x!xM7hd6E;$L+N)5ZxB?y>Sj(5J|Hp%2BV#wG{@=RCK@cW!vHV0MP;a_WyX z18oWI7Lc}LYApZv7ePz?IjQxk0K&1N)&|Q59^31%cuJ~M0htb9ep|+TG47Dr^=jl* z@KO`;!yq5Ilh6*M((p7sFIGV`PYMsT(+;7c>~^_gB;L*vh8u~~k1zCQrzOaK?&42f zDJT%&uW|h)k)RK3v0~!|4IAcbf^F&w2-L;G(@e^DMAW;;!uJCvq%2gi1;2Wcn#ad{ zV?m(4tt%L&M7$zAmQ1-|g`g@qc{KKqp0HS|0q{hq@)P>;_1UVu-M_DKV72-guw-42 zIdG2}#U^V5Z<`&udU{s4d8c*wjat`0I|IK8&$`qnjyBi&jTy+??gQXEUe@47&&rsX zFK4jzVok!+3wehIm13Gx<7T%-x=Cj@B(?*7VYDkoiC1^mpoCobX0bFNTs_MOb@-oLGmqcjBgw}r#lhl7awIIOh<-z|FMsC@B zB8}%efWIDlQVa9q*;V)e`}#izjeV*I_Tz~K_GH7*+U+P|ai~_rd$j?nZ$aCiRl;nk zlPQXu+ya%|Cqk#>${KuRf{l{72iY)91s+{#5WaE9*gK?wzVf#LWojAq*BOCxBhmp2 z&wEjCz2KGE+_cvr16(r`r*z zP54%^)jpsc*AD9GgjPbaja2@~*OKCR(KkPZU{|a~4+?CPtxd77TAH#Y2)44RC`)$M zA57a}Jn$C5wp=BVK~kdcgW$vCt%KFBE|J?(&AZ)mJ=V6fzR1n5l z&u>M2@AV7X;_^$Q{52Z(Q9GRx7I^HyNY{xfVolqz$Mtl+K4t;q4YZ_X3S1?b9XF6W;DudtPP}8&Nj3 zbU}6N@CET~Dh2gzYXqWbgc%SBQSDZTLUtxhz^&$&!RrYNV!5^{mpDxxi06T!Tu10k zDK%x!V5pPD!gU1ij7|}}iR9b!zr&i8x>!?wyMbL_(~BLCBJMX-1S+h&{R0dhJb*Fl zhmM|%p5((SOUEC|kuZ#^Lw!;QlyEHo0{apL>Qe@vxE%me&Nq8Q*ySr%o@QuvS1%ar znn_3AAmbD<-Lx#oFaz3r8v#C@zAAZ7!tn2TI}A21N3hbK~|Y0jQ7KLyaze?M#}v=;ERK2c~NG8oGmtDDEdmbq$0I*;P9#QLIJqZ z@!|yVf#}r_LMC74AoO(2z#r#m8G1QV4_)+-=U*Zf+LDURZ3fY6gC5aunzCREo}}^5 z-3_UQ9w0AWRneF1Ym{zjOGGtv*+TTPxTLjB^TzeHblo#y;v4wuVvS+EMqe$4cO`GO z2QLd;tv2SAuPkNZQRouN6|ZUj&mu@zoyAgj?{_Nm)-C4DEsx0WcQfsKIo=otOSo}QXKW<- zaih?z1RBH-r3j*lt9672j>zZR-#dd-T?w*{+EKo!q9hkiu|bjF24LW?>@XW}mlM8^ z+z&g=E)hLrdZhfcpc6DK;x~DIv>!F8LgLkkL>T#&gdv1IB!I>@R^EWl!?YdwJ{c`(?d{q!l(NyRm2;)`s znGl%YQZ^m>>eQyI%!Cm)b+CPxO_slKogjdbRRYK>VVUe?{xzVIiqP&YP7_^5bf_Fa zFXUSv4e+EoZzCnc6G+eq2^e4Ly{9w@1tU9>U9Vm{umw)xUm$O;c$(FI7ac8fVIY3du{iPsGv)4*5{{InQxHf|xgvakGIb?fP4NqpXm304U ztJURQGo#xvOa1%niLvBcywAz+W(%`p2(}j=9CP=P0MgyIQ0L_z{Pg>e|g4fvtJcW$7=FF8o%gaO@$ zs*qrj7}~Ahi~8Jf0D$F84sokbJF-mwanvP}0o`L(bZ;u4Zmi|y5sa)2`A$d|MNu67 zFg+~#Uj9`mPo{?`4^$s~XAI*$Dh|Ld6Yu~kx_bFX{S)QRhiD3@@dwp#*%>w+Z)RXK zN1t_;t9o>3r|?bGiM9UGFm3N~5cBWt&Suc*$Z?RsOF9^WeiumDl>pRNH0%9*h_a!t zDb<`A!(`V3DW1V$8ocj1wBqX&czzA~kc00u zICBsYZ@|8}P~W|>p}q_7thhn2#S(3>2E*Ohp&Z9PQ8>`o;ehtQ^%1R*Ml|T&8Q4Lb z&H~sWcxqsH4{6@wm9da=_TSUz#5cip^M2VvU`jOSAcC-tK;?BbUxfFv4)P`Rc+NuL zD^H}b&kL%U=L|GiK8RnsDHepd3#F&$npf$AO}YS9!_4>_84x? z_ghKe1)<6!%4XAm%ryd_(4*>rR^Bgwp~1L#nSvJsW4s`4`^`c0{U@=$+=GGYH)Zr0 z^@9(-2xi(XjhSza8-NgaImf$5R0m{UZO2MlIoQ+ds0M33RPF?r^$- z?ss_L+n4g~0+a2b6tNre1-5jf`iy3PC}>cD0z#F$eMY3dptJxgA2}0U*qVU`eU&kv zv1%~#SJ8vHrpkeM(ItU{@#KQ~vtoJg)s%?1eIOG0gMl)~Dm!i(o$sgXzW?r>k@>H$ zV{BaGA6IO3f_<=xfsaaE_v*y!_TCDo&>tqKJill>`hDU_7>5*jImO`xoEP!e9zAFb z+r*LbuN#3RRfuj?{IDEy7o2s~xeye1_u^(_o4O<7bibO;{dp}p5$Cqr!GBfC!?Imx zAemphfj+e{_BCyA6Mp&B6a!e1$KU%zSZDdHEj`iykk;N0aky6juM<2v$(T|((zC-m z@)>}}Tu9vO_C@Sh#e{)Ao*wMc=pVQqs@V(+lIGQBb_dZ)-JSR>_w4b2%9`Xi4{@&I zC!M2`JLW6IyExyl=`NP3KUBb&_MRGUtmic6yo}`R`7gLsi*EkJL@zgA zVkh6zCx^evCWo(p3pfw`T+VOJOU}`wg1_*edj3264u0peJN)c1kS|>u#dY4i%b)MJ zfs@`<#_t&i_{!sQe50AOeCbqmcKd>G{+>S>e0AR<_OX07CL!a(Rg?tt7nFWu?;Vju zBR!t6c>x{#p|*BTljC!L?l|#nzbNp}q{#4Vl+VgJ5(~?pMGb*1gdi%{z>3u2{*nw zVv4P>b>?6A{DGaS3*pZjJjscq=y893N%I@JWVywu=eQ3N=lJ1N5fFc70Q#NsL>Z16;+C7+0t}fTpl8B`jzh&2HIPsFfaa6h~I)i3Nk&xzBK z*XE@B&S9@2#hCTRbL^?pf7z(j#+a@IiB(F#{55=l4fCC5_U>HAsjR(;83%6W@1DGb z-FR!l%`Mg99-djo-@i_tKigr&70O&^Z7RM$r70O@PlUqf(Z`YP#S%n|4sH5YUtil&| zd6XDu!tZ1IHxEO>n&oK6#8XzO$Pn{y-O6I08n9PhQfS_R4@l=+2IKc=GkeF-j2+HQ zV@uoLpu18C*5}4N&S`xMCU35VF0-St+ZHcTJFf!FS_#DBlEtx==FZslD2_>TJ;%xs zF|6NA0=w8M5!;+4%3s>@npNX@U?=5L*x;Oiz2^nRPSWNM7K4$i^3+oK~%RFCoojtyLKek8Y zJ{xc5jUAOw$JV7ZV3Q~AupK)uV-p`bv7ZBlY^=mK`oZS`M#Lc-BjIZH*ws694ZyH+ z7PVNIu{l$(w*r&jN3c$Q0qjP8JNmWGjD5fqVM*s|us0iH*eiLBnErzj?6X=IJ#a#d z?fe{oc3L;1+gqF19hnTu$1pl(p)c0yy$5sDRl#&hKBH3|nOMYXDWn={LIXdiVa*>d zU?)S4Q=26sS@6d#G-O7Mc2V$Rb4M)MsvB9Xxt%R`ve}+qZTa@)^o(fg$$zq!>F*8L-uF z7hpz+I9ALXN1HPHk%fuVz=Hw%0 z>iK^3*os8v{H;T%_EbA+XuJeFAAOy9pfkb{5;k5brcOMmo$4#mvvk43GY zR-u583iA+O%v!XbX6}c#qN3MoP(Ax%N~`e^=_0O&+AZ#ZBh5d;^{L+G#YU-Zl0X7t~XSoEfo4;q@D zk0e~b0oh3sOn|r;TQXS7D2B;m;+MUdwU`T9x;=~W-ZY48r-aOt4u3SF`Z;6RqR51f zPg89VTM=VzMaJ&?V|tNp1u9+C#7r{hkSE6s7{B@s8XgZs_wU$8Cs{8>l}ZmXm7~>Y z8Fdc*hzaRVXHmw@P>yLy%S8@;>1Gmwof&Jq9euQ?k=EaygZ7*>V-Cs~F(YO9NL%;| z=7K#RS~u3p2$wogC9~P6zlQ>?aI27UQ%i=eOFL0rxhTZF=sw!0;ep0~O``7~w_%!$ zGEf9{Vtx}j<#t@y4U1IX??O+ry6XeA#!Af&$&a~& z>K^`1I}|q|CY+Etbw`|*>CHf0GxC^?>!O+I1Br}!|9NE6b_41qQqNTWm1Iu-2GFgV zHQST1m>%q@^NS>sW$nT~tJk z+I$$Fh7z<>RhAmP?nDE1(R7-5DcY>I0X_65kvx8wk2sx4M=eyNXzQoPQPsuz%(Fg4 zuxm(;X&#=>?9ZwZes<476s^D0>AI_!fTt%>qeultUPF(G3|h=4Q#s+Bi{m|VW(V?9P@X`k2kW+z{^>x|JF%b@H-5VW)t+I z-hIpglT4CI-wbs1$QHX#>;Ry%RCE_)Ry9zM+*ByO7pUT@+Yi zgdUxLln%F3L>1y=kcfv0jJ0hCy)L~DaZAl6?*tu&JwH7{_=hSHnJ)vhXp9E3*I*Jp zjaAUgT@BER3N6OhwhcM8ZX8j`ltmHG7TR2X8}iSLKYjm0Cw;AeNBbY|M-&69=wjX| zoHUYxxTd^EN`GPW<7{!}sjUNY@3$NyCm5l9=K|>x9~nk9Mg-O0$Ra0L1*Ux86z!$u zN>j=c^kujb`FS~)KJ1f6J3GtJBhHe{+k#H|sHz;gV(V>Myk8U@75373;-!@0PG!1w zMw9`cmmt=q7D(QXJBZ57LFi}dRd_{TB9b&ROwZX#FjWngC=<X5+b)s*zmQ+jmXGq~l98P&eH zg-%)S4M(}?P%REUZKXt`7{mrzSsqEFF6k{dvX@((%3=!oWBc2wDalfr7@K5o%=LzhbhwJ z>`8?uwNSUqo*=(6%8^4{Kd6353t4U@z3eM)FKL<*$yq*sIHBERV$60l$E2M9deM60b_dcnfY+h!wWQ)oC8}05=d&sW#nPgO(gsC z3Zy2i0deZqrH?&Xj|iPe#Azal4zqtt+4Q`J}j@vx-w-!zdUabZCM>g7j~%9 zJ-3$7Ez3p7ldTV^rIP26E&tzTZGR(V*r9`*4)LYBr{>d_oA<-}RxF}Vt=W&LJ7?2~ zZ8uc@!VD4n_Zb{U{e%aju2LtQv}j5_9@*Wo19@$+lXhM+M@OeF)?!TQ`@xgs5|4wlhm0ES9p3@( z)YGE9oo3*XOWWx~FBxhmW+63g_#9r7{RUcm!-x)UQ-t%DZl&klaHF3c6hmaAj5K~V(BsQE6r+m_+vo9Tm1ANVdPMt@o!3Sc9PPii6zWX5+=_N_`jQ)mOzpbS6 z<^@wao<-Dvh&VjIAB1itBgmOz1GwNr0kt7`5Z2f<2|tRSfUHltz~giA)XEAY*tisd z?K%>uixs2rc#9r2^YjQLI<^WPy*WXZXlNt5#JXWk*bI4JW(7~@AA_@eohUzgC7mj) zAQC6^=?xcC;75+)^fJ;KR`qnCz7CV*-4z?DnE5!hFd_nua_uDXeM=FuE!SW}i(izz z+#n46TLuH*fpL_j(mnb>Ckyrg$G^*fT|X zUO7*RG?-C$v}>uhQ$H!FOrCyKnNIEawi>cHconodD*}gHYk}7-P^WZ_3*kqR6YzFj zSMu#*6~wyrJNbT2nm*(R!mn5nq;O716eNlvLsUL>R$QDm$}5LkAE&_@UH#N>t2L!x zl>#5SAd37IyF_jc6~alU7gA^Be#4I^wo^0y_S7IgNu>swfl4cVsU-L=RdD+z^fV8H zMcay?XtEnVJTDVc@Lo%uY5W1t{H}mqMn8h8?_NUR!Zj)5^Qy3_pb=IR$^E}aqsWkC z2EGo2QK?pq(4M%pl(rSm|&CL4Hz(`)Ev#C$l_xE6Y` ztPt)$t3Zu~rBSaUjp19){g9$;50#FARO#^ySgz7lxc+Ms(JiM&Z5H1Pr|V|J(w1cu zyJ?op$+iNev~N>C=bn(JE7IVcXMGU0yo}P)pCg(cOTg>h-h<8O&QKi-M~K;#6zs58 z8m@ek2S@+d1`i!qhn~N?3W+$}gBFfOgIOCc!u=c#uQJo2GPrv9WM(P+UMmpBe3M|+ z-1GnMd>wqXA{#F6(4^3f$KVEC4Jz9%1)hie0~_SWQ|ms55>FcasLK3yqJM1_{3Y!; zeD_~nNN_P9K5c&$-fyP|_YR+>)Fz`TD%2LHJ;tDH4I6TNdkJL!@*Pydx&qmUMWQGej~v{C~OJ{15d$UP+2Br8HSRVv`Q*ilGl)d11^ZHW*+YzrrN`@+_b zJ*Y=(;wZ)Phve70&ET+q0XdHmg@3-#g;L+ccCcwVkvj`-BJb?vm~{wa_x{i^P}r*CCoH zg`(Ky&rlK%0TC{zmsrK8&qztM+Mj7(EFnXaHX3L6gC$hpc7ax;;I$D}TjFDoxwUR{40iV@*~-~At0 zu5PHeJa^|0G^jUXS!OCtg^JWc_BT#~7iMphhF0x_kbVe7$K#}Vp^C6rstWoUp9_KM zQjnhVdQz1AOnQB71tWIafuD(NC==BtF20b2^c1zozB7~LnB*I>b=f7z0FHN(isP!$f{@GU+1Y56!*)2U7JYfpX8^A)fLipiP6GkjwW3B2!L; zSKPFlOpH|~#cbYz!6Ls2m*u(SaYaX{GRqCx@Gymx?C24yUXLcX#!Erx!!^m_X*G!0 zH7C5CX99Uz9|7+@2mp2U5ol4@|E#<{ONx4~BK>6+LfJPwp`vPU@VTT9Gy+{Be(%*F zw@i{yTTPqIbH)#*uufGz}=4%`KKk4K5jOAW-p`FWOC2F5K@ ztnY)ZaMLl(rc$cx0NZzS2@Dr+g&m`x~oJ3&2Tt|TO#L>$)Z0*Nd3WaT|E zsO{Rnglba>AqG7XNPg9U{t0^nI$hZW{n;Z8K_+6B*Lp$F$6|tbPc0z`aRX@g|N6(P zTM7P+R?sia4P2S{mDoH_hvdbqB9~}vgOvVmA#>JE5g$Z0lkHA6MEeV0qV6&c%B!q~ zK9rf0_SYoHd%XZO*4PNvOZDTH>x_u6*Sf)tshXrtZH6$XLWT_V(5hZP zN2Z8RsH0`U_gj`PfNWvh*e#*ZrasdOsydr+S3May*3W>&ju7pSYc_My= zI5{(al6Zi(5~qKNLK9!VfG_{?1|`d0fYw93;CZeVG@Gabs}1^z*Ty>F5~)ri4NL;v zv?Pc?(>}|uZRy}W1!rPesy)$ke~R!v@t7!l+(7))uOke;=7GsB5|%ZU@?iV;AYrh3 z8PWPL2(D~51oy!oEN4$Q5kpZlXt&@FVYg(SP_tMI47WW?L@Z7MPx$W?Cfh9q6M{TI zS*VzJx4DT}>U0V`SaBUZwOS2)*ZT{&E4vbuR9i-bD_;XINUIVzTRJQslQKj=aVk;w z`4*_J&k5ZFtB4?GlhCvB6fq^$2{edT6Rj@$EN6f10l&mF;WO97@KwPV2)WrbkZ8F@ zn0#FfzHvA}Y!1*C?(LgL4E(>_$gCr%8xn2VKd?Y(_%|6OeNnJ->1p9;Zz(u1Nm{=8 z6blA=$XHJ9x<+gi^C8HoWa2*n0096043}p-mTv=wMfS|f-jbDhU%1B;g;XT6H8e<6 ziUwtbY}q?2QfZgxzHpBxMT(TrK&7FCv=pj;^1h$XFXxBz`2CJ^eCIY>=i;;AT=xmv z%k*N~xzu?ZCxsB(=4s#N^QvmLelL1#G#^#k8l>vk%Jnmw_H zveh>c&Hqdg38n@XF3kit-9rI=e+h6rf0k_5DTbQ_&%zkPA@ZnSEaVTpLf$>;2DRi4 zARlE;gAWA8OE#$63j=(BnTDLK54{MjQ+sy(rU*44#;%rAwy&Kh8TKoKsL z1aPZPAy}$ELYgkP!?$A5*0+zMP{3dgkuDSfi~Ji9S6yG=z1;?8eyt-9Zc#$?Gc3r< zTv;;MVKeNk-2+!Q6@$_VbAWwpLb6qBp`1nzloZy1+m}1Rr}aHeSKDL3qt}^$6#Y!z zuU92Kg;k-+wR!7$eoYwg8%OlEJ|S;u4j}hJ^1(1)HZpC_2iMl8!tf2%VC+^3v|Hr> zXah^qxbzfh{NxGw^vMR0Qgex`5fKD0Ri=UN{aeAdE9=0jr<>q|thG&j1{TPZ5mOjz z6U4>E)&olHFfb^rMXuZbhkUr`4~N|&NQ+82VAbvbvw~fqK=ew`Iu;0Khtff5atdVT zWyp1Q(}=Oo8glsr3fvrH$Tc<<Ol6sDqw1s4!;tsVNcI-$p0AuT<>LIogbHs5G^L7%yr4m z(Obw`PfO4r)d1{9s$kj2a-h&ONZt`FhrFllOKq&CMe3YIz_R zH)g_45bU-E(JtjM#Zi)bWL6tk zRtW(~-CiWV*$+xBRRhnuRH4C{?MU=oFqE2-x7PAJf`rIEK@JWdh9cYTtlt=MJBejp7`zaMv94HOj*5AFH5&`cClcMIBJ0|B#o1b|brv03b#7k-rPGp=-)% zGCp(=d6|fUiI3@I&`xJA-?~2X$@NI+H|z|n|12jDUGxC5-_yW_d&bZ}ejhwsemF`VYZ8M8T&tmW^$CRI za1No9v&fs-Ze((N2l8#&86@=9lHMPC$eTlzH)p$&fnRSCS-3F)f~ToqrmGINB`<+% z)NYVFm!Ai|UvhzT#~HA-))%I>o(8PGClXt|3@kd*jl9=(h7-H%fqYLO@`6Pl&<=6@DoG_J{SGg=8e9dc=BXtq)A^*Cem5&W5vB65&y^ zop1!rBrS*YU{+cOB4t<#m#wb4Ih*eWLu<@PU+xwdvC_o4@tYR>RhSLz%x;sy5@X0> zi(y0UA*du|47eOXyk{B>R9&-X#z2zqN);NTk>$nea6rZ#y2X)*l>K$9`rHs$ z-`Q&Y_Ng2EcZq@8kn#uFJFpcpWDQcNCxE<{-vtuI>w(t)dDeWL1GGOKg@4oQz!I1Vrha`w6gr~7 zy2=6|JX#69+)IU-2`fqO-I1^@Bms~~(cCmmez-<+8*GbEgDb+*z^eXqa^OQblx^}u z;`s8xufmNmr^Fvf_7{Vl<^nL9{}ggssTQfck_&2my`Z+E80b28nhbRRMy5T!L2exC zL-G?h0R@kI#9jRc;xC*Cn3imCa>NW6Hazmykjr`qv*0oKFLtG3H$5j|zaNWCugWmP0k!6F{tg z4jJ5^2oE1t{J&Qi0_ksy;f2mgByeF5*z+z0xSmV_=_CL5tXWTz`gazfK4Iif4H@ua zt0B31>Jaj}Wj&NsLE-4hTF`glua$UBHVETZA@!TJp&>U7h`8y2t_=;a^3VekKGuRR z0;15duaSHx(np?os{=ZFl0ebRYd2dGj)7LWrQp>uJ{UA?1*>l=k{!yjKW zgEzmCxI22lIjV*$^EC!0acZHP!!&v1(Ne_KcojM0=1=YzIt@29r-4kKw6!HS1IUL; zfVCIG;KX4q_^sR?iM4WpscORTe84lZaYY6Q(93`mkqw}Iy#e%-kA>YaX<%Z77Sxl= zggUV*@PO4m@|VYHP;t5ld8`r+&fk)Un^Y~}z@MYAU~3I52wY4G**QYp;HzXJXFnL) zD+VQZg!&{h@j@5+HAsW>BSqw~aXB#Fm<`V=bCK4) zOW=Pl62^et%?$K_h?g@}oZ4H$htNG_z_vDV~+)&V;T;H)_Z#`Rmm)4Pfg!RSJ8Wn~CF zdeae6H`7EY`7JOpa1c4zqz)%6BEX#&HlVpd0_rc#CBO7ok{YS9aEqTP`8p^Kng3G< z$7Rj{*MDAQu6QW0S+NxyGnaswswwc1bSM&QodO=3Zi701&vTX4R)CqQQqm(~9T?B9 z1WWUUz^Dk9RD7%pZqx|?ktK(Lko+;^h1wkAt|I}%_!f|(+WVn6utVex!-3|hZRGMd znoxYS4vwDH0TPmz$Q8-w$c#HXKu%H_2@=$RMb$7egWs{S%0PB%PvhA;};HYwJjq_BpgWc zi$eLkB8bAX6mZ$1#X1aK2fcl5!C#GYz;Px6N}5K&GOrOb-J}}!HatKEc{)%bq!K2$ zHzCWur6N-bV+h@cf-UpkNW|?3lobeoeZR86*z=Rfa2W?k4VQw9o;tYJQP}$O-Be`U z+XafKV_-%3E7H|Z1U}ed4XiJgfT)~UIJ2({ar{{g9tKr{qZ)x$yB2jJLcfGSS+Xek z@vkm;_o{;QD&&XEg*w=93I%B@J0TjO3nHS=zzEK1z~}1)lwUQGORBTs0p*Y6$V)G- z|Nj5Gd{j0lt=x+QJJrLFKBk~({y1<~DFwP;{?9{2GJ0~20Bw835gYoFMDJLr$5gv> zc^T=ic-uB#;<=^1;rwWn*7re#wlZ=@ofoi^P{JW}@K*{c?Cu9Wt1 zw9bh04!KY~^?NTku^7V(yD-LcAAZMSi_Li@_p^zc3fnmnUXa(yNAXP9Vj}dG1c#FK z;TT^(KvZ_0<8`(_;oKJC5OSW*oIjwSxA|8H@75|8;+R`K`#yUyTbPTma?=}mJJ$H{ zU|$O{Uwe;enOM%rEIG&t5UXTg2Q;(qw&rpc)U`RA`LFZVX`JM|c-_cj0xGEc{?VL_ zJ^7pfi;KLE2lEL&gBnisXd|V(@)Y}w{|7$2se$u#;0=4SNrOk!2y>37-DGR4?y}qL zEqT$7GCb+EjXY23EFNcn1(Wx0Gw*}tIi|7XI#2q27EkbgCQo!YgX6TSoEI#q&TBgC z%6ZfLl-OOC#H!B`yz<5>o_T`|@1rf?z4=zban=9GmfzChxwVH8aa*_WVu#~7g?Y(5 zQ8tnNlAp%8*>ap$WF1Td$s+6`^E96SWxxqN=_e}P)=}{1?SMQK32*$hGSb|$+IiH$FxRPF$D^`obtOKoHU04BJ2>sYuDhh zS7T;bdB+&eot!ClCx0qo6t2L`x%RXFxka#w%i?)g%+uND-vxQ{Rc+KJtA~WhZ$qBr z`OQ3^l58S!%7s@DBf#r>FHJPH1##4Po0v}P9lX1JsvI4mMXaJ)3&YEG6e ziGSmL;qFh_0TDw8~UlPM5(;Z5(gU@d;+5vog-SXCDnUUZ*6=ahmW=fYwR6KbZB8mqB7uG_cSR4^;Z_Cy? zsB_M@O7iw>*oiyI=&;q3i+I1mQ+7|S9MimVlASHl!58rL_>Ma)>zJX%5lD(*z00NX z2a$EGbM6TiaI`q%<@0Pu!6bEgJeOU&MT#Xy#dt4eGl<*O$%Mqp5vE-60BiGEg%gnc zfR)QiBKCDg<1Z7A5XxvGv5sGkz_^$26F477=IbhA{pbK|V8Jq!_3dm~$4S~Z z{1F>z-@~-8ea^-N)G)diYYA;slxG*Sm^Za=9*4tBrxNB>7EZur`k2r$-WZ)E^h zN7a)x7Vu;>6?};N0D_=D)UldA{KQ73VpbyfEAe_a$sTN{+2!r?!8Eigu&;}+ z6M23wSjv6{d$1cNzB7~Thrq)`<60Fq?5Qgo`sfv}fXV zwt3Eu-KBJeNwVuEekAG<9Znb7yucz>#8`~IXcfsG&VfXf8bRoae#9SAv23b-H}>O$ zD(hTb!9Hh?5yxK%vArSr?4w>U*6I0GB2DTx5%F=FPW+X}TJJJv#pS=E-rH>1BAlNv zddE-fT(*svTH{KrZZ0Q2Y(7C4r-T#DOXqN@!J};B=4Hg+&LCpipd5G3RAXMMII*<> z!Nk8J4V3OOw|{ER((u~{oRYRP~!{Z@n9$W8icTX#|v4X;spA&W;z=be;5B0+(4bE$B1w1 zw6Va`$LVHqX}o>(B%$#whB%_<`YkEFt7$i*UZSZ}^`zMnt!u z525IPj>-7yMjzPc#f++o5k@{JVe!(JIq*`9SXZ@xX$t$`fjX+}smf7&iB&It^|uIX zB(F{vxXZB7iaUscze{jGdu2jN?FWrI2eNJ%!+6t~D!fL`kI3)Yg|`dL(&g@f#B{JK zt9I83ml)WLXO9Nsoroq2b>s0;Z&5vNf;q$h4#M8K#zLr<<4^ z%Zu@I<eNOQwt*ScHI^*#WpW*UjE5d7-V(W+W{df6z*`q8v zu)UMXb?agtw=Lj%U^AwDl#fX#7nrS8yK$H7e!QShoFQ8>@jA{ZUN<2_M8=FTocL$V zf9sd9!EIj5w{^?$w+W)`*^kwjUXdpzzGa~G1eJgry1mEl$gabHFF_lP_b_Z27jF(s zV$_K>^r}o)P>yp+py2)iL1CIc%%&Fs|I4gy%lF z#1MJ^XeZ}O%)8Kb`j}ZJzQJUWzSDM@alUGUXSr>|+qV3|ZpEcCGmH#lSdL*hWzKxF znxpk9a?n5a`MC7UD*SlRAU0o~gbnPxi}@y~VJQYIgGACy)|+jNl+G0@`nM~S=9h^- zHa^Q-GbZpH&ok7tfjf=<6vEo@OH77u9>Wd4fNQIn;-&o0v2k`4ZoVTM9~8@A?mNWN zM*D(s`*-O~rmH)1!9kZf@cAB<6LuEQG_s`2kD1^~Hu=nmc?xbRD9Pw}xiY#-^H2$V zDgG?`Iaa5=j%myML`(c$!&Fo!Kw)dPQ0fJEwb?g z*EZo6A5P$&xA?G$E=gLgcLUaNGoD#ro{TqT6yU2x5^=IBi}@9iM`_J&z%}t~+-*Pt zH#>L-`w(P@`-i@7IWr=I&ll&RzAEF`wpumXug#fyBHctiI%2|{at&jCT}iK zS2S?jsB74+<`DWcu=Ox{hV{cF@A%kFX<`)3D6X!I;2Sg1I7? zgz+zTVop2x(P0`+^gpr7)Zo}0^;o+d3w^wV`I^x{AB=iSdu;S((h`@Uw?bQ~oGW(t z_bEr_X@M5*;cA9+zfNO*A6@a0>(x|v?gIVP+m#t>4#c-8e!xN_meQ>y@)#>=%UCMN zFi!pgj5gm>dcQhFy*IPN{}NlUnG8YP^N%>*@2G=!%#F~!hqdvH8DVC9$0;lV6T^Fz zzS2L0U!X1D^jdqjA~YgYfj-S!hqun^Fu^)~sJXVg-3ya!O5i9T?koEmQ_R*ETyYlwJ@Ybhxxsbj@fNr z&a@i=`f^Gf&O@|fE$E-MH{_kv93r7 z+DKFY?;Kx)dK@UA{_A$7IhP{o;K%h?G5Qtl-5|tl-4#F!O=n{kJs;@7nHg+hE)l!3 z@jQLYSPB>V{1ycn`gG%cDSTK$8b4XE82?aGNG~`<$U@O}UvGZlU^i|6l zdTVqweOudzvRQtSPPjRQIrv?%)sdaT2GL%6TB-(1(u3G-GhyuWonCCmwk+)A-#6IR z`V35PUjbUQ?hSp@Lz3PkbC`O&VwX2kCG4Bzo?x7VXttg&HN? zp`BWI*cL?S6dvs{7o|GkJw3O8a0i={E!3u`bn)gtUkKm2xdVS^s_NF)=yD_kf-X4FMR!`BwGE76z zS&N(W`1c(8+vG8-dRHZtoN^!C|Jet7CqUuN`{ZLXGn>yk| zf6F(+h+6@4`Hdaei?Q?;Ls=fW_n{({8ETFF{IY|-?G;Q9eOgPK?ma~L+Q-lp=E}5_ zlQp)~#Fw7lxfh$13c+;m2w@*OMyMUv^jlgp-LNX9quAQ1_h?n@B&9AthE6(#U@F8S z`rH#6OmOTXdU$Lam0hHO38eebCl?Z_fxKioZ)0@pkN@X&lR4<+Ya1!WGs&p8?n%39 zVh`Q6tKQZvQV`>xPHdeswMLhZ1fvH=yV0)NMVOXS2)0ogrSGtgl-DacEIUXHy-*v9 z-4%^T6(V#|yLaQQ`-eK4Zzz~z#yOst`n5sSC@~BZnAnFoIbT4laBbQs$ea2+zLO?? ze?nW2FTq|}o}h*U>ZpZd?^~63bf9Nm&!ZhhCU%8!r_h*?F?4a7D%I%VM$KBTr0*># z)9UC$JLEz^%K@v!c678mnl@ZQ$%~&tv6;816jF^o%r*UgjV|<6aS)|!UDD_Tdb&%VUi)sd zT~K>J8WC4%dxM!q-&xL~ne#D}zSA`7v4TfAY0_xK?qq7t^dai4v& zrRc}35+<_rd+XhWB&wC%-rD?Slsb|Xf~q)FqVLYmp=}Eq)TnzeWfi@Wa&nKRUOtnd z^3M&SP4~W_l~q0zcdt0TeO4UPn|_UUk>%)A-yVuO+iX|*_8=vm^#D~HF-Eg?+@N0E zTtq$b7^0G1jO`;u#j~iUWG9+Y zy(5F@#!sKASlwM{=V%goj}JlV8-=Juz48?weVby|M_KWM+KJ1^r{I=gX`J*qkKd#o# zZgsm?9`vWpZ?Ns1O}F>pgMR(&MLD~pL0tiONuLl6A*tv@S= zQ|&zd3iJM_le*8JZdFgd*u{SDILjwTe>olV0(h|9IN9#^Ha@S8x4*A8Hh=bZPG`IO z_wSYGc4v+QoL;-*_k7R{U2C_0t*&1?D#78-OXpj<7~tzMu@7y|rrWvq&}q2nY~~d( zpdb6km##vS^IvqpLhc+!={PD3r`zWS=bwJ|rOyd+i?jUckDXih;B;F1kl!1y%^eqV z2aMivfVYE>oBJF=KVse6_@{&pU?E<|fmk|lV&BE!bB5dl8?}4>@^Ur1`BsHIqMle| z-_;0y{m$bzT#hdHSnpnW7j+GJb#uq(ZJfRywRHJ)acucoVvYg5t~M3Up#yn?mTU z9(9Nb`y3m1p~3Yda$dj=T*wu$p&RFfH>RL&k_T6;07nzcrL_@o&qkk!xp^oN6v_=3%P~PL$3uIkyl`$ zXP^`I*bnpl0&xN#xj-$@58>Y3DoIxt^76u;13%OXwY%tZ$@vQ3k~!?X7Wg_Xn0pMw zLvzH5{Iwybc}w3LpHsw-bMPaFumKkQ5gTmj{*~mzdf?aMO?B+}u7{rV`OCvyKJ<4p zjy=+qGX`T*4vO%eQBG!DcOM7cr!wS)^>Qxu3)*kkmh#tf zE(KQ3m=s4VbyD2%TuQ7bWpPc1OSm=EF(HG8_LFs7L<6o_%T&qx>JBKAq;Y5g3|p zc#pA7b@HC-t6Y5Bbvy${%$J9_+cH1ezfoEzdlrZMW=IZMCrgXOn$}>#Cp9<5yH+tZ z0!POzW1Voy{3;I1m$7RVo9gAs_hFI~Yi{sO@G{YNs(&*Tx9l+)XKJq&{MNp#$K<Ej@y=7~?Evy)(woRR<^eof4Y|p~7{9^bP zM$FIRvNUO*)QH)qo)pHad=^(skM>zbr*&re>iPfSjOnP2wY9$tXBoD%vp8bDF}qd$ zV`nK1)_F=(t=<}I+1$d6@mMzct>>1mn2s3lsv0a`OVfG$EPJhM*<<19Gr3OT$Ik1_ zj!E-cHr85m=rij)_QbfBJ+*qLG@S=q>(RKDt)-{V(y>{-8as8S_pP-xQW&vEYN@lh zHf^VER`pGtTW8vC6Wg25Tm-(eM#Qe0z6mH5t zWvex$_AI~HW8r9x%j_&XYaUz2<|%CJEQTldEPQKSuXw(UCpIiR&AG~cN{gjk^TxPh zxGCFJ?UrGlhj&#z{lqk-d@|!Kg>RkdeLw#f00960V_<-T1~eLliLs9IG99mW=*0svA%1B)Wn BP=5db literal 0 HcmV?d00001 diff --git a/tests/golden/make_dsp_edit_ita.m b/tests/golden/make_dsp_edit_ita.m new file mode 100644 index 0000000..340df19 --- /dev/null +++ b/tests/golden/make_dsp_edit_ita.m @@ -0,0 +1,22 @@ +% Golden data for tests/test_dsp_edit.py: ita_start_IR, ita_envelope and +% ita_quantize of the ITA-Toolbox on synthetic signals. Indices are 1-based. +% Run from tests/golden with the ITA-Toolbox on the path. +fs = 44100; n = 2048; +rng(5, 'twister'); +ir = zeros(n, 3); +delays = [37 300 900]; +for c = 1:3 + t = (0:n-delays(c)-1)'; + ir(delays(c)+1:end, c) = exp(-t/200) .* (rand(numel(t),1) - 0.5); + ir(delays(c)-4:delays(c), c) = [0.001 -0.002 0.003 -0.004 0.005]'; % weak lead-in +end +ir = double(single(ir + 1e-5 * (rand(n, 3) - 0.5))); +a = itaAudio(ir, fs, 'time'); +start = ita_start_IR(a); +t = (0:n-1)'/fs; +am = (1 + 0.5*cos(2*pi*50*t)) .* cos(2*pi*2000*t); +env = ita_envelope(itaAudio(am, fs, 'time')).timeData; +q = ita_quantize(itaAudio(ir(:,1), fs, 'time'), 'bits', 8).timeData; +ir = single(ir); env = single(env); q = single(q); % am is rebuilt by the test +save('dsp_edit_ita.mat', 'fs', 'ir', 'start', 'env', 'q'); +fprintf('start %s\n', mat2str(start)); diff --git a/tests/test_dsp_edit.py b/tests/test_dsp_edit.py new file mode 100644 index 0000000..4da43a5 --- /dev/null +++ b/tests/test_dsp_edit.py @@ -0,0 +1,135 @@ +""" +Tests for the editing utilities of pytta.functions that follow the ITA-Toolbox +(ita_envelope, ita_start_IR, ita_time_shift, ita_append, ita_normalize_dat, +ita_normalize_spk, ita_amplify_to, ita_zerophase, ita_randomize_phase, +ita_quantize). Each test pins one rule. + +The golden data come from golden/make_dsp_edit_ita.m (MATLAB R2026b with the +ITA-Toolbox); its indices are 1-based. + +SignalObj inspects the caller's frame (classes/_instanceinfo.py) and fails +when the caller's first local variable is an array, so each test starts +with a scalar local. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import pytta + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden') + + +def golden(): + return sio.loadmat(os.path.join(GOLDEN, 'dsp_edit_ita.mat'), squeeze_me=True) + + +class TestEdit(unittest.TestCase): + + def test_start_ir_matches_ita(self): + # ISO 3382-1, A.3: the last sample 20 dB below the peak before it, + # here after a weak lead-in; ITA indices are 1-based. + label = 'start' + d = golden() + ir = pytta.SignalObj(d['ir'].astype(float), 'time', int(d['fs'])) + start = pytta.start_ir(ir) + print(f'\n start {start}, ITA {d["start"] - 1}') + np.testing.assert_array_equal(start, d['start'] - 1) + + def test_envelope_of_an_am_tone(self): + # The envelope of (1 + m cos(2 pi fm t)) cos(2 pi fc t) is + # 1 + m cos(2 pi fm t) away from the ends, and the ITA envelope of + # the same samples agrees. + label = 'envelope' + d = golden() + fs = int(d['fs']) + t = np.arange(2048) / fs + am = (1 + 0.5 * np.cos(2 * np.pi * 50 * t)) * np.cos(2 * np.pi * 2000 * t) + env = pytta.envelope(pytta.SignalObj(am, 'time', fs)).timeSignal[:, 0] + mid = slice(200, -200) + e_rule = np.abs(env[mid] - (1 + 0.5 * np.cos(2 * np.pi * 50 * t[mid]))).max() + e_ita = np.abs(env - d['env']).max() + print(f'\n max |env - (1 + m cos)| {e_rule:.1e} (away from the ends), against ITA {e_ita:.1e}') + self.assertLess(e_rule, 1e-2) + self.assertLess(e_ita, 1e-6) + + def test_time_shift(self): + # An integer shift is a cyclic roll; half a sample twice is one + # sample; None brings the start of the response to sample 0. + label = 'shift' + fs = 48000 + x = np.zeros(256) + x[40:60] = np.hanning(20) + sig = pytta.SignalObj(x, 'time', fs) + rolled = pytta.time_shift(sig, 5, unit='samples').timeSignal[:, 0] + half = pytta.time_shift(pytta.time_shift(sig, 0.5, 'samples'), 0.5, 'samples').timeSignal[:, 0] + auto = pytta.time_shift(sig) + e_half = np.abs(half - np.roll(x, 1)).max() + print(f'\n roll {np.abs(rolled - np.roll(x, 5)).max():.1e}, two half samples {e_half:.1e}, ' + f'start after auto {pytta.start_ir(auto)[0]}') + np.testing.assert_allclose(rolled, np.roll(x, 5), atol=1e-7) + self.assertLess(e_half, 1e-6) + self.assertEqual(pytta.start_ir(auto)[0], 0) + + def test_levels(self): + # normalize brings the peak (time or spectrum) to 1 and returns it; + # amplify_to sets the mean RMS level re 20 uPa for channels in Pa. + label = 'levels' + fs = 48000 + t = np.arange(fs) / fs + x = np.column_stack([0.2 * np.sin(2 * np.pi * 1000 * t), 0.05 * np.sin(2 * np.pi * 500 * t)]) + sig = pytta.SignalObj(x, 'time', fs) + for num in sig.channels.mapping: + sig.channels[num].unit = 'Pa' + normed, gain = pytta.normalize(sig) + per, gains = pytta.normalize(sig, per_channel=True) + spec, _ = pytta.normalize(sig, domain='freq') + at70 = pytta.amplify_to(sig, 70).timeSignal + level = 20 * np.log10(np.mean(np.sqrt(np.mean(at70 ** 2, axis=0))) / 2e-5) + print(f'\n peaks {np.abs(normed.timeSignal).max():.6f}, {np.abs(per.timeSignal).max(axis=0)}, ' + f'gain {gain[0]:.4f}; level {level:.6f} dB') + self.assertAlmostEqual(np.abs(normed.timeSignal).max(), 1, places=6) + np.testing.assert_allclose(np.abs(per.timeSignal).max(axis=0), 1, atol=1e-6) + np.testing.assert_allclose(gains, [0.2, 0.05], rtol=1e-6) + self.assertAlmostEqual(np.abs(np.fft.rfft(spec.timeSignal, axis=0)).max(), 1, places=5) + self.assertAlmostEqual(level, 70, delta=1e-3) # SignalObj stores float32 + + def test_phase_edits_keep_the_magnitude(self): + # zero_phase and randomize_phase keep the magnitude spectrum; the + # zero-phase signal is even (x[n] = x[N - n]). + label = 'phase' + fs = 48000 + x = np.random.default_rng(3).standard_normal(1024) + sig = pytta.SignalObj(x, 'time', fs) + mag = np.abs(np.fft.rfft(x)) + z = pytta.zero_phase(sig).timeSignal[:, 0] + r = pytta.randomize_phase(sig, seed=1).timeSignal[:, 0] + e_z = np.abs(np.abs(np.fft.rfft(z)) - mag).max() + e_r = np.abs(np.abs(np.fft.rfft(r)) - mag).max() + print(f'\n magnitude kept: zero phase {e_z:.1e}, random phase {e_r:.1e}') + self.assertLess(max(e_z, e_r) / mag.max(), 1e-5) + np.testing.assert_allclose(z[1:], z[1:][::-1], atol=1e-6) + + def test_append_and_quantize(self): + # append concatenates in time; quantize rounds to 2^(b-1) + 1 steps + # across the peak, as ITA does with the same samples. + label = 'append' + d = golden() + fs = int(d['fs']) + a = pytta.SignalObj(np.ones(100), 'time', fs) + b = pytta.SignalObj(np.zeros(50), 'time', fs) + self.assertEqual(pytta.append(a, b).numSamples, 150) + q = pytta.quantize(pytta.SignalObj(d['ir'][:, 0].astype(float), 'time', fs), bits=8) + e = np.abs(q.timeSignal[:, 0] - d['q']).max() + print(f'\n quantize against ITA {e:.1e}') + self.assertLess(e, 1e-7) + with self.assertRaises(ValueError): + pytta.append(a, pytta.SignalObj(np.zeros(50), 'time', 48000)) + + +if __name__ == '__main__': + unittest.main()