From 71bf96370331acd58177e506cbd0834db90a44f0 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 03:47:47 +0000 Subject: [PATCH 1/2] Add the shared FIR substrate carried from SampleRateTap MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Five headers under include/tap/dsp/, extracted so RatioTap (the second consumer) and SampleRateTap share one implementation of the layer above the converter engines: - kaiser.h: Kaiser prototype design (bessel_i0, kaiser_beta, estimate_taps, design_prototype, design_prototype_compensated with zeros at k*fs and droop pre-compensation), plus solve_dense. - sample_traits.h: the sample-format core for float/Q15/Q31 FIR datapaths (coefficient Q formats, exact int64 accumulation, single saturating rounding) and the tap::dsp::sample_type concept. Format core only: engine-specific extensions (e.g. SampleRateTap's inter-phase blend) derive from these and refine the concept. - fir_kernels.h: dot_row (with the SMLALD dual-MAC Q15 path for DSP-extension Arm cores), dot_tile_frame_major/dot_rows_frame_major (register-blocked channel-parallel kernels, bit-exact vs planar), target gates renamed SRT_* -> TAP_DSP_*. - quantize.h: row-sum-preserving coefficient quantization (largest-remainder), refactored out of SampleRateTap's bank ctor. - analysis/{sine,multitone}_analysis.h: the measurement instruments (LS sine fit + tracked variant; pink tone comb, joint multitone fit, program-weighted SNR) promoted from tests/support. Contract batteries ported and extended (45 new tests): prototype specs pinned at SampleRateTap's shipping preset numbers plus a non-power-of-two L=147 design, Q-format/rounding/saturation contracts, planar vs channel-parallel bit-exactness across the 8/4/2/1 tiling, exact row sums over a designed prototype, and instrument floors on exact synthetic tones. README gains the substrate section (including the fixed-point roadmap: Q15/Q31 as first-class embedded profiles for M33/M55-class targets, traits over raw sample types rather than wrapper classes) and the second-consolidation provenance note; CLAUDE.md's discipline clause is extended to cover the fixed-point profiles. Also deletes a stray committed a.out. Behavioral change: none — every function is carried verbatim from SampleRateTap apart from namespaces, macro prefixes, and comment provenance; SampleRateTap adopts this via submodule in a follow-up PR whose icount gates prove the move costs nothing. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Ldeq57sBySx2nFTQsGcQB6 --- CLAUDE.md | 22 +- README.md | 90 ++++- a.out | Bin 106192 -> 0 bytes include/tap/dsp/analysis/multitone_analysis.h | 229 +++++++++++ include/tap/dsp/analysis/sine_analysis.h | 116 ++++++ include/tap/dsp/fir_kernels.h | 168 ++++++++ include/tap/dsp/kaiser.h | 375 ++++++++++++++++++ include/tap/dsp/quantize.h | 83 ++++ include/tap/dsp/sample_traits.h | 202 ++++++++++ tests/CMakeLists.txt | 5 + tests/test_analysis.cpp | 79 ++++ tests/test_fir_kernels.cpp | 130 ++++++ tests/test_kaiser.cpp | 196 +++++++++ tests/test_quantize.cpp | 101 +++++ tests/test_sample_traits.cpp | 111 ++++++ 15 files changed, 1899 insertions(+), 8 deletions(-) delete mode 100755 a.out create mode 100644 include/tap/dsp/analysis/multitone_analysis.h create mode 100644 include/tap/dsp/analysis/sine_analysis.h create mode 100644 include/tap/dsp/fir_kernels.h create mode 100644 include/tap/dsp/kaiser.h create mode 100644 include/tap/dsp/quantize.h create mode 100644 include/tap/dsp/sample_traits.h create mode 100644 tests/test_analysis.cpp create mode 100644 tests/test_fir_kernels.cpp create mode 100644 tests/test_kaiser.cpp create mode 100644 tests/test_quantize.cpp create mode 100644 tests/test_sample_traits.cpp diff --git a/CLAUDE.md b/CLAUDE.md index 20a7d86..344176f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -8,19 +8,27 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co portable C++20 (standard library only, no frameworks), consumed as a git submodule by the individual libraries (TapTools pins it as `submodules/dsptap`; the AmbiTap/MuTap lineage is where the FFT came from). Four primitives today: the real FFT (`fft.h`), the YIN pitch detector -(`yin.h`), and two pitch shifters (`psola.h`, `pvoc.h`). See `README.md` for each primitive's -contract summary. +(`yin.h`), and two pitch shifters (`psola.h`, `pvoc.h`) — plus the FIR substrate carried from +SampleRateTap for the two rate converters (SampleRateTap, RatioTap): Kaiser prototype design +(`kaiser.h`), the sample-format traits (`sample_traits.h`: float/Q15/Q31), the FIR dot kernels +(`fir_kernels.h`), row-sum-preserving quantization (`quantize.h`), and the measurement +instruments (`analysis/`). See `README.md` for each asset's contract summary. ## The design discipline (load-bearing — every primitive follows it) - **Fixed numeric contracts.** Each header documents its packing, conventions, normalization, and latency as *numbers*, and the tests pin them. Changing a documented contract point is a breaking change for every consumer. -- **Double is the golden model; float32 is the embedded profile.** `basic_*` templates - with `using x = basic_x` / `x32 = basic_x` aliases. The double path never changes - for speed; accelerated float backends (vDSP, CMSIS-Helium for the FFT) must re-present the - *exact* golden contract, so the double test battery stays a valid oracle. Cross-precision - agreement is pinned by tests. +- **Double is the golden model; float32 is the embedded profile; Q15/Q31 are format-limited + embedded profiles.** `basic_*` templates with `using x = basic_x` / + `x32 = basic_x` aliases. The double path never changes for speed; accelerated float + backends (vDSP, CMSIS-Helium for the FFT) must re-present the *exact* golden contract, so the + double test battery stays a valid oracle. Cross-precision agreement is pinned by tests. The + fixed-point profiles (`sample_traits.h`) carry their contracts as numbers the same way — Q + formats, the single rounding point, saturation — and exist for M33/M55-class targets + (Bluetooth-adjacent converters, eurorack/pedal deployments) where double or any float is + unaffordable. Per-primitive fixed-point adoption is opt-in and is a documented Q-format design + each time, via traits over raw sample types, never wrapper classes. - **Real-time safe by construction.** Geometry fixed at construction, every buffer allocated there; processing is `noexcept` and allocation-free. Numerically fragile recursions (e.g. the order-48 Levinson–Durbin inside `pvoc`) run in double even in the float profile — documented diff --git a/README.md b/README.md index 076dd82..a7c992f 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,9 @@ Shared DSP primitives for the **Tap** family of audio libraries. Header-only, plain portable C++ (C++20, standard library only), no Max/Min or framework dependency — consumed as a git submodule by the individual libraries. -Today it holds four primitives: +Today it holds four primitives, plus the [FIR substrate](#the-fir-substrate) — +the shared design-math / sample-format / kernel layer under SampleRateTap and +RatioTap: ## `tap::dsp::real_fft` — real FFT with a fixed numeric contract @@ -125,6 +127,85 @@ per analysis frame, with every relocated bin rescaled by At ratio 1 the correction is exactly unity, so the identity contract holds either way. +## The FIR substrate + +Five headers carried from **SampleRateTap** (where they design and run the +ASRC's polyphase datapath) and promoted here so **RatioTap**'s fixed-ratio +44.1↔48 converter — and any future FIR consumer — shares one implementation. +The performance-sensitive pieces are regression-gated in SampleRateTap's +instruction-count CI (Cortex-M33/M55, Hexagon, ±3%); treat measured claims in +the header comments as contracts. + +### `tap/dsp/kaiser.h` — FIR prototype design + +Kaiser-windowed sinc prototype design for L-phase polyphase banks: `bessel_i0`, +`kaiser_beta` (Kaiser's empirical fit), `estimate_taps` (the harris length +estimate), `design_prototype`, and `design_prototype_compensated` — the +zeros-at-k·fs variant with passband droop pre-compensated (closed-form, no +FFT), which turns branch-DC uniformity into exact transmission zeros at every +multiple of the sample rate. Runtime design in double, deliberately not +constexpr (the header's design note does the arithmetic); run it in a +constructor, off the audio path. Also exports `solve_dense`, the small dense +solver the compensated design and the analysis instruments share. + +### `tap/dsp/sample_traits.h` — sample formats: float, Q15, Q31 + +The family's sample-format substrate: how each sample type stores +coefficients, accumulates dot products, and rounds/saturates back to samples. + +| Type | Coefficients | Accumulation | Output | +|---|---|---|---| +| `float` | float | double | plain cast | +| `std::int16_t` | Q1.14 | int64, exact | single Q29→Q15 round-half-up, saturating | +| `std::int32_t` | Q1.30 | int64, products pre-shifted to Q45 | Q45→Q31, saturating | + +**Fixed point is a first-class embedded direction, not a legacy path.** The +Q15/Q31 profiles exist for targets where double (sometimes any float) is +unaffordable — SampleRateTap measured its float datapath at ~19× the +instruction count of Q15 on a Cortex-M33 (soft-double accumulation). Expected +deployments include Bluetooth-adjacent conversion (RatioTap) and M33/M55-class +eurorack and pedal targets running TapTools primitives. Per-primitive adoption +is opt-in, and each adoption is its own documented Q-format design: the ladder +of headroom bits, pre-shifts, and the single rounding point is a per-datapath +decision. That is also why these are *traits over raw sample types* rather +than `q15`/`q31` wrapper classes — the arithmetic contract stays visible at +the use site and pinnable by tests, buffers arrive from codecs and C ABIs as +plain `int16_t`/`int32_t`, and the SMLALD kernel's paired loads stay legal. + +This header is the format core only. Engine-specific extensions (e.g. +SampleRateTap's inter-phase coefficient blending) derive from these +specializations and refine the `tap::dsp::sample_type` concept. + +### `tap/dsp/fir_kernels.h` — dot-product kernels + +The FIR hot loops, target-gated the way SampleRateTap's optimization campaign +measured them: `dot_row` (planar; routes Q15 through a dual-MAC SMLALD loop on +DSP-extension Arm cores without Helium — bit-exact by construction), and the +channel-parallel pair `dot_tile_frame_major` / `dot_rows_frame_major` +(register-blocked 8/4/2/1 tiles over frame-major storage, coefficient +broadcast across channel lanes — bit-exact against the planar path for every +sample type, float included, because lanes are channels, not taps). The +`TAP_DSP_CHANNEL_PARALLEL` / `TAP_DSP_CP_MIN_CHANNELS` gates encode which +targets prefer which layout. + +### `tap/dsp/quantize.h` — row-sum-preserving quantization + +`quantize_row_preserving_sum`: quantizes one polyphase branch to a fixed-point +coefficient format while preserving the row's DC sum *exactly* +(largest-remainder distribution of the rounding residual — "the coefficients +of every phase must add to one", R. Bristow-Johnson, music-dsp). Plain +conversion for float. Design-time code. + +### `tap/dsp/analysis/` — measurement instruments + +The quality-measurement harness the converter suites share: `sine_analysis.h` +(least-squares single-tone fit, frequency-tracked variant, `snr_db`) and +`multitone_analysis.h` (pink log-spaced `tone_comb`, joint least-squares +multitone fit, `program_weighted_snr_db` — the program-weighted metric with +Fisher-weighted ratio pooling). Instrument floors on exact synthetic signals +are pinned by `tests/test_analysis.cpp`, so a consumer's quality gate never +silently rests on a degraded instrument. + ## Notebooks `notebooks/pitchshift.ipynb` measures the three pitch primitives — driving the @@ -175,6 +256,13 @@ would silently miss the other. DspTap is the consolidation: one wrapper, one contract, one home for the next backend. The unified wrapper is MuTap's backend-capable `basic_real_fft`, generalized to the `tap::dsp` namespace. +The FIR substrate is the second consolidation wave, moved here from +**SampleRateTap** (its `srt/detail/kaiser.h`, `srt/sample_traits.h` format +core, the dot kernels from `srt/polyphase_filter.h`, the row-sum quantization +from its bank constructor, and the `tests/support/` measurement harness) at +the moment **RatioTap** became the second consumer — the same +extract-on-second-consumer rule that created this repo. + See [`third_party/ooura/readme.txt`](third_party/ooura/readme.txt) and [`third_party/cmsis-dsp/VENDOR.md`](third_party/cmsis-dsp/VENDOR.md) for the vendored-code provenance and licenses. diff --git a/a.out b/a.out deleted file mode 100755 index 67cd4ba0df0fa52a677d7a8caa21e480caf8a446..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 106192 zcmeFaeSFl_weO#SfU%Y#bc9;1^xXE8dwXv6_TlI`wx^z2Yr zq7g;M7mX|`D>}WXKlR1@<>7jNukS-8|MK4@E?d->Z?zkw#1Xy+yW*n!cfNOhtKHH5 zwSUUhaVlC*Jsf@7np9!@o@kD_ zp?_CzYMAsN>2FNf?(urZkKW%)!*S)m^Zk9DcH+w?{gQbbCmKe zM=8JmDCOeXk@}r$yH0!BPXAdp9w6->WNn*mVv(GC?o?TfqDw(*1I$(yo^n`QIKKJY| zcc*(w&nrj*%)dUYe)qp0y1g%d{rT)ybWBkR>0{9WpCTXOf3ba*BNTc08uNv}JFbZO zkrEXCMS%3e`jO;ce(c+*u49XWuJW7j=-Q0AYVQvEH7E$2+5b{^Xuoka_P#vjud@u1?8d3pyHxzL3z|l+(&ysx%7Dc zUuQx2iTOgl{g58@uJPsn)rRoLP~a$YH4_WU3-$e!g7Um=QH$vX<+2<3e=`fp&EE+t z<`k4WyWuMCD<~ID`F{%v%JV!`-NJ%$CanA4qJnY+sr%of1?2;}OYtuY$`O?Ae@_;a zpU_=eRMcEh?iG}`7L*t2&(#Ix0}J5)EBRj>_%9Cp7YF`}1C9f`z0>vujqjHPP5s`P z9VrUtH8=F>S{*blEm>j!U0?nj#a(B_`1DR6StR%=s>yBa>gt-j!16L=xtDtKGGMvI zJ$V_f+`snZWw3JpvnMY@m3yEkF9VgkrzbDNl)J4bFN2i(W=~#*D0fv)UIsmPNl#t| zDR*v9UWO?5xt_cXQ0~(`c^RHupPsx7POjrne!N|0{Nr7bqRHOrRYk^+t}{;PmG9px ze~8_DA^e_R`5nFTn|tNo>Xl#DEB`{T{PJG;KlRH0Td({-_sajGSNq-sEj$l55rnxc2dL{)!iq9yApSPC;>9?~#g~{2Oe0ZnN zrO@#+KK@U5lST1s;^X3D$6X!V{eHQ%Y;rAwv_mzlImpb7&7o0{UJ<0L=;PT^`K z(r|<|W!Vq#6NAjG*vMhc|1R;W9`;_4hTx#-*Rf`xcmGPY6Rj>;|Hlw5W83|!+#}TI z=}T6AGmN(<)K+$u&LG14-hZQcAeYPYaZN? zgu@OhKbmzmood(>os7J1_sFg9Oz~=Pl?86-^c%CGPtnk)2kGD1aBS$+5CSOmjY~@6 zjY}f{*LQ08Q{s0YnDZ19H1~JFXzU$v@64qgUR}Qc8auq}?K|6epsOoLXT_K-d=|6X z_5EX-bcrj9M09x9c8W8a7Z^MXmR0NuT0ZI%?ENUn)WDrqiZA$4Q2*PQ?q^G2t9?#L!*XX&PZbVOs(eo?_yXiSZriY#bQbmz)p+b4VFEt3Q9*4>kn-mOd2YpiiQ@^F2>>A;Y0i1{5=UX%f3 zckI(t?0kv5Paa}b>7BWj)Ts%;gUpw%x@+zAlT6-a^Ha5Ln0Z%E{HX5GeTEbm(iqXf!vbpWg(sP z3Mm)FAn|1*RwSRg2ylqG&=E?jLXuK@Me>bS2YFLJ5_!_ew_!Dv8j&3lrpmHAu(VEL z{gVk@FP-eoTZG^TcQ2K`9oBr^H{z-IH{#dFCoCzDXE+$7UkuW(C)35T>~r*&Y$}cg z=@udOlj$9r>H`zDg;>|oHuxC`*iJseL-h=;Bmo~jU!Z+L#$ zjVIB#JaXH$-GY2N5}fo#wI^CO_erMrW)F61S+*o1%X5oS0FxIO{8zi-y&m|!f`7gT z{&VVMLE=v|)4$o<9+Ur*_!AH=|6=e5m+bO_%GH54YH>0%d>L|5wC%UuRj2n<%~Mrw zGHpI?42J}n3ol5ZQV;AsQuP3;PgV5>s%06oXE9C|Qz>1R zzs~sht3|BKUuT?eGS1gXKdUi54^(P#=Vb4kZFj)RbJpJWREL)mTXr){#xPz=qT@*8 z^U7asT*3-pytW@$_u1gQYxd>tw6V?6{QsHpnJb=NlKbC#G$Xj>q@N-^`Rx2);EyVX zZBL}vTpy%22958vj~QD%tQiB5xv&M;t`m=Ih!lO$(}eP2lmwZ6XX>l6v*AQs{T(w% zWCj0LnAj6ZOb$P#s$Va|G+oHnuh5Qn}tAjM& z(lg~kG_R|{3!0u7p)~5{OCM?}JmnWum4#KGc2%7)J>N~Ol>P=VG zT2QrKX{uUW)#C+K&0*CqT-CyYsz<}BCRa75plU%_HQlPdJ4NKdof$LL%XO}BOulfU z3NLeoL-U0rRXEHQ4pJd%aYp%dl;%#M6nw&0Rkan2 zDbK&3xT+)Z?`~IhB>qiyRY&6A6|QPdL0@_P4YjIX{5!=J9?rl1uJCaF?ds6jzuSv{ zZ&2Ece?ikRS}e%@o!FxlLB;ai?+gn5F6min-TXcR{Jg&se%zmfkB`q>))#Hq+`U6e z?`IERK5tL@D_&xG`cB`AjvCa`*@v(4^r+%=drS6o(jy}6YdX<~=}{#?X4D|tG6ju? z`g-3R$Bs8he;8z{yMlRHZ(gXEnee7ybKzg)f_Z7_x0xe4>=up-sxPhCL6TmF7oDmo_5Jh%q(Iuy(N3i zuCA_Od*JPrnX}50>5bWb6fieK5oS5cJEm5}VyTx(yCS}gj zhgx==9yzq;P`X`%h#dNHeDA@a{-xb<@7m%A;WwPEZH+f(Be!pW$NjF4+l~fj^4ezn z|4V~pB|evT&mBxtvZ;Q$x$S7SIdlXYgNnsL`cNXhHqmmhPfg`Pukm;Q6O~J)Zt%~C z?aBURC;nUdBPL=RUGy?g#-{XDfu{WSeb^OsETv+mZpZPgC1@tcCF z?6U6uZAVvEB5mE#?FZdG9?1TLnq>Ot!o(nji4ivItDw8zTX#Qyl49d(eprdh?Ox+o zi9RJN4_H4}3R}f?<7UvZsSj*SMBdMS#dV_gbkb^HjqCMQ|7~`dJ>%M89X+wyAe!3v zN~R0gX#G84_qaaDE~SQX6<}gq?ZR8p+Rer-uDh&t*Vq|$Ck?kM zoH*-oA|p&z?rQjymgtMwn~&P*XxFL3_muF_sWZA7V_%!Q*IlPtQPgThX^dAK1{T1qN}YnOyX^>D`Al~Ysp^d&6@#{LDLV-GihpcPRZ(S_BPXR z&~!Dc2OCwnui;Y+5!tT*5smIM_PBM`Q#rs@x_NV=(o?x}o3#x2+a?4(l`pBXNBt4P zEic%>Kzo>K{+1De*~H^OC6}R^ z2k~7l)r>G7%Kq2^pt~;}$)ZcFdB0Fabl662#)!iNyQXr}ZR0M_l>RSaMKbx#{BJ^7! zs(3%&M@=Sn5_lVaZ5k;1BWex%Z>%MDJ()&gHO1?Zk3`EBG*nmiF~E!$Xq!ZODaL9< z_BzU#;4!dcA}Gedidz0Ac3Z0dCfrT+09-Kzmds4npC{oY9#SGnNv1!n#`3af_N(XZ z@$O$pe+=u5>S0|9WW7a%U~sI*$+Yn(3CY zxRjvY9}|r$A~ltVZvR8>j{G=fujX!%y^=wIJx^;2;7aDKle4pPj7>VY3_<*8fHp=lV@G*@F~@~W!1pr0+Eb#ff`^CQ>KxJlNJ;Q8#^gK9m*o^EPM?yxxeB|-^OoQ>p3mm{!M!p!*Ue0#Vr}-9z$BZV za5nFUszengD>@VgpZ^deszHs-+e&#&WE0&95GN>^HbfN2Oy=GrhGp@z^g%?Y;nO{0 zlZdR%_NNbQf;abhVj8rRDCQL3SYUGYJMk@3x}7gn6+Q{tDwjGrNse?(ER^Z`ccO?#Z;oH_`MkQaArkNpeA9j!=+pa-}Qs zE=G-eZ!@@ z&LjPzOSh?fiA%qt^zU7|Rq0>3^dFV}nM*&e^bcM7pOn7Wr57swU6=m8($idej?(p{ zYtnmj7qBx|P$ft&+cqFfOY7{TH57Sf-FZh8hRYAM}C+s8I zA==B}1yZ>s`-OM99Sj#QSlcEp{F!(`F)a$?1;PKf0fS>Tpq_xiKWoT;Wx$|2Sb#g| z7Za*42o}5}|8E##OQahAuP0D&BlGhY0tJ07P;dZ>4hs~_cw0)fFi`NkW+My~jHC9G z1q%9Fpx}0+<70sW`X*3NCYKF0szAYkw^(*Pfr9+LZxh9slKVt~IjCAP=;BB9WgBX@ z_4|$WbF@H#+W~%EFkzq|1Kz^}1+K>f*%_|KMb^*90|l zb8;+D;JRDqx|8+E+ePi!djbWn$Jr&~WOtxol67&VpQKJN#->^zIX}zNaD_v#C+&0ZwPn6)-syu3xVDf?^#RsMp#s|J5 zP87rk{yPZGht7v?FSqXhOnl()9Kgro15Pva(MG@N@c4k#g2e|mF_^>R0~=ls{iB$g zDL!xtiN6#ds5kF+jv))<0}l$3#Rsg9Uh#oDs5R_Qcl*T_AHdK5lZy|qoVv2}0dxC_ zp7_AmD8n)xJwEU);XZ170KA9C2k!bS;sakCl8+Bm(}^Mm{j9`4{@|IaNwszm7Z`9L z9~W4ho1Y(_#RcA$SqkF<2gR*0E^yOpjOp;W!14RUA>snXB;db`pTHl5VTlVEf2Pqc zlob~jI8BKDTwI{$u%2>{eAGKGV6z2I#RXQc2U|~EU^Bq^xWMbAoqwNJTmZUs6@s_` zbZyDLN<*V79~bzZF#Q#AfnNx|en@UcH+{BGc}9E9?tO~+=4SzEae?0hicVStS98BI z>gJF%$H(ZqM7$68k6+2t2MUV|Y<|_~b4KqI#0C0?^j*?JpY4a?Lx*S#`*{(xPJY9F z?sEOeH)WsN%oik^b^3T*;0kanF7OYu$z5b4Vn5;d;Nk+?X>TIL9wi2p#RcvKGPFOr z?^>%l0(0>Z=)0$YzH@rh_trY&L#rBmlDNPj!E1jT?WQl!hZ73=SpnE#e7M8))8_h- zb%bP7ut$yy{3zegX&>w79}4#C6MFY^g6pT#^`j5~{g~N0a$Ml-d_QM=te-vNc{s0o z*Y+CEFDO91KXF`Ox!}dO@_auZj|)5i%wLWRtY=JrE-ui~>HPFWPh21r8^r~7ydtR~ zF7VbKNsWsOyhNH{;ar$z^44&B&zj5cqkm@&i6dwBnz(81+{WC+1^O97@|#%21zx7b z;c)>D&C(VH$nB!XydRX05A+TUz={HY19#S&Yv$TT2Eup%2M^xdE=bD919p?{jt8vj z8HeeP_Lsvo^A$G;;m-}R7FuTw|I!1WkBK~OO^At1vnDPkQg2Nhzk1>U7ooVj;{g*t z$m`2q@qp7l)GWJrz-@;tZGnt)xzE8893hdd*pvGMRoUruo{O?yW}9F7`JO1iNimB8 zL_Wv1`T}66zERa-6ySJQ{YzD!YSkAW`Ww>kyd=;IJCt7Q(jU<*yV9j!SNYQ}{esf} z;?j$i{$DQrZ%Y5zrGKe(#-;yJ>APL}0i|zu>F+6hvrB(R>2Hu$6yT3mzJi@KIziZ?%7|v^c6Z*TW)_VXsQ|1^!Hp8T=ge9hLGx^**!d-s<43~1x+cfX$PN;dsYRA%bB;(rvLc=|={@4>*8%bhDZ@VGX$xtLp9 zT!b0l?;#MWXuxKMnTnP1nIvR#{c5poo@RhX@A0=?hbtt?xpa1I^W~Z1MObtNf`!;?Ey3JVp5Te(|8{^k&`3<#Joqpw2tavhVgg zV)k8f+hF_lZz~laNw;|SE~BkGk70In1v?WzSXJ3}S4B94ZIAvrLR{SCp(p7-*_4rp zCo0!83`$n4N~T{Nwuh52KU?}@R~Hv*{Y1JYoLRfeP?O##N$Q`dSe$6t*N1d8qJqnc zm5dix6Lh}+E=*lKb5$2dX2O+hn#9ps(>%vee-=WWo*J71Z!bKpEP_SlCUT%Vs_E7G z8icFLl1`6}4Ns((2Y}WxpFF{^Z5n>;rA<}8Ev|}W)5}`6pB`DZ=Ac$td7_@nFx)6E zva8fP$F{~lgSCz zKQnjqJ*w%f-&UpfK1EJJ^oCM(et^z(XJ<*vw$md^*6at@@X8h5ym|yccr#ND<))Ah z<9!g6jE!ODn#I4`pni`8xePME zYnxf_UhNu&!2Q-$-oQ^~#zY#I^r@|LK|Y3oxHS)A-!J!U~NokfrN#lPgUuK3a5#|bXf@(==zTN0V#AMEyO z!5d|T9lZTzZgQ)if@=UM3YqAJfz|1+uJPkxZj znCfyGD@e~1VRkH@dw9Qi=q`wJaFf01S=lr$m>W#8qCH@tYF{ewLD_OuGZ%xl{UZ2x z@t=R(yKd!5c1tpI!903o)q`Oyz}r>{m=p+6@7yedmfRSc=;;_%+Ki_X&%h1+o4u`- zsWQ4uJ}Z`QyAKcKho`U@t_BbK!Nu4AFYIgIKW|_A9>u;|(Dx-i>xwZF(c%|nUlWzB zN4Bpg-Ew66+Mkge!M^gWA9FOXw|z}!ru*vm`VxJfz*TymwP*JA%Z8}HvKFtDQ@-;( z*8EvK)MF-=R?f7-_coTG58-rwbZdIkXJPc~tsl3hS5g7y{}mYwn2l)SQWth;@%W z#l)2K7)9C9hd5f!mbgI`+R?wK{7826_#RL+5ZLi2Wk)B;jyf9}L3U+{u%WP`mksS> z_VcEA<>I@(iOqD&wfuce`+!%mowALjOL;z{uEp8R6K{qVvd=lYxfQ#)hAht`rEq#` zS55Z(Ph&TEG{XPXWUqGhWb+;{o6#P#dHa2bo6KpS%w$fRg5C*DCWAcw&XanZOq5k< zGW)_#c`ar#uRMau>>HX)=eJ@viFNwLTluUj&UV{PjT^hk2x}SiUu`!p{Qs5Re5{w< zJZI^jwVUlr{tCPKgaQ9w+0E-p5&l!a({r88_NP2gvW%)AeOXCkR~dU2mTz1U^xg4) zM->$%>Q`1rR;FKQ*<2j2TX<;d_u_T9LTtM8>5Vtqw%@8PA~-YyyGJ-jSF_kC|}i`t}@1r=*m?0xTF$X8`X zb;TRGN&mg9Y8$r}Pxfwot0q(2m#f~H*ua+UeI|R;pXFh>WMm`ZHf>9$6$5lzJGbwP zsXWwzVlAd=vmaD0yZ!Bvdsg&q7}mJu0AyBzY+Vr9Poukv8$q3a192FgJ+*fwj=w95IzmK)GjmLPx_*Vo*lr|+Y7d&Ngq^3SNkLlw77HPjEU3^c_ zxVSi|eD(HR9^LGJk#o2tpY`T10NavbUAgPwnx4xs|AqX_4_S?tF5N2V`Yaal1_ zP}JDPUR*E6N*|Py8>C;3rw^<-sFi#0)8<&9Fl>=(Cov$FeWt~<`Kb?6iz-@8>7^xn zoYHu(*qd7cPH%oapN$9m-hMr3&bS7^VL1Jz?`=E9t>f9(GtH<`I`%wds`z%8>d$()jD6jC{KRC%7*Pp zWU9)GlbPZr!Y9_8e&6HV$IptpONRq=9H}S zko;4>5asO+(M$c@BTR324zi;vy~SZk(so4lHWmmi&f;m_sj;CjL)v~>_PA!yCA_kbm|FNA7aiIk`uh<=0C--gb3R|F#~?Nnc;wxIg0E|65oRPrs?fLU+fn_oVdJ z#@8IsuwiibTH$VR_5ntasHZeqZG5~OXz_0u&r=_WG<=(_)rO+#^n2#sXbJ~u;$Z9? zdSQe+n;o+fjXN3rnJk7nBqsaR)5e&7`%f?{&*n!+H-OsjwCQFWn|O(R5LcSx=|kb08yerAj`l3^?$^zLC#5q~A8A(M><`#^wN_=$ zvi+Giw**zsQ||ZgvNc6jkY2(X_?|Xs3z~{DWBOD@+JpL4RgtZU%R!br? zKt4(TU?GYRUExWhC25p=LP8R-Q3Vpi&TymMe-urL{EL&q*&5N}c^kvJR2L}N>F z(72*ZRIP_L3=$8Q`p0?EW4&l!%E%P&YHaQkRtzMdx2vb7qpBj`12o({T50hY*E$WH)iwQoaaRwI#4!?(SE9@8wjb|$IbC|TZ#MRsxrlD>x+ zR#j|Ea2EwxVQx62sfw*eUP_B3IGIaeUxyPGCpN3{0>c1@dC~JBi+TSb5{}^GK$Twf zB7g{l)USrERTVj4UJB}YGMc@VlR3l-M16ozU3yhx3ve9*FbOre{Ho#UlPMh$)Gw}h zUE+huvijkR88MF>`^SX%&c+Zd^rGiD*7c{RK#t>$T?q0@*p{f@kf>i*#Tyl#3({}t zbr^oI_l+PkV6Ak3VlmBHbV25V^MmyIAibM#$X8)hi{5=g1cZP;d-1PfG`|w0^$`DR zU7iA6ZP2p0Z&0y5sCYrKjR;$bf(BescLFKz`CwGLZp zMX06pN8)PR0{U(z@KQ$DkH4j50pBQV()w#DQ`6Oojq#>~6=i*gKgr|Uw8@IG60qWv zsv^yZ9uKGcC!>H)QthKCdR613T-H|68djf<{=c5(Z&gTgPs1;IQQid85L@IbO{Ez| z3FDc|K+-+fF9I?t*6#YHNp;Uaq4<_OR?ok*y{8f6&d4{Rzf0Rpty-5C zgYUaNt*TD)dq(QvIh+)#+&p9BZQX-CSwnm>XhZ_!4}QHZU_u>X;1KzvY+(3-z&u)} z0ZI`lw}{!UiY+Xo-30Yco$4)dOolV1APkkBX-0g~9vCcT4^)LSv(z_kWqp`cE~pOGE7ht?w;(7L zyk|j5ab*bB(@%*sumPg23y_3;K*@lLReUn7FNu|<_{rd++hSbzyp0zdIJab|%?usu zau0jr+)uf-&@d2n*lIq=)ACJo4|5Do(FeFDMR6-4*U5FIqC8;k0ooOKI|cq7j}Xeo z(T1mMT)gE#A04c=(-vjZiGzmEGlln2n$M$sszBuT9B_WmfgZo-HxMJg=ZAcop5TC@ z-==<%B4oKUnHh`&v?Sg%5JMK(7dkyH@AQ>*Z8)K`CR23+kCC+z4e?LFd2PjQI3cL- zP*A4?Pf6csEGct&N*A0jn4o2&;7`W`UL%0y`3mSp2OEM!CVBoKaQC+|Jb;q{JvKyr zEVcdR;CC0OU+V-*IN&65pAl#`=93&|DyvQ!l&8;&Mghc6H*!nNEqh4UV(O>HI-O(o z2h9_XSx}J;(pdMMxM#RyKTM`4#TMeCJw$cUFkDW8Hv^VC?`(NIQ~a4K?8JlSp4DWE zmruq&E9Uhb{%3H^aL}5|NGD1n*PBE=`WTYfF0hQ;6=TlUeR{{*n$PQMGF!Uy+|;&AOT>7ge&qQH%0jsc@Zn}f|EqmvHUT{o(1@^ z6{!NX_DlN%cX}SS8}*TvBpVbh^X?^CRXH|xz?=Wy6y3agsnpux&iMhH^UYP^0~w)n zK8Y1#VCbQM+7k=^LS-d+RQ;Pd%F(V8Jt8>aSzL}nJU z;Bct%AV~DEgj7h^NY!}$7aSIp4H!^QnuzFWM(*p<+ak&=C%iZxZDAU_^nBq1_e%98 z4B|imX<`P%sJ4Pls1Xm5>KwUKR|x6GzUoZ_HZo~*ni?3{Dvp`>QG7BOi1U|j#P3tv z!5YCy;+IG5v2bvV`fwanm7A86f`zKqYM^6WJjkGrZZi8bK~??T1mWVTR$FCRk#0$J zd!0L3f!w*MsedW{j4VS{dOLf8c@hyt33P3P(-Cg-V=qq$_P&+Gv1SXSeg~Wqr1vHf zttkQTjl54S*p{?xDZVl@V4!vt(>Vy{1$Uq2%}haLJ&zjO-ic%PWcmVjL{)li-RdCn zDl*xT)OO~Nu{P`X?SQw{Rjf}`tWKuCA6so@b7Kqm?%62)Hnu)2*J0|D7`Y;F)Kt72 z(4;x6TBOpn_DV@eY_(DbH8(cLD%?J5zwOL^5_{aL=f)njI^HO0TQ=?5=EiixCTL^GE+kUT22oUgL+(0Vum?5seb5;v$$i%*oEUs)f&sR;U2uMCbBK4hkSL$ z#^@);veC3ilC~+tl4b@k$Eh617LRUVlca}RMJbpzHTGatdQ&p8UfEipE>%O#F=>;+ zHZ@tOt+%+w6p zTEYjdL1eXdqwNy^^`71%WtXaic_DZLX@|>8F6}$`k4k+uz2268bF#iIlBKiFqcf8> zI~^eZ!2crErr$m%?cv;y)0O1ak@aSTKO?XN2hs0|*e}rj3r+=4%x$9~iV`qZMe-bT6 zMqe-Xd*)q3Bs6x4veP5$9f6{{92|*3CQ=~>Qv-dWIgGf%>12aH&{xw(q4n!Mqh92Y zx`{nb0E_b6JiQxPp1_{be%?@pgEhVmH6q2xR);drNGPKKP7eCBve$zZ5o4`G32Vnc^xEgS>!tQ|U#ccaVm zppdCF1gF|Vq`vY~4VCL*spx(@j|>KMBV+nbQD2ACXH1p4aPn!xezE{<_gdC1VE+^Z zQ$*2^6HM&_LiHfUR)C)=|GdzOo@OB&S7KWCMvf98g@qPU!1rH4lx=ZY)8Lp&7g@+n+e+mg<#Ppp%lT`ulERlhuqw^#|y4_O&AV$f}H~_C0G#Hz^`pAj&VX`MB@N!R^|A^_L&N&n1?K6xWK}Y z_trt8_=-fR#|@T`8oI$>1i>sJ8+<~=J_tO4JX@caIYcxMX@)d@T1MMrC5>2YX)|$F<|WRl3#k zPVv=OW-ic47AHt-6+BKs*6~~_K%8z5RlxcjNMR@If0o-|n>by6p7hz=4trjeXHxJA znGJj|qk?>j@8!)b%+zWsIEycTm+yDTKTA5xve?SsYMdh0*r)g-M*nF_rudq_ktHN& z@)11?4ed`EG`ZRf4uQfv3LKeJaV0=Kah-kq+1r$@S4ORb7v&TT&w8 z*@@Z7Y0zWk3HoE&0e?5qm^CS%tlVxo(NTGAJM1hMl98QW1v?7LuA)up8{#N;E^6tS zizc0vU#9Z+yoyxpNlof?WHou^s=GvWdqA8cb7eb>!PjI$th9(gH+6Mavg! z`Z*3L-7X@ms6C?Ov#NbIA{MR$$tn`5HI%t_<*$jYyuR6ABM3-FW?M|iQBoJlQQ>NtG^ zG|5c}FZvx(Sw0I=s=e9+q*Ie(_o=*{74Y3Q+LZ5C!%o#NspcmT+J##$rG`E#4yl5o zD3e83n^n7oIpX$GEowCNGF1`)|5Y4hhR*=q4Nh*__}rObMx<&W_41)~@9TJJwogON zfqzw*J})T)wI8UycCi&JX_%y`PbMRliBCX(UFw0^^UBw=;@8{!inYsE8CTyCmUjfQ zQUpIIf>(-F&rwnL0u|+)qNnug@LB}6?s@e+rVMm?QNDDcHqo&}fN$`T;@VEF!0IO1j1j&7 z5s;q$t^l=Sc$ql-XQ03=Bnb<=I@Sv_b*&a(`CEL{zRiFRsqHQjACjoe31BFgc94kc zZ13;T;L|h_MNTX{t# z=U`Ca= z&y$@xgg)ywbTBPnL(6_%xf~5Vm7L3%v}#h4huU2#!c{`4!t8P}J4aFKMLueGYG@k- zyH(lO1v{=tFdGH4PKBG5-L0(ptNl<-1RXW$)l*Wb4Nrvin=tHlY|3C?Dy45{PU1}5XwY!D?3zuhx ze;yq>3{<->&kQhI4A_T>iWe@=BtKt5T3a?4;H92dujdit(LK57r;-`=&akVYulcE> z$^R5c;U6CtIf6+IQ)&(EU(nd3K5N$rb{HQiaglH<5Ie-7cSN4S?G$jxHQ@Ft`=M&Z zzvy`a?ojp>RT#tGu`a;ASrF%nJYT6E*5wcZc9FDtO%=Wvj(Qn?;fL9xqRWv0GUXo% z^~==P@?L(Mc>2)by}G^pNr{%e%tzf~_;(BpVshz4{(an&E7;?hp znX1nZQU}V(DfvhVv%}gUrn2QE0N*1ZNWOIn&Bs`I`gZdVO0IKa_)w8r)&Z%cs_ z)a~?8B|j(kGMb~%UtA-9u`51~J2g*h_1C03<72r4gTHt&;qq}T3k#jIXUBDv&&XlJ zLR`ix{h3@j+6uZlT?s}EPu$}iO#h~= zSHVpmEea5nc99srPR`$>YJ0gE0NPJ!mnGXrt$($-IP#-H_qbjDuAJH@q2pU}Tktzs z3(IWS*isQ**P;cq9tgYi%yVWX0_2BA>^D5#1U7}ucLJkp!~WzWH`_&kEah?J9kd(mHceq=B7V zD}M^_5g4g!T6?v7)p&%E%4aT@c48Ui3A;k5QbJ|TkbdgT|HNyijcQUW)KH2GpA^$C zV0oT+q!NrnjMD#r9d_t^Qej&S|DYpnX^>JA_J&k{hGCYdZ3-sgHA;NCrk((W&MshiR4=_ke_v~MW@+d{0`gLX&HPO68 z8H!~>cPnR7ami3&FrmiD^-2$c25Dm2nFJ)EbU6jPG(ItVd8_s_wY&Jh)rI9~CwSaR z70+5zXd9Fy+dt~q9OFTXee$`zI=R2{q=A)(>0hO;@rMC@l{`OALso*v5GhpewQJSp zTjuO;6W|6umTxc#cvpFmi9buvOO;hwEj#y=HbhJx6P?nY(%-3x-seSc6{>Ri)sQC9 zu)9U<&~Ih%z&^@UL6_V0c;#}-ZPQl!0u|pW)0AivQ7&3juc-5PndOA(O#E335TLzn z?*9h#t7w+cl=k*(DrlAb3VsWR7s5NgOMV^b zOVpov5>cR*k{f?W{;>0Q<$ZdogNP3oO(blRpWKyrBJG;5a>v@DqakITz`{JE}8N&jiZldsc(sM$DwwJX2esU zwz>F%vK3^(bfOnY{00S?u|ufED`af_8UL)llqKYcMJXhjYXD+RD}|T@wrrFm5K_~~2!iQ< zr-Y9{f^I*X1(1HzPTqijCL?0$)$-I=inlP>P#Iw@vu*JZi{R)9^^5pBm%m9v2GRUo z^$m>w91c7AR-^=Z)7ww{HEawDAJ5Xy&vpLVguBa-qmkS>oS3)}IzXUBOp#IQGX*9H2IXTSXoHeAUI+mWl|iH6SlGO4j^8CS!X zd!EN)2*)w+{+T@HffMO76tu99L%Ee_p>pOWpCAB2JmoD!6@e(lgr4{wfk)+rJV!NPD$gh-3JRza4!b;{yo)vXXj3(LZVlvJc35)ugLJA9;`WkuKH;zbEd zQn5;9>Y;WqtJ!a%;xY|_-_?VpW)UuI&3G?biGD%@_P=AF0ps#ViAM(QPg4?mkFcAc zM?E;qnO#jGizBX0i3%l78Y@9@h>%A0?UYO8SR9gID?@LAeuc&Aq+rVB{}IX%)oRcx z!jLbl5W7#+wFy)zBDEaRnh1I2gsr2H$nrxN=@CtVw@n@m@AK@Pp9M1E^YJaeikGz9 zljHJLb>xG~GYyEpH~z%LYehwuW{MT=KOfV`4iF{2iK^0<;`~CKL<+3SYwD2-EeHRr zsICu-AVP9bi{?;lCxv??$x3aeq8zqrrloe4XC+x5kodqD103#~U>T3hr!P$<9|^m{ zz+@kN*{^a7?=l*coLoEa#AX!HQij{#i&j!}fFgW`)LznP-r9E3ZKUCroVB_oJi$X2 z7N@Yj7Mq8rI_aZ=-E+a=RCpiU^h+YWI^X~$|C`12Z^OxNl|VO_==_VgHpM-(4VluP zN~6@&^T2HBcJ44$XU3iwNn~d8`>0$<0+LHagxR?0M)ScR@q?0F1$wLFGAh=v<6P6} zBm`06mK!v#&~b$r@M#F&&wU?GZYU>=L%^5e$wkrj#zgiwDil$xm7JvLtsLK^me4}# zRmfS#Q00VofR?%x;oxBiL0^hm06~Xquce*0`L2=}xMOEGI8Q!d3%?cki}?Fg<9@W& zl)gGyNk!Uu!4eNnOAj%&e#yulog_tgY5&+rd#3u}nk-lG4vH^TnTdTa&!mp^?di%y zdKZ2aGot%Pu!fDhUwk*4U6crVi9w&$tj`OYrFBBxgt|l$v{d|tW&@hU7v}DDq>7UU zBtXl>++F7WL}aw|y+q?SgF*ek7-vC)P|e@u<3wbZ9Vk-0UK3z( zW5sz=8$n{uIfn~citISigeO{n`xU|4D;E?x(vg|dRB|=g<$?MGJ|y`lp{;!2!Z98V zO4?AKJ?9qEOKSzZL*?d(Tbx>88YqJ>{dZ8k+0|o%maj%E8Kk1W4q2)*%fCtZ-qnX3 z#5(gA_X^@QH4s6yTlgp!k+rt7Q*_R^*q|`kri2(BTv*F$8_XOv7n9(C1={7K+FK$N zt{dOZYNUxC6YxuaAgr>XI#OBwvT44JNH`}$v=JXuFA1cOYv?nb%AhUavZjjDABdKx z%gInU87<`Mw8PPs+=2(pD)5T5Nobv`0lk9Z@oL%`Cv4KM&=iNJNJka2V^RvsTa2?# zA&awepMBfRnL#p&5Q;i&^6KRKmCJcDrA%|oj#JDyWVl`CYiZv}jMaX3*w+t_qV_Ox z<~1muRgU5_V%K@m*bnsj0an*_(#`S`uJZLR^?C2;mNs`t+N&blC-aVrvpTDJToShh zcc9en$byL9_8Xq#X|#A|j26O&c3biSwBCu^b!UOu@VLj_k2*gUVRrs^@^v2M%qlJ_OBG}Ffo^#`(lt`N!hK!OpAb-bB zXxyoY+qTi(C?wFwwv7He9W{#cQl(1>5Z^|{2G^D1x*RH42aLOJpOQ)LhKj1J;{hQa zb#3Pxiapesx4l&@6$(qePm0i`cDxu8{&d{V?yFs2c3%MBRACVuRsN)RQF|iKzq1Td z^pX&}V8euo_oK>1w9c#T~oOpA)i9-&U;ebnZ}RL#9YaHKi>L}pA} zw38vw1cQn0iuyVQKkS6jHgEj69ayws;haXt#R`fL)@z27Xk1et0>sbv2;CtLZ18M< zU;!vQ*W3lOEgs(@I_-p&C9FfCa``s3T0G1mG!iC?nx1j}DLAK#HPJ^xokY-33zA`! zQ~GeYM%kKkY$_GPwJ^mRAZ2RoaTyWCgKO;^)V6GzTrR((o#r+ z#|=Iz@R*cT(C8wMb8m{6{$TH$TLaQ_dZ1AM2gyS8&>Ef zE!d~jBMU_-2z2hs7ca3;zsN3zTY{i?K9(KH@~3ICOh4Z1)*pd&*tXhXP^2N{yb`(; z7{~ZfFYeqqhLmlxQB$c$q>!NI_okQ_u{`XXYQ-1r&!xKSnD3NbSrr}yg^4ta-(rxh zZygjNIdD#wFgT z?C#S7U0hN!154P-Tr#t3?X2F7z}4zJ)3lv)$wA|fVzpsE>6N+mpIKx>hLN4Lhw+l3 z)CU@OnFz4klPRZb7$G|@ER%0j4oErK;W1+kToEd2!iu=AE$}Xi%eO25y$j1G98a8_J)H%ksMO9>o-lO2dQI*lA7P? zdVECRj$z@IZM9*ZOh?xtSE6v?(n6`Y2pTZ; zFv~_=8xE84JPNc(X{80`R|~5St!t;{(&i~yo`V{POqHS@ns`EJlpF#l7U^+=uxMsE zD~E84Gs0R5csV}TLrF0=Iw6pgjHj?b?Ozb}b}DDP0J+Vq6b#FK;)A>_SubL7kWO=& zX|d^MOmwtQUWfLGksa7&jZ^?JoAkIqfQ%(NYLI%8txHMg!GBszD^Ujw8T@ce9Yncq zP8l5Ho*`)T5glF5C*7JZlL`ReP!;j25q(WXnS%WdeZZ%bHb@lf$Q9bxs~E1O*i4CJ zJ2R8kaa%Zv*GMyH^qrMNwxC$^O=5l(UKQ92nv zMxB}~Y+(+wF-Z8o`h?{6Bo5DtxB~OjJh2-Rm4SfnT7)nix5ZM1_n=W zk&rkyMR?E;(fr=aEf8Zxi8zvfv!1w*K-g3+S68|=i5#@)Ypfh#`7FAv;wHfYbvEbQ zG(frdB1MjxmWhl4{IprOeVvBL1D#;Q9nZ+469TC3cAf=lQ4L#|NNaTbbOCAoYLQDgqQ{VM8((SK z=;)cSL_rharU0kaYF8`d%tg%=A3X>V;~?um&WqH9ED5oJ=%_LZwEkG{*5uBil#U;v zDqgAHecr8jB4~KB*m>ZiMv}NmWLRAQy}T(|a0#oMs*a(=2U-65SA|%I4f60`1@uv+=(Gk>!w6tX&zF^`Jha zW9#;)RZc1C9l2HhowcOo(aK{Kh4P3s_vB6IF99PF1DI?MzpOLjr>+SK(9i6lO0K7~ zW$H%<%;oY`_4ari3}_0Az_8IM5yU${)Jn8J1an=7;ALJ@y9Hs=GGV6rlrr_fKoUtFu(#B6HJL)EK|-O*Zb)gKLPG_3W;u|)S34kvN*&J zhpEt+3_=EW3yH{8bZf)T59}5+ld4SzI+fAmN<;3e0+XDd%+O5R>y_U8WQ@d|35l}d zK&UkZwamYky&g{YQf_C0@;22#pbiL(TNd<+5z^dzKs^pfzz!6)bgg zn0`C^EYWj*3~G+m?8pwb!yEf4iyhg=!a6Vjon1thpFB;b4;iBD0)E4oa{bceybzKe z1-XtZFM6wouguM1ozcD{`v>ZvhnK-REL-1Hey0#+el-W`vm?R=8a(aOA+kZr3=ZpT z@H?{loqLX&%K52}aY57Yn4yHm><9Y2E78X9Ph0y7cz4Q<>;`_$O`|n_?Z}=TqUH+7 zQkcA>=kO%=g#GXvuP+(SdxSPRQZI(o5b1+dZ~Nbq7tj}+9sZ2tA1`{bS54xoWk2nJ z#p9c{ICk)BCY!UL03cOLZKoc#$EHir&xcIq32cs%g@LEuov zwQ@E1wlI9tq0fCi;a zqKC)AMRC6#4B)lElNa!idr4KgBhGK$=rQJ&eSP`GqpFJC$@*$?!zGBH_XD)XK#ctl1gAmC# z8BZef+q0zd)b!Ku8&cDBAib1o$pBh{2Kiy-&Z4+rdh`CR!W8D(%CqrDOWorjSbUAg zm?<)Kk1of1wIdwG%l*;jo1)nL@iO-onw9p~{SpoiEY8BiTinc{5ne)SGur%O@^TD` zkmx*v2b%f7N5-U-$Re>_PT4MB(Re+aC@W*#Qsd&E=DW@wr1!dq{!>Ucl$B>YhCvxEy63}iAIOlAF-ASdGHvK7ZvBC7Keg- zH!?Y%ICVVf%}0I&&Q4O0C*IAc=;!SF82&_=JqWJAxvuB%n$b}_px7ha;dr~>)}P+5 zwbE^~hebdT6C znYfD_DCtBXU--HYORN~Fo?>xD~WbgiaNvCrxE?XlA_%$*M z)-&&L+MIqn?EbNz>+R}ufDYHy-mu~huAn)(u-Dlc8=`fVf3Han zx+JIHlR%ibbo)HFBalUrB|o&;sx*YI*DLO@>s z_LPohnnU@)*d)`7gEX(NMjT~sq|wj*DY{Vdx$S=x++Rr6zsPTm@XEbR@%h#1*7%yv zAIJCd5?6lq<9`F+t@dVd{k4xqgx}3Si)zb1{8TcEwIb8zRHFsBkSC?*mApi$O0!^ z&SUs?IhN};U_{^XGQ9}CG`F4XP8m0w3a+J4tqV1mqg+$1UR2+{R0e9N@K<_D68x6i zyawbxuOI)D#{bWsaQx6$IDYOR^4bVHC$sT0CLR`NaO{@-s2e~2zFw)WM-=SA^VPt+ zu?z069Y4R1n%*lHrY5~iLirp*xz-8gpu>dH`@SCLtH$>#Mp84Hd3%qK3g>mryrOY` zu{VD!Id5(i7MgKY>!XGbyDg}TWRU9a`a&tcnN(AtUgOtu)jgC;YSq?`M)WGAy5leZ zjpFbcW@(Fl6T&*vYPtPse*K8mg4>MCr{5g!g){dRR%Z9Qcn#NxRTXdfh*u|%UB-+1 zUerr-2V~ev(_)QQ|CEprF%qwxYs5FBWT{w>P+qv^0qp-+^mi?^M$^XdX zk$p_?j`RAC`1nM6Rk}TX&5z^b_Ojg@#|*C7&#%RJ_kUO0r24FGrOy4KrC*4rTdId` z=lv&m_;KE20tebMS8&n44sUWLKOIz)e!s?h>{(roqwf7AyN>nR=e0`-%xkXku6&30 zpCt5VWeQaNC=TqfMC6sxnc*Y!j;B_=dg{@d$Sd=jRoy%~GiyXA?Vh){)7PX|@c0XN z_a#zAiO8Cy_tO|!$WS7=d@(A4nE!D?9H##$5sGwS_|K{EQ!`{1q zNmX5W-(65bqb=uDxVKr1^wV*+%HWk4%ipem-XcJ9dz6R3` zp#}>`OpG&=M1@Qu?U=`O!6q;Mt`}?m`Rb35eGMD%H zzUTSgex9y#&OUpuz4qE`uf5jV`|K5~Sx*>}u=7xdK2!pPUbKguAFZkzva{}zmtr-z z^93zd!#VnbQ#a&|y0({`wmW5-ayYEsb~MzSQl`4n4NhBwG*Tg<6DRjEhU$}s2;_Rh zUw`zKn6C@TII=t3b+Dd|z-*9f@V(GbE!U8bA##nj?@9ooHfg-@+lcHM=Fh}@FR`g{ zSnOjO4Gf;nmwgS@Y=@8!SG(a&-q~6gyCT@Ts%^ARBLYI-|4Q!JN5409vtP2_ZKydk zHwBCDBSOP^U+Yd^bCdIn6~r%SBE0$9YMm-61pJ&tPBfn~ox7+d1~>RNcwXWSM;+s8 zo^uJGH2+$A!pVrs$efT`Fv0_GDya^-=JV^OmjQ%wY#*054p3TN;kFSkI|2^b5p6oE zPS^cbK9tUVI~+-Uv6rLlwZ(G>NL@K8mDuXCo0|9xX%6jNp%+`Jpk4Ee=qZc?5VJJu zJH)fRg6@}l>U9=Q8ewb8JP1<`^^k8*KB&G%;6?4%fSV{BAinUnxj)H_W;|-Vg{Nwm zH6%M_%$zPRYdBEo_hM;IPJU28C*rg0!YSV{@2b960y$ zt$PPXQ+1_Xd!z6l-y~d&0U~>iBDIK9ce1`jfy0W45oFekc1LSg%$3(Sz?tD^+Iy+% z?336Jw)POu=OCCzA4y$LEU%#H0bLm)VqCyELOWIG1>Mo=^?mYObx-c`CVi9Z^nunt ztXhJDt(c2p&?oir&uWoYZyQY5FD#|?ic;Ap=Npm)6oyi_VE>6qM#8hr059;hq_rO< z+*GQhKa`&QeMulCmM?Q9nG8efp~I}VJPPWl zfj(q2^ZVux!!U<(a(iFjPQmslV2m=}%t>a}$yNbBXX)6?#|_H`P`#L8i~^(4fD5$I zT-2Qov-P#83Y&YD7vp1uiX{n5*2z(5u_a}py$^W&T`zqgcf|kAa%jYt_>jCVy5v}$ zf&z<dRC=fniN|?4EPiN7#T~i zCn`&ivuq7*cG{-kVg#To@+EBj#U8oC5Fh-MOD5H4&wOAhw2ji^)iH_j$`%ktt|i3B zm+r01(@Xh^3r?`U^ok5dNW_bTxwS$^0FOoMVT-k%N2Yf(g80BEn0Xs5C)wi-ua2I9 z$u%@=Z2J4lQEC32SF%lN5p`qnfl;1Da-10EN-W7BzoZ-UjF(j+k|^ApzLrrGgU7{w zYC*pBt$Y1ldsD%SL`o&+s6c4L`cxfX{RQ`9>S@u9X657}xr`w3$Y7&I#G{m@B?@UK zzM?f7=RT(nj!M-S{(JbJ0gywCc+CWsOV40D!=LEG*O<(nhgu4fMmbKB)8Z>cuTFDH zqPjvnX1xO8E#Wx|4|Ep(CB^u@6Mw6|4k|+hmT&epc&rOO0~h^hbEF|>iEkvnlIZO^ zsBeup#{AmK%i3!#8cOpW~{#vg&WQ)7E|mbFgC!`?f9zXfpECJ8j|cB;sfzCDf=trE?+svsWr z5`2t~p1GSb8ZNa4^O={m@`zRS8F}iPR*R>j`sO0j2rJ~~4Kai9fk#(kko1M};fmm1 zqWB?w70!oZ0*gx20g>}`;{?^Fu3~x{spEsfsYvLQxXdGbYfrhe=>JC2b}2t$&=}>h z+^Oi$``migO^k0j$S&=xuij9P3dLt==)rIK>lPB&C2?)(AJtJhgSs-BtOyC0m=kS} zI?*){^Hsg%TeXo5z7w&U3IYyR4065`2@U7DzUF1=lu=(#NO2elJ%EJLf@uaC|z^q)tD+HQUrEQJhGKFj?> zEdy-S6{~5V`y``%OtQOmt*>>rj|gPVp9wp^$c6ZN5j*EasUCyi7hiIPFHL;XM-Ycb zb_Akc#65w!#@MmzxKJC1xEsbH> z|F8psq1vyAv~@Xc*GQ-&*A=~g@Z*+rD4g8HTo|oB7VhFubs|@Yd-WxBeR3TiH|N*W zXtaEf^tJT!cK#3Y6`V9D?~CR&j9w#+J}GtYf=HIG4@$Y=G&8rpX5HL>Pi}18(R*Zt zuX$P|xpREVUs>JJ(9|AIjV+tR!9Zc(^Z4~p!uiBRV+>lD7&nP+Lq$_W+~lzTS&QSV zb;*GLNI3a?tZ7fC1q0WUV}JrkahAO=tDlP`>k&^uV$_UGKeGqPyBa++-(KH-3q50! zAafG1MFp{hC2A?>LC}~tc^btDJ>HJfVAD(^?kPU&@O)gDshFp5a=ok1$z`4{MXsM| zo)TKlvrti~Hp2P^6ICfjiH0LM&-)EmrR~v5qPui!N)5kF{ z$UeOaoyODO!>lm=SqADeC@1F6D;TI7oMn%YXa;KHg`T>Fo1Dd(&G1(BS*W&%iQVww zyoJiy>RQ=!mOX9S&MgdR8r5f^?w~_Dko!A9SFSr$o14gkl}fa)i;w;d752AQ^YREN z(SG)-pkOpk)uupLNxr~Zvsv+*PSx&qn( zdDpRyQv({BcA5qHAjvTTxs`G1(k|Jc+J%sgRm*kNX~w9>PZC*dhRfTb!??}cq2@i) z9#!m6+rDIMO>UUGY*?mAh&&@|6rWt>1GPoI!IDjWu(cIS-J%JfCFF@ozL z+|!WoZ2l&c5Av$z(1&d4aFcmn=^_HO4To0kUh3) zPL*Pu2Up^a7@1xMV`u&Xwm99S9G#B9!kMdi6>Cao4(KY4n7L&%p(~BdIKS?DB_+OgH6h=oq@Ln~>Z@0=6=PT4 zyc}Pe_kU^5`Y>4%kq?`p5xGql(hwva{n8OJdl4mD=WYm_X=(i8VF+%rj1Q)Nf_X{M zc%kx}??UorqhULz?dP7AStiI6*aMun3$i(X2{7$*?@>fTH_?ly7a+IvNgRY|0o4i!q$8~ zd7#d>zpJ;kmngosP&JHNpPE|Q+B?8mXiF8zHIdYKS(JP+~^@T*2RIKdV7hQ#rf9Q}QK@xNcK_a97(DuKb7f08^Yu!Ygp zxFY@3_B=;G?<89VdebkwuBo-<{Dy}-2_Q&D;5@lVNDPKvr_p*sDF=M;Ev4fepnbu=s;(5e0 zBzD+fT;+bO%L%eB*VjAIjaW!#T{60z7;%iSE~D6fvMybFUa3_L-wJF?o-zGbINw>( z;Ol7c9g5X3h*{N$3|Zf{*=ZXHbu{P*K)fNEBorH~xVW_2_R6|UvQfpoxb{1kvA*t6 zo4O9fd`DSe8E%f>j@{>Q;t8|f%k|i|zFy-Gj3cu%C+F?V$vHc7s+KCWjlbOQnk^Yn zR_rLP-F_gSOLJ_w|3b@99wkz0P3PRl`Ia9`Gu-;7nVO$Nl3{JGfNi$0u(w4sk&lf` zJ`_;nNn-9Dkin@>4(SHH&PzRiXO1Uf#nf%tu$JDpiOOz$@)a$3u})4JIA@P?G%fkd2UGhFs(v}SCm)_G`6c*yZs%?jN$)V$>U2sWra zcOD-8OIQ&Ip4Y%q6gF$UZ!MW9)?PQHBhvPq({>cq%G=VInm5>EfSMDv&JWi_Yc^1# zJzU!3v~5x*@kDZ`P#~gL8Ms6Sv&@u{kzAi_Z>V0Iv-2_P-5R;ywa(-nR@DKGJ$rSH z9p{wmq_04JOU{KU&_0`b%=_Q;QQT?!XX;KKtX{)uTT-S(_SBSfSP$@K6dCkrOnHMo zPwi-vTwOT4GGC=U`x*0MnNe1-ORLhKe#IMauP$RMUqa5jKB!=eHj7LcPE7~{LGujb zL7CsMG^){%8IR%03-*B%HShFX-z8U3wYcjDq^;?|>QDdXWl#2Jvg9-SHM-=SkAX)t zIlh!lM9$JSoB3IGJskdblvvhs*VPfq`&ndTxMrQR^aj0W{Z({tAi{ZS#1NZ@0ec=SXS-Zbmk0Eb_j|*yF#`x%)LZabWTk;eN_FQuj`0&W4gv?J-5N*yTSQc4qs@ zbWisHdV!a{IxF38vO8L{dhYhxBWngWk7|93H64Eqz+Q~{PEe@D-wL4h#DPHbr7$D$ zMJ4s>sF@f2e-~ymRJiqcWy_?~sIlxgf@JTGeE)Io5pK#Q44&VbTv#%oW56qg01DZb-eB90NjJ^hqdY6J?wm9^cuJ9VE7`uttCVrVs=ib3Qam58>9JPtN z3(C3M%Os!MpFRH;R^;b~3(wJ}uH>my?1B<`x;muWj9v)F?Ar^3TKP#_+_fmA z53yf2a#iWxr^Q!jNrfY=MTeL|M~6A+&O7l2WXoDp6HXi+$9G`mrQu|IxFbD4hDg<{8#QqI>4RnnA-&|wY;F`8a|jc7^^ z(uv<8((YZNMW9>To$~8SGyhI6dwK<1)o`a<>F^&^WmRG4w4@ua^WNRf-hB7`-oI95K~y;}A2`HPGlWg!Id8LDk9V-kZO>3| zBvb)f>zwlIQRQlRbC7OTDV7m&Pck$<@gmunC@H*xJDky0;Y~HhH|%}bO#w?KFhd1p zzyNiVU(iQUu`hUu0HkW?|K-API^Pf|j3xX)HkL4;a${+w(c}+fkabObKXEI35Wx0t zJZEpxnBCvnddk;)FI!0dA-O+n?nia+aCM6dKM9~m@%q+-l476eq%^je!QK5CquM!_ z0gMHFddI5Bw_j@EIp^lz4YzliHexmJ&iyGPsjA4{=tqv+c}7N$`s53z$2;EvH84cF zp$+lQ|C#ONIRmfu^mym*L~`TWX}b>1Pw~#HZPkPL>!{eAO__^xJ}^a4dP2Mt1^8`t zd8QbdYa;}r)C3#jY00(P8$RWGICMETFS3wD22RMUi7q)ex`ef|E2Bv-p7;iHsd6Aq zzj)$D*fMcOJn`kH#S@PKR584!EG)@K6knW+C{BGjeOhtF;lCldxZ;d8H|)2#;$VyF zLCnY*=a;bM<$QIP#yGzLI{kzE->5VK(&^VnKyeC}J;_hxNxzC(#L2Ult{_`r%9Cd) z|4WtsrIn9VcfVKd6=Vxc`Q@4A&$S5R0{}9PEsSTc{6mzz>CB=NOLI+DVf^`LFMMmM zN|*I-I6=cH6#FF$4?A@i`Q>|-PyA8+m&QT-XsGx>V>XuG(En}Rtz>P`^4B@WC;)8d~B*21R ztRSh&z4+k%YBBQGJR{FLosnZp1{WBapm4^>PAo~7K^S>VG?ILoF0o+#WJ zY(?V~z44vXqHtf`8%fRxMX;wTL8hmG!yz~9!Q-X*| zKe~e%+j^tnK^Q##)jI*^kidof6Nm5Zh;JKrJ{JEwc56sPs=EO9X#~Q$h9yL=-pTgG zdFPi1Vnp!QSCgL|7Q8TNfs7EyQvF0ejnz(tv7}ST@8QGwoCb#ez}O=Iw3_+|6+6}` zRv73-;>C+dEw2lE7}@xnk4XKlB2gEDQsZpFT*pg-ZgKQwEtXHNPi{oJdDW7%C30mx zUUh2hZ!Fw2p_=mmaPwazXXEtx>Tc}-xQpfH53$+nLtq#s0`*$KB_zDw#feXy3rwxu z$}yrYD@^r+p#{pMy^U=aYd9B`S~06*yw&&8dZiD?PYX;f>szlL;o29Nx}B$bZ48Qp z?vD~IoAvPOZ3gVSp~rG19_5Pht5r#G5E^Y>_+{7?c9|ENE+5ui$(6pBX21u=qr`_TKEzgHsDi#bNIsTYS+ z!B0k5uORsIudEGh>w|9euNZY(+O3e<|``o^{K=UX*35hSU~FFhjowI zNsRgiRiG`LkJKf5!_^1sFnm(C3^E9mYIPVZ-wR<+0#mE(qk=AaNI&$g66L7~ab$ zyHOT{V6DdZx#5BSU+89>UtE@t}Z_EY!@6~5xM!XH#&Z$j@^SaQd$Wto200I%w^mL9IslWk7zD~%LG-dWX{ z@WnOfFm}ORae|veniHp;#W10;RZW;Oj4fGgf787qc+zi+;N=klgFgwdJqOTW6~NFPs$sbB>faZUVZ~h{ z4XY2KN&})j$(p$9)ZH!tSn2w38Haazd~qhoKHf+SV_?*g>)VrYH5DWuAcSc!5C1CJwiMDTu+1wvLt0pTLEP>7*jt{s5EVrA_R)G#6d8?myE9*qC>kp33Omuw|<9FpaG$YT;dL=MI^MG z%7`azzmMAQWA^)${XT(N&8Y`99J%2vtHO(yd=-Hv*;wmMV8x-K^Y+nuDQdnwmLspH zvH$CTgPKOW`gFAFm%wz+X=<8p+i8j#G%__!6f*umJwr3=pDVAPku1e+35#xnscWbT zq3=;ovk2MatKP;v3hyQ?7}FXaL0j|Yyu3xIS*xc8jo4XsFZ+ZuKWAsKMNEDN*H4M( zSINv_107opBc+8To}yZk<7yc*Kc(Dz*vGle!noABvNvlx?iZpD?At-%KKuA&;C=1m z5$5-D+Lqwuk$wDETuqJlgqM5Hp|9ICxCqTLyreueZh!(~-pghMvhD8j8g})TR@& z7DFK-enZ1c`1j^`Gvot9Jfl7Kq@Rg8XS-HBIVqJZuO>mtejq|cmV197zx1BX8_5A;($oeugVf^D|K4IaB817=3Sl!(Q_*psW`h4j3KH=1gD_Xf$gCGV>x!cp^B8Qj}$d=DAs*!sA(Joz>w=;kI{FxboW4BZXrh_B>Q^k z;30zLi8?^$oD*=Ks^bmPFDDI$`mLUq=z@x?`9XRXc{luW9?aC#8zvtK-cK2k#3&I4fY0MMi2jHS^7><#y$ir{niKLY=Q&Z9GQe1%EOMgk z&pSZm9!ai|c{e8SCSm_G7U5*C&a#X0?iHD@)8%Ja01x4mAeczGTrJ61=q`5Y$Z?{m zqbQuI(Wl5#60Nq6EPDi8z7zZ?zroFtbD_w(CgnixiQHK#r^#K~h=L5mmmtLHz=u=8sz_=Ut8m!-a2nYj zG}p_x)RZe&@R+5KSs%~)QX;R$sSH%krSmleSNCpM2h)`GWUtHJf;dC1w9xV-^n?M zhY>%rhU+bkkjFizQ#IOU#lvz5;5lhmwP`i%2}(|)gd8T;OVSQ27I}UX)*71L*KHyu zAS__o4JJ3^)K@j|; zR|P)e%QaURGvT@QWyV7}%Hb3MNI%4i*ue?Rz{5bnM?uqbxflf>@!BiaB0`-qa2y9c zL90JffR#h8ubm$?NM19IoK=RPx}@C$m?xm;s(^xJ@5gI>kaJkFChd{1Of`)xKo^(i zMgptU<)tXdd333wfH1k-ex05Jo8zLK8>O#P(2dp2F~MnBGFX@oPPc& zK!wXw0(W!Yuu*HFB0G5sO8j}i)+j*F*QD;Ty;eV!1OCO*W?jF#m&{n^s^8NK$Q)W$PEVv!eQK2JFvn$kSoUu0WyAs(fDQh7;ua5YrS?T9i&trTsXin+2(2X!UV$8_KhM zMbeJ7#*e0Z7ELxF=SE)Bh2OM3S!`1MY5lgy<5b61<;;)yz`V%g^n2~($$Xg4I13&p zOFF~j^aIKjd7QlO$sBiDoP~kmjyq|LUlbz0| z^{8xjIm4#%k4inYisT!j$q93KHHJBEGMgKxvj?^4{7hOqi3M^t`0_I?%Fom(POwAC zF@9#p`I$L`pZjJx!IzZGjvse^vomZf*9q=n9jZ1D6l!y&CbJ`lH&ZK|VQZBGi+6>d z6HaiCu9eQPw9-Qk2bVfQ1gbk>x)WT*b>;|XnB}cwP1gVXt$YI@;lM zH7EGC?&mncEnMf1afY!c$-QqZQMq-!pIE%z&ak6S;J#8^vtuB`{Z3$N)Zy%Jnt#`6 z=j~E^6^&vRn!qv_*mZW?K(^r&D}{0+i;rB>_eNCFN3W?Pc$_vm*YYYJX7T~=gvuNV z+#Cslxn_=Zf@}@qaOPYyE5K@Uk52Zm`ukXE>w9j0v48(4aSLU9ZF1;3iX_9;$-ofUeAAh)0GX2;!C z-Qm>S6Xajfae$Q%eaSTIDmw2)^8a0h`yIk!A0e5B z2D4#=0+YkN#!E9oYgH*E0^NuqI3pyGC;5LDu&4MxEYSP8Q9B!7(GTJ&4ivT0I;imt z%D>71m!uzbf`@rp#a>Hi*t<^Po*sgJ(lcR~pSxVDyU*3os*=X+{?%~Mbnma=+#=}o zF2^N1o%?U5{MR&;67wtQV?C-`_v;nPJbk#+`R2_ubq~VA5B|ti1Gq7IY3eC;OBq_FAyS+?%}9ukjMCdI5&bI(v+AD>M}T$j=YPs z3`ivK6UEaMh`|Nw{dQqO4e&za^CHf(`q3f%h0b|FDH(N+*SFT6iqfMz01@YtgGi#H zs6En{V;*?DQE1b<#kGyk(1)NL7ng40{utM1>C@-5PNY6m?`@#u4D~6aQN6oXC6TRO z4*)cqbK3}IUgff#|C3G-sTOC*D&;jYoN$6$w0Sr%iO-W(I}3F_EEPbcDmIZ< zt#3%aBT`Egu^^i!=?&Pg&UZ@PtMpeUq{==;_3;-d(a^MCw=l}n^-UdhzUQPCPfRsm zb+TSd+)Y{DCRXc*ZMZJ>WaA4_WFHi71R4t$li0Xg zsQ8eUpVgBSrgW z%03HgUP)L5945bF8}>lvqQ-Y8LKEhM$iPyDli+0d2%D=BXBS8Xpv#o23@Ldg+)E+e zwYU{8rv`VhQZ5I5Nr|_VGD<&}@)Lhf?`n7_9ye0PMJ~65r0q)bsOSks**5WgLVp#6 ze4QJQkej&+eAH+Q-#;U`y9le0hZwk2?*l*-&kBL~7NvL;1=@PS-9r+c4yW%Yhg`XV zPXYw?x}1ztulUKAm^MXzlWe(hPwsbBL=s>WF}*)q%8eSB%-)#cW~Io zIb>jHa|QqqxkVZ!=*^bPAe*ViEf+e$4SfF!fq0|YB+gW@39>S`xI6hr_4=q&?ofr4 zGfBYmh7#z*D|Xt#khk{fUAHnNpn{=y^`duuXb7LYr|vw_(8R|$vqI?xTxzL$9|SrQ z?%0PMxgYiGUOvMq1B2Z-2LR+=5&NRH$<3&&j3wLgyxW2eYkhnvfd_)&Qm!eKO5Q<6 zm3Jdj&g(Y8$yBmRdY3d9&j~k0nwD*pm@VOT19qP_V01Eh#ZAk(v)fcI+ECA(bh17C zr_rnSxN6yW&Bvx6+}I)%3Ez;J=zo9PCYXZJUV1%4=(`UJQ_4!keU)J0uK>esGFdpC znl%|Sx2=P*uJ5jMhHcR6FG}Z3#CsTA9TY>6imxRh4IUk6a2-fAn?KeoTPHcXsP%=d z{ABdANk97m8@NS;c#eOq8YNC5Nd7`8`zaTfgzM4`nA*fiZCadal-~9(xg^38?c_i^ zBU@UOVSk@4F-XbFgGx8uQMp8fMo~mMsTzp9A!RRjh(-t`J(vP0CJ>%CYaqvQ@{$mw z1y$>5|8vTSW(e~2rnkjVRsbMf8#UgvjuVh~gD_pEN~D8`JbWobhpBuc_KcKGIH2)G z<>MHe;5sTY1)EgHRIrD-m<*!oNq2ggex*^l6D#z4Q;<5I*WFi#L!M}tRSt|8HV=bPIPGm`gt67i67&o6s>h|!_CS`y(o)DmlwL^*jYIFimjz~r@d)MOJ7J5jy3}eNt|r7WA*!AgOj#pL)u?J}ef$`$XZ9K{6befJD-xLZt#b0a^KdIMzq}GCFo#V8|{S@1!^xdk)#=F$G zTmtcjBF^2r7#QdZumcPoqzzU31o}8^0ucgLR3=u@MV!G!97ca{q4St-8l|K?%avi? zD8`@<=z9`H?8c3NUrVE~6s1Ns!Va(vFT^3w!IBZwL1`!%X_5tk%DFl*lvC=H64(fPoY$y8G=PY0qlfurte+l>3bA>*6c#mFy2(drr*99 zIYig46#}LQ4y2ig&78&vb-mF^ro?KV#9LA(J+(&N z7Gjr>sfcblx}p$3K#9IgA0@d9NZH)4W=W|t3ij9hKA}w7&=BBbK+_XVnrJRC72AqQ zD>Omdh@w9A;-7rM@GU5Y5cwV19;IL_71*Q6Q&=FnzK^WMAWyUvd_8 z@dKx+i={>CVv#=RUOfFK(-*S>TeUYc_Pz8)4%RCujH5r_&~z+(bDi&yeI#dRnXTg- z@ptQ!$G<=Wht?n6hE|l~+>y~%2y9zEh|B8J5KEJ43(M8PjpkqL&?8Q@oyY-<0{n5X zB@d{Bn5w9}XoRbr@|*cYx_kf=J{=+DN_$k4#0Dm`KL~(C-@&JOoDq6HU&FCAA5FF{ zH7j&Vo-Un%hk-Uf+_gWJ`jS+*b9y+E_Y~*eul9|48THTnifpO8Z~$nY&!;K#KFv?; zd{o2cPxHeT{X3Vdn=euOT352KZe^*SwV%G_I&`x=rRZnk3;LYNp`{Af7>^jz@r`)M zSfc4HWlE%*Us3;~lI)`2+jrzValBt4vj(4T%bXtGNC93eTEHQaB+ufgVhV_h)x`@a zT3B+h(OcMJ;& zyk(k8Q0;0b*eN5ioWVoo>!G&Vt$CsR+~Q%q+&d2}nmpuv(tC6AR5ZDz2m&SndK0N2CS=s(J1yp-4`0q?kJ{x#6u^W?^adN%Ds zLVRuJiP+)=UA`3_Sf$}^ z9vyhg|9wEkzF(>#$v>zAP+vC8>zPBO2Y?pcd}F#y82eZU=|Gb5R|em0Na< zyP3CJNZE%EL31W-C+)B0GSJ*cI!<6))gW)n7bIDZ$dTo0^2H&_M5`K=7d&6&wWcCm zz<{@{)XKLPxch($x{cj20ut&mgM%o}PI=jl_yYVe_BiBSCyk|dc+#YIV`YP|3IRPg6gJUZXprG;OXx8SL8s+q$b7~vd_!<_g} zq$hmLSuBU`CIgMELWM`x;elihV*bM2j_*v%VlRNB#@%Wuf*6yYD1Xjrg!fjNlFiZ@ z{+m4#j0t9uwPeAdZb^)y2D!voH^blK%>`=^_fV3CZINVw9GY!x8?nS?I~}xABOmXy zb?a8%+R{NghJXw)>a2R5t|ktpt>{TE@f>343Z~%x_ZLZ ztC&fk^xvnSVe$Kmb|^qDELuLk3CE>$yfHvP-UUQ++zS=evN^epGj@{OyN*|H)c4WB zMMQ+#+Hqol9{8p(Jr96h-uch({mGrnEjQm@oh@$(`!sTSw|;R{dmjeKwQx1|mWQL_ z6Ybre`>sE?p8FZ-HM^FdY4Uzeh={|*u+wHBd`%F{oX&Gxs}$0`W`{`Ms+4x+ujJJMqHggAx2_CPAqp?nE*;;dTV`CLQeo#OE*F|E@g7Eo*{%pVZ(TsEa|sDX zjdjBbqgj5Jgk&3Ji_6Kwj1Lf7T}r&_RigZja%}_qX)M*_x@Cz8P?Bx(l{ey@FwKcI zQfk>p4ipJngM)k|QbMa(Tq4-qYD6=Xni2Z3xhq&sYGmv+G8xN6r|I{(;B`vGnG9k^ z+8}kw+Kfv9c~`GS8nNHxcG?tf@#d1{nH*X>$Yn&uY|;(@1~(*mnHl z;f5op7JAHgS)m$j?RSPbCZM!z(?&|aig?i~p6So>;IZ~1cw-<8qPOuC@T5_VxY{s4 zX?{W-OFhjYaf5r`JbuA0N1U-_vxScXl-sTh^Y0htOrTYVb@P_0z@^8q;K{g1jvq~< z_zSp^7NIpJWBC`9$vSXW=p-*RDaUQm(QNA&JkOM+a*fiUamn&yCHDzUN6Q$v{s3Zn zLBjfKmAD6`fdAdJ$VQWlDvc6lV`OTzkRFaBHK7r+wM>`9A^zn#GIy(7uv({zYu!RM zi0|nkx)jweFmt-ql8`1xErHp0DRBlrQyEBz;`mF{GD{=Q>rriwfLyWAJLH*mKQGy+ zM#Z9<76*t&E#(gI?W#yR2n$P;P+vNM3G*ls3H?qPn{7ry$}{L7Hg%~!h8^pSY0HMQ z?{y7lph?%J{{S>ga)A!g+PYk3bbnxGtu3(76v-EO8 zV~#M~l#sBX3&Ev|^p_|hxhL1YBou`u1OpR7*(3%gjpj!~B`|-ik?b+6GtIBo znAMw)n*%`vjUdL1${L?Q*O-ytIgXEk`pjQrjbT>;k^)p!uebQOMtArsZLKdbm<*tr z(1Iw)flMhs)5_3Q0S(zaaX!IM>eS(SP9~1Wp3W3Ne(x`O^2AX?YrAmZoc%v;KnCyFrm(P+H8U7S9Lfx0p zoEG0DM4_e7)2to@LN=Vo>MSx!F0;Xb_c6{h@O;Ry3p0Y6-Dnb4cDrPNdtVhl^Y2nQs!Pm{iYI|T zi@NWd1$wUtff1sB>^gIX(^UzZhd~>FCY_}#*nrdKCb2zo1VU?u6=90jd)d`XqJ$?R zfNJR)cRq{~pKooNV!m*iwv}Ag6r)ul^mw7J4(dWNwMj{y;vPm7bfRTWIdTpFptkR%%_ATW;PkFPY4)X(}<6P_0nDm{* zqy+ED~)a}8(6SZ5{fivtp*)KcNns8zUCB#?_)lLRSSHs8%0 zLt(@aXeTDk<9_Bm4{vwU{mKZ*O$4&81Bvff3z-rdO_P;r64!#1R0n!q9OtQFy(F?g zqJ=)zaILp@?Ra=44$ zn-5BI2R;+s#Ielp%!`k62eq-5#4*f9k5SYX0%g;(w8;<3UpMa}p&x^$DMFI0>Ooua z^oE1=#U#cK0{zFlvPF$F`82%o6sWStfF_qrU_u$SnCQ0UR7sA&Y)|Xg3t1)s!b0rhIr>XbKwOw!CqmzhpzMyuHlfvHvGHZ@P; z50oV9KtkLgFI-^K2&0wZ#s^Wn4NM9tkJKCTNL0~0DV6PXl2L_T5`v(_I6-Lv1?GP$ zGDE2YsMAJD%4fAe@#ac|I$#Yl6`;7AfCL+ccGf* z^`P_zQx>IuYA~!N_YSK9^4vuN!z(ap7AZmzA&pLIKEpT!#xr3|1NG?#m5WN=WQREgCyQjNIH4~|qLvC!cR<966Olv&(G9q2o& z5OK~hJ-u8a*ReZEE6egB;Zp0|`CAsG_x_ggNCLZ_Gfw<~d`O*#0*sy)Y_{J^jsRj= zqo2J&aeU8JIJz%5B!?8adz{WwTb<4{n(hkxDr}DBvR?#5MP&2X+qN~Ln+ou?bguX< zLSNj7D{B+pKTsfnolJ2_FE#}LlZiq(@VIrp!Gl@tJY2vM%gb!#$mhI3EA3h8u~Dbl zg+j<%$q`29py7S{qip{= zpSpRwbIf%yGb9mCPfi;Rp0=@Ks%BYulktZ#7_$jbf$_(7bH+UpL^js}MDT3=*7D$(R!04K372*fdHz&gjNAhnUqq zLI9W}P-ei+oq!^2A}3-}5qlhgVwR)4rJ9DAA>R=QsgCY^Nn)FfS5pTIkl48-?4x2+ z1H2Tv2_LPw9%TN;*(Ca|Mu}}_(R3ous#&~=X=Bx?ni}}pL3yW+5eC}fg4zUVp+Ma< zbe}TKRHg(IFpV}db(-Xuyh^qRyqy~+Lx^Qc;Yrmhp2`t42N`;j&gh|m05+Rfs*G*a z)}oQN6sq%Vb;xrPNDHxzgHqcA7v@B|%_JgqU>d{^%r2v_2yR<}po%_;bP@!3749?4 z3&4bv6K+!QJm&(RFsQVPQU-wCsw&NoxqY+|!S=ohY0n5kzY5F_N;@K>4G2}mZf9wo zCY<1WjTp&8jT&!jD}G@V047Oes6@Euk(!H@k|TX=C_7%~r0E%8lk{ z>e%A1IatrM_@->^3*VxT<%tQM`c-KOC&z}TG^A#fF1VS~e@b_JfrDRhMK6e@ep|v3 zGO^UZDDu7HYUpzc)dl-#>Fcn2~pc6m?+f{4PA z0Vwq^?Vz6a<$zi14y1pX1ICvY2hzLnsS2gsyBB%nYq_e0?Q?I5E*Tv%kl>j9QXaR1 z>3gg^r?^r4r(W?|qtJo$)Ewwa#m3-$pQk!AkDS_HD5#+8;W|};?A79f=>+fhKGt`{ zFO@>y(SKav(bs~W=%|3I^rc=QQTk|x@}e#$?xwd9b`jbM9=*+ZK>)CZ;L~!!-=H_@2I`=CG8Og4AIR&i8v++n8|d`f*K)uPrhiC= zy+1AbqRks`-r{UMk=K(iC_K$@Hp;@Eky&MPV4|NOv-bbi19fOwsgt zj3l;Jv1Ign>L{rjMZ1Ynz}cQ!GolJM_ZIPwA`FB%Ln!%LLvm+Wao3M~2&6X=zli8W zRsmrb0Ilx>^dVqry>Ioxm*;bu{h9HEmB!Ab`s7-}=!?lrZkZxHEi;Gmq7Tcn;ag&< ze^ozO1Rsvgx6K)BMxMn0b1QqxJNVz+#;1;L znzjOy&}I5Yw6sqYum$N%OK^@ZN2zQICu7`q+i9ervr$(3Q^gqZ~IT{$N0Qot3$5O@4eew{R@J)Anm} z_p#N=JD1ldpHCi7u9aLXWsNPjinbuVijohU&`6N5B2d}CP|{l*bCz=Adr2|_^BNDr zKt~8_7m(X3jV6!gUJOIjdE=g$6LjI5#XHhZNVxaI|Jv1se8%_r z^8Z>WPw-+z_&NX*J{EH`U$f)%YhlfBlA`Nt>UgPvE@^fT!>ojIgG{x_^Dpb@?Tw-e2I8te-l1>jDTBg zZdm~~itf;wAtG<^Z=;WdIqRpollU(ztNony zww|-~7pO(7!Iwx;b`mtCy_Cm^D4*^nF<&R8FjP?Y*sF#)s=ZBoY^kk5E}L(#o3X$6 zGFB_^>uu=c8E_|`u4+%F540LGT3V)?_~1D`DjkI zZ@t=Z;0 zNch5%x`o;MBg6%FmMPwe)Ang*CP>XEG3uRlyCjv>4DF^7#_Xo5D-xXH+rraW5wZPF z9Oo+A+)3ydi^lV>uM#Bu6l}JljMedv)G@cCL9_oMlumRYFLV!MFh2sAzHeu|Ehj2< z(~Ia$9X%YB#g0jfoLrXv9*j;i$}-pGb%TtG>Svm zPH$9tB6Vb*u#s!@`S9`Z)Z|9uEHWk)RB%JME9=D-#Kd|0Z;oaO@Zv1}HmW;qiZ*TK zsILJ|+mEHnH?dV9xq^_=4NcD;-4;VJUaL={Hb$$rGrIXo-Bd6orNgQDC}M}Asi|C8 z6JS1!I=@)Yt5;$ia^EfEY#{-qZ#Wb^x?P735+aP#q+^^6<=Z_mHNs(VpL;Uqdu?I6 zQlH}Vs87^8H>_<)O=Jh$&W7ak_0_z7DVAI#c&i&7m-xrd zLw1OLwB|sq^P|ooha#7B#A*(Dr8oo2i)h_h+jhul`;KPln&$u%<@A35t|(paw9O+m zoSNF(dYni{gGEZL6Lz$3<7_??dVepzKQc4;1f1ySchJfeG14De!-o*8xpO~NItl?atgxJ zayiEFaq*QG(&ZvsA99w_%3Z*ESuFh3@<>ODPQy5t#9=bAnq4EzT>_X z>^v1re=CCF@oF^Qwd_zMA=nzl7)b@c8A*0m%TFfHGY${F_h8qP|C`(mmweOuK1(I? zHZd9Tk8dAewh6qz2k>>U1$<4efcmv@pwjvT76P!AG1?q@T!TcTZez6RsIF_7Q&VL& zKI&2v7}PL0U&n#dd6QC+ZKM(e;juJ2wn*1_&KYeTh{hy z)a@suzN38A=QmmFHOXtiUfPmtNdHEN5tvm2B@N6X|CsrdI011mYTc+!rnhE>Bcwcg zGTHs-h5!8#tW$>o^?TueACA;)c9wmGL<|2L``$?o;$1BEi6l4mMgPTgw(zsvC{D}= zia;BaJx%yeiaKpoj8+=TY2%LDKm)hZ4@^3?l4QYu=Uz?VUjqL6yY{k;{?#%sUQ{jf zjXRl5@9m|$sXPNvue(E}XnB;!2wTCE0Z;9a_NnKA2GjE!d zHWaylu)d1qH{}9O6VI2pM8FOl8I$`KAvyVA4)HHm=b=X`EJY1%g;Awi|1SVkmbva zv;GD*>Y#&co-@kf_U%J+?_%-w>2ZAXeK?B)b_!TxbUgU3G*I4Jqxm_r12dmD{v!ker%)1mX*xSW;8Tn$#vQYX(&Z3UtX)tyD4Er`rJYV#s+>z53 zqAI?`hbS*YFBqGP=Q~o0b!`pte9K@*4a<>?mNLcjJ<)Z7cP=q~dAq^5op`-sAC9cm zT9oVmK1aovi-CF3w$d^vqsfDusQ7Wc=fHE@Sj2KBEa`auSHwT#nCGGT@jq;v&&L9u z=)kSn!7ca4k~{N({_s;&ci6tr7}_=cBdEb|LGKX(-IiU*ZRi0{Lq&#ZcYELB8sN#x zZyAnjs@1A0t9wdV#PIRy&W|3yhhTc0bAK!iDcW=jzYs?de*7fT>5tLm*iiKz>R_aE zykzkE02pA!{-8ua;$G4bUwCjOa&%tgSZMMsDOnxz5LKlQ(bc97b}|N>wx8k+1A-k8 zkeJ|qisNnc9N@4WQrVD$b9=B7?A9G*@0375V15(hl_M9Tz*6Gl{d+6){`684l;8F8 z4@Mlsi)6TA?P1F_ITCtx@!FFpKmkA2H ztb^-QWlZLLBZA_#RJ)!vqjDtbPqB|QeQ$!##Qnn@Vo%lF=z7jU?nEQh7 z6g#oiKL5|vGzK8Gw4Lxn&8D|v$z7Rt1W!>9LqD&F?~-$`(8IW%j-(gjMJ^^}o+f~Y zpcd!1KTXH`3G^$v4q|MS=1JQ##94F~2ykM6FPhvy6wd0EI(jXZ^5Xsy%a^FVD~bK9 zgIjTg5)qMVOS{sk`w8y`{Zq}oek9wNXzM{Asu~*W=+Gb3JJRM#$m8t;NNSctzT%*3jB-8ac)A-2Z_Us`b~J{Lk+%6U*zBCYTHgYZRbKKoklpk zKFK+`ME~8cr%KgWUww>Tzo@?YrTXd<4K*8_)>j!e3jQ@FN?h%fzdq(WKt^l~Tq6X6 z+c6;Akj zK-sg#O@k=T(NFjHM9~^m`Z-FkS&Vw&l^&2U{n!0VueZ_zqQ2u;;LgJPq&gAjC%LtC z$4TU6>l*`8K}6UV(Y1-)n7>#7J-ss?QuCG)0LFi{7^@^HHHLWF$8_?A22{!hlxBm$ z<5Fc~`jS%<2h?Pol#KM$#DNjtMjdXgC)qHol6u9qj&8TU5)Y}jJ{azg4B_0G|)e2r<-9N+Zr5?86i zu)|~3?TU0<{Wf19GPvG$en>duwU}{j7-LZ^>=^4heQzK9G)Y*h$zx`q&1aBno#_N- z`OtS5PO<709easGvZpNC@diPOH`iCcoYTL+lEZjpTfq_hXZY{rtfuyhw2>+>zT4rf zD2|HK=W!Ql8q_5)jt|Z^&t1v)TsuxCJDJAY_JMrzSanD1F<&$8lealC?Xu|6O`IG2 zBMmsH-f~{OuM14zLp}L`v=g?&9P$6FC0HIyg`MoGod==k_`CCi^Nhm=k&xxYx zC;bBJM4A6VCWPi6B{MzAnRyzO+sdWa(`g_h`Gsdd$;H7;PJwK{tiM1^f^0lhH_*W$9zag zv$?!;(?x7%z|?gb`96-1#8x)^N9ZJMNaCj^(@^M>_J(j%VfU z=xE)Iv!mPj#j~kfPIGq1fZ@~Ky9>^anDdL(_%i;P>$CPHYRer&QQOMg$Dq$UJMe}q zF=t0ZaxH!iy%}Q8j^Ui>?HmUK0)F{O8M&$h+M@FqPH8OJn zdWz@o2x|rTId2ucfiLt<^PdO3MVsT!gR7CHbxSx0z3D|fqTh|6!7TeWY3A=peD5S0 zm9tm}$48Ru`y3v&=8joD_(_*`=8iS!I)ESKe3Z!DI7DI`Zcu%wfQS55H(`*l-h&EeDB!me z8n`_s;u7haiz%9SiP)C@r|n(NC34KjhtNKrcYL66M3T^N5({_*Pr%uh43=v4@#aM> zt@S3in^R%tLg|=6JR>hP4ipvuOhN32>|T6y z_(ec1hpUXcd{<9lyH5dCEXZcoa;?v;<^0Gm(!ZZG6;xZqH`e44yOdU_fa=oqY%?hW zo#~-iq5q2U$@xPnLg`NS=ERm%gfcA*0;3NYT)cyT)0w^K-$fI)Hy-MfoQzwaJmRfO zuz2bDM}}ixL~{O-G4hXKnK=*st@t0SV~kXO`9}zOj&Fp6zQ4#DUJP&TVWn{l4@BNd zuIf6W_+HVb6Rkaiou&7wbQ6)fF|~h>8c1t*N0Vvdb|AJE72;t_1_Y%R<@_7NxXt-D z%t!0_H*BAsF*tc)^Ri)Hw`rHQ=im4iYNDrAT?mI3{2RmTn^vJgu{rJM@^37PRKHFI zASUO=!<07X#$csMDMrVtkc^b0AR^TKf6&YtRh@(Sxj9CrkLTSSQk0#xk76UMgt;FQ zgFH7!!EYhkK5jaoz4u%f{z2h2ile3lw|>3~)$SQrQLUJa>0(n!uocaxBHjE=l(7YH z#_8UQNP5V8UR}>CqBmx=HILaYBCiEqdITP($D--)zUJAqBh^$f_L20HuyhgT@|ZgA z`F~diUf2$K&r5L&1l{VZoU!5_RWSQnj<}8W_lFi$EZ;nqU+1vq_)k$`#bA76Li(43P#wGGN(;6osr$WF-np_xq zk_)3l;bLmK=RObqb-w+c_oA&Aj|ogypPJ;?0^h+--^nB}f*JSCQ|1amBMrd3cMl~b2c{cjfJe-Pz%(F2;lEGXX z{zV%4i*^=!Hq0sUzV3}ma~smhdG`jrl$w>5d!x4C-grN+#^GqqtI<@2IWT~+pT-E!r?J<38Y^->joN}w_pK zj(I)}hW`W zOfHVkXzbzQh_Fg&qhWB)o#9_GVB+5ROsVJTu*_x^k6o7R z$bkZAq1<~x^M`AmZn+j8$7fu89k@9lF{dkRVmG1A=nQ^8Qhgx$hzBt5>A-e3M~5fi z5&XF^%Fv9YCSedy@`q7GQgevmzP=t$Km=2%QqB$5_1K&|S{@qIqKb#3sar&a3$YVE zW1Di+q+qP3qaihcvko6x86$#bDL3tPLsmv>IM3su_ImbT{0Jwp7Q$=mF4>A-@2ZBH z%^V8)>5?!`lUNOhQa!Y?ZpgN}wvA5PE12j~BHRYwOIov!)odVFI}R2u9d)H0PTRke zMPdY@F)*xocTpSiQ4Y3BjqOs98r+Rn@H5>xznUMfhV&u2j+2meID5UFzpgVy4|t9g z4hF!WB1bKrlrMXp6dgY!Ps(_fmDkVh0IAC4ODrGD6?P+T(>Cx(_EtX|NsZ68_6}&i zOnwS;ZzyHYmIY8FXJKTc~|IQqMwv4EBI4btT1Wm=>gx$ zY&EBdw;ok24-5X3G3kX{Jn8o{q;}?1QT_-E+fZgg89}nmM{a?oNJ0*0>Ts55GRbJE zE`sE$G{a-my37G2u$6O(-zEXC61Gl3v7tF1i{?J$a-UykO#0M|%-#3#UV0jO8AkT1 z(v@Pt|MU0%3I%-H5?az;TT)WXKYNI8pL5MMGw$W*%eUYAm8NO;+;RJ?X?Nb)JoXE> z|K*%ncU*gVdeyyO`Nkd7nyV|nc;8o=@2d22RerAW&RMf>Z~knh_m<~bv+tf+Irr{+ zzcTagJ7!dVczII1b^Xz*nXWf3!%scL>Y`*LEyDR_l_WP#^VpX*P=jZhH%=y~v zrqkP^wEt>5Ad%Iqwp{HsrjSH~3itgQbr$^ov-F-Zd+wKK-SMXh{$HiXpIgtz?!Bh+ zp1-{Jb07QU=sPRt+;RI~cz54uy5i1z?)WPnuDr7{`i;t}JMN!DH&*Jo`rMMQ+}$*5 z?uD1z5c)EntSfJ_l|qMc1PiPT7z-ztl3S|zH&FOnoEos zD@#1w^PHjUj?a5fJFXjkMP$l}OFmfZJ$L8QOY+Zu?0ERAC$IXO?|IL2-pZx-c=_gR z@tz-Twkv*y5)WP`2WNOLf2poY<9|W!pC=ydmmky?{msumb6rq$U*fxH;AdAcWT%oi z=keZb5BCdx)|;(F{a(&>jIKZL%`QjK@8CL?Yspi+*&XET<~o`8BY)GI)#Brv-}Yu7 z;rb}o$GL9gsxNy=e%G6QhwI3vd$SRA+>x*JX1~jI`5V32G}n!Ldb3Xfcg{ZQ=lbYD z>cI|K-rbvR|*_9bDiFp&Eh;NdE)kL_QzZ+@5pBNaUFRVc;l0pGbfvUl&gDBHY*>+ z^tsT%A^mIM!F4(AoZoRRc?3LRmg!t8b^V)cb~M-JT%%&5zs+V7Tqpkn^>AJAGwR_w zl5OK7aiGolSLnxe`oDoc*U5j#X18)(&aU=-dR_~?KEV48*=&&O#%HtHd0e}n%VuBU z`ozX;wyGTbyRzA_T$gj5%(b2C46g3xZ1yDA>CZ#2AmzC3UhiK=8yyp)v_UHS?O(nkjUEf7zgG-mA=_$Q} z|KY#x&0eM)k-?%wNr?Z^{LlHP-mFW$lA-SSp~ELS=iOJjpycx(s`XRj2Qa99X9Wb!Ucka;A>(2`hh3!fz`G11{ zU-6m3d*xAG&-1_W-+Hr`aMhpE(){l`Elqo-N*Rm!d}^=$eBYQHGinU}yRTQ3lw6dj zNsM=Dzn?PI0QJwNQ^pTNI(+HCaYHLp{&7RA9y}*JbmXEz;h|$%2TvF}ePCeF zxkJa06dqbdmT^Oej~nWa8(La_p5TcVZ?2IY4ZLxG7QDe{hj-n1lYzC0|L<<;&3=h- z?7tunk4n$LFXG|X?i;x8oU96EJ51^9S%y zIdZ~3HxF~dS@6u!amS#6a|H9cGYC0Rc&cC8f!Xs|Z#K%k|I@%6x?@o5;6;NTJcoe0 z1D7FrOG?J_{|H&GJPxPwcSU=nCa;|wGv4Q-SLT2tJqOGyZAIzuxg zA;tLH&wKmc>F%pL*VEDoe{^T`-QCZ=efxIzz1{cS{hpA;T8S!MSd)w2Kj*HZH5y}a zR}04CkT2F2JlKj!zLVyE-=b(+vE`eKSz2(ePJ7gEqNp~By0^V}&V3%~=Hsb4)be;E z4I=$H2|lgZ2c9DX;d$9iV>r-27o%;*7Wl4gX+GD)FVJet(I1_2zm4Imkqdnj-R>pb zj#=F~*5dp7*5>@k?nc1x@`bK@57HC}s{oq;_6m^a!J~U6 z_9E=rH9Q)vE53hQIG+d2-z<2vuExt9(pEPyKBAjZ(oNXUzD2pkZEJdKl(OI7vI|g8 z4bt-6a0UEbnJD*bvd?{=Q-O+h{@AvyM zRD#*}EVEnuhW4Q!pqM}dvQ_gcjW>LH#5_F-Uy($7O=9|vh|yOY+YeG7ZrB5?>2Z3|-2vpdli|8^@@ z73uM}Yy*T09lL4*Gd+c7pF({HW+P;)_7@l?B!wox?_eO#y3-x$d7=^^(^!ws~I9 zSTA|q=!IUU8|o#@cT>XxOd8sF9&L>Dx^9tdY>?JW-e+O$hg#9cT7Sj&N7x^{CWNKH zkgRL5PVR?KPwz0+hN85_&Fu~KvVNzbNWK1+;|&Wovp;kgZM5ROpa=S5+u-q-Yh%v0 zinbkI;9Gv4{fzp4hWHI%Sl^E{=zE{^eGTHZx47;v)UQ)>a@F_!2MxQ3zCYZgfY$DV zXyX*#X|BJJzGE%?B&~S!9ZUAx$j>o?m;tzU9)i!bMEJ~qVf z74Ul!?}Ld8@pHiMq5sCu>BoEL4?p6%bRIx~?EN-i-TPhlCyZFE2<-9$uDiMh`!KKw zu%ck3xAy`&0;~_VYi*Y8fc-l5C;Di=+Cux4@5MLIx7h*AJ?lUC9Q_#f+^UZ=SoII* z)|bu9=SetAUfU<s8Sr@d=R`))7|%hAN&Ifge?VDJT7)9*cxC*PQX_rJr?X* zVEuUa_tBUjB=I){Ec(cK{2A6$_GjMH1{)r2?9UANeY?)<&peEF1|R*u^Jm&I5Ep;N zb-zFzYwFK@(E1J4mz!<=iTd4;$1@3DQ{a^#UI;YClX#r=U)u~o@6YScj5J)Ou}{KLZm_k%0Y_UWXw7~MXM3%Pc`x91-hT3DJm)XIFJT|>u4yr> zX(On57JPQ%|8acx!u*-PY$y_ZuVW1hdfA^D!4b?;-^4i{t>Y}~`Td#K&an7H|*kOf6Idn3$}3IZ(EA>{u!Lv zUSQu3oTu-Q(!h&Y`SgNm|6ME1*|y(V{M#$d%gdJ_b8y8jq(2kHwqi*R`=4K7e%SV* z#m}xVf3cEeVm)U*O%y-?H!{mze2`_APvM31O|*(eCvB;paR8hvhJK zNMhM`&08d)`4yi08c%+M@3jPPm3-ik%`cNYx#AQqvyo|+auLgD z`H#=J1tOUdU!bZ);#WvyPDwerV-!9v^1&=Q(76f)I=-(BQq4p@YAeu^MPW+n!_+w7 z^`Cb~VBe(>kT@u@ zMKAoxmoGKseLX$z3a+_nbRu7w2zIOwtq+CQZJ6L}=N+A0p>S8I;|eLLsWvZ3YSL={ z!}S)>ye~k>L7AAK5MX;~NKfT0CSW482P5>`$XtoI5&x@!6a5)AD53u~!2@bYf^QaF z?OtFrg2V!|%(Tg zrUHppJV*3)pdZYa^pANO(f14e>@X4tymuqsF1R_tu+ej^Tbl5pf8MYP<3BC)%lCDtx?lK2UbKAp z{Ts=9n((<#=w~BHAh557__*K!5o9%&8bj}Am!RE~U$Ow7c_H~1zynytReyxs&<`<) zz9{2!g%tda;F04N2qI4*iG@8#yMgGTMUDA8lrNzX9kGPE4t@_=;y-)LQt~-BlK(Dv zMC30DqgMqFUS|2oFqwr|@M%0FeU^~!;waGm)@UBrf}YwvLyAFkvLS5{{N&RN8$R1X z@_OcPvD)XYLT^OCd!#TS_+e?B?(ry$3qHNxw)+JOn-b%asfmsNr+!VXwBXyMjDCII zs9)a(zNF>bAWy1F-NNy**WrdpB^B>Zy@NZh2?+KB71USv(ndKH& zdnbQ~@EMt3x;LlrQ{iJiVe@>4gruu{uCh4Yi&IzugQ9jL!#2+k;S-htI|KcoaI@e+QE)z=0RCRV17Be>!{^vYCImNc zvUul`f;z%tukJzr~;`HK9mrRTd&Brg^` zKtJz6(D}Ml@R2TyC%Dv@<$?!ApJlksYT&Gg*IPpR4uZm<&`--cM*BL2or2GL_>Tgo z{ldG>7li)Iowh#T_aeE6>6e+HXJ7aLaFWNn5B{O>nR(6f;d5Rj?-P7l#^D89Y(6iz zS!>5Z&m12SJoo`i$amRDeob+5NDy?szNP%{v$)>HeiwM4#=o+q#{84gf8HXiq~HbM zJkB0_U$n^LL63gE131~~&<@MsdEv88@RK*#e(`w(%DMy(ZngL{^ohc|1rI!9^L(#> zWG`@*$FpvvnBLOXKKC}kKQ8=drN3&|pB8+o)5@d%!2N;;uClo9?+*bddAxS|w9rRB zYo*cq@qYv!tQAnwKV$m&1p(+ajpxXp+A?v=OuOI_a_|szzE%l7(_?Wyk3#tZUWb{q z?s0|CPsu{g=XszA3w}~wtaW_e1w2^C!Da@8-aKgYqVDE~CiJ%peRzqLlkXc)bwcpj zCoCVnhePtWfYUg8^^M*;X?%(@4l3tSq4(}bj|qL4{!tQw>f4iopZONUM*WBH2!8k) zORxHDb|gxbj!?lkPJA@uRN`YMIXdR>r?QD+(d6Qpd?-=KCJpGyl?nW0GGUUbQfe$y zuB1v%CFdmYjdrRG2FZdmmMx6p!|xzQj^1HxtU5mX}VG#gFyJ+w>;+X=eo!_PVbEyx5k{|Z`?E(qZ-8$zBEreiSgSb zeS-s=dYn$Y$aa}Tx{~fFS4_54n8+uqMV+CmO)_1sIB8Q#S4PVdqb6R7=R4~NoLjf8 z$LGoG@p*Db*nT5VgFKLqMLS}<&1j}l+TEEbrlHWi;^kWXQ-2nBTtsb;bzW7m`>@^%Vh%qEx2 zcA`QqetT*UiR3tm$v916Jd4lg>x8z_c5Gl3;-Q=D>fpw9R64_sQyDK6_Bh#8eylR? zq)MeiX=qEra@e`kq2Y#DJ1QNW&dvP6zJoza4P3AkvR}#G?R7#~ojx#nn=@e6?a)EkX z+3V~MO9LH-wHc$5&5R~6199*n3##&2-&LI~NMOplO7-%WlIaRD8x41btFl6DJK%zl zmYiqKqMgvEW>cDu^>oH^Iv`ESl9tSlmPwFQ&aL@9Iopxy{3NqOyK|AAojW&f8gM#b zBi9(xRqg@lXEH4rSPe}*4N2m5va9;x#axqZlQ{EcJc$`k%8XlinA7mq9yRa=l~he= zD3WEi6&k`!k!2&b*C_)tj5p~uchhFm82-6b(@{#-+eeLV@d#KZP)XRRRR(DadXmZ~ z&%>6dkd;40mGbZL<0#ttB-Y_cX>mt(Jzd z!%QYIlxWBfC=%2*uT0I^_jEwpV5qsm?)v4Skr~x!67`SzOffBKwtuvWA=OA+RTG+V zstd!%GB~`|ncJ+OGPG3<$MK67_LWdEN zWm+LAcx;@sl2$oPYT8?NO8+Jm+Tst8xHJLovGL5d?3y3ZD zAl-qGHKFp}TqQn=xKgt5xTdhZrAkE;!f&Hep|Si#s8}jsNvrIwfkr1X+2p!RQXrW4 zI&u$~d^}zrH=*R-Jesg^rDRKXV`0n`@^u*pWu;U$P8B3o%vMZ@Z9ar#Xsm!NCCaIU z2~|>)$dY3K>H=?Zq13pz0^><=(;S*fB%IV_B2}!wd7+7AuGUMiJT#5Q>krM|}DT4JR2+(zLI zNYvIJm--rCu7ncH>q!+I@YL7yb&Z2nN}KT3r*oK};Lb$*ulGtCYx{UT;NX=%4wUYc zwEiiPUt{{;xD=?5G|pOI?|(ADsEXD%+K|M0Zv_$#)(h60F?3u-cS2fU?}0Q938}JG zI*so@Io&a7UhkPSo|XC}KaHKrFR(czxYpNuFOBuyOWW7_+W$Y0`kREI-lJ(eB^{>v zba&{j|7l=EqWtkHVuQq<>k#gPP;hhv&)WJysi$#3+SfWB_U3b*`t&;`3K}cDMdyy0 z6rVvsZTk_auQ3;DWyD%f<43&p;0|Ld+q1zQeWp^ z&*e2%JJR~z_W#;bU+-r$rr-Y6257Jf+U+e2VR#55(W&Fr~ zDQLSo5BCEjNwvP-dq?P>IwG9p=iEc)X*zMh}9 zOMNfB<*K!fc{>2Y!gCO(wEj8|uJvnUkBV`PerbP741n(Ge-VEwkM`>#U@ALTl*bct dPyVb9yx`id(vrD1s{dKKs6$xmso=?){{k5-z;plr diff --git a/include/tap/dsp/analysis/multitone_analysis.h b/include/tap/dsp/analysis/multitone_analysis.h new file mode 100644 index 0000000..b60c037 --- /dev/null +++ b/include/tap/dsp/analysis/multitone_analysis.h @@ -0,0 +1,229 @@ +/// @file multitone_analysis.h +/// @brief Program-weighted multitone metric: pink tone comb + joint-fit residual. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap's test harness (tests/support/), promoted here +// as a shared measurement instrument. +// +// Why this exists: single-sine SNR is the worst-case metric. A filter with +// transmission zeros at k*fs (the compensated designs in kaiser.h) is +// deliberately optimized for a different promise: alias rejection weighted by +// where real program energy lives (predominantly the bottom octaves). That +// promise is unverifiable with single sines, so this header supplies the +// instrument: K log-spaced tones with pink (equal-energy-per-octave) +// amplitudes, and a fit-subtract residual over the converted tail. The design +// history is SampleRateTap's epilogue chapter. +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "tap/dsp/analysis/sine_analysis.h" +#include "tap/dsp/kaiser.h" // solve_dense for the joint LS + +namespace tap::dsp::analysis { + + // ANCHOR: pw_comb + struct tone_comb { + std::vector freq_hz; // log-spaced + std::vector amplitude; // pink: a ~ 1/sqrt(f), scaled to peakSum + std::vector phase; // golden-angle sequence: bounded crest + + /// K tones from fLo to fHi; sum of amplitudes == peakSum, so the summed + /// signal can never exceed peakSum even in the worst phase alignment. + static tone_comb pink(std::size_t k, double f_lo, double f_hi, double peak_sum) { + tone_comb c; + double sum = 0.0; + for (std::size_t i = 0; i < k; ++i) { + const double f = f_lo * std::pow(f_hi / f_lo, static_cast(i) / static_cast(k - 1)); + c.freq_hz.push_back(f); + c.amplitude.push_back(1.0 / std::sqrt(f / f_lo)); + c.phase.push_back(2.0 * std::numbers::pi * 0.6180339887498949 * static_cast(i * i)); + sum += c.amplitude.back(); + } + for (auto& a : c.amplitude) { + a *= peak_sum / sum; + } + return c; + } + + /// Sample of the comb at input sample index i (rate fs). + double sample_at(std::uint64_t i, double fs) const { + double v = 0.0; + for (std::size_t k = 0; k < freq_hz.size(); ++k) { + v += amplitude[k] + * std::sin(2.0 * std::numbers::pi * freq_hz[k] / fs * static_cast(i) + phase[k]); + } + return v; + } + }; + // ANCHOR_END: pw_comb + + /// Fits a*sin + b*cos at a fixed normalized frequency (no DC term; the + /// joint fit models DC), returning the fitted component's power. + struct tone_fit { + double a = 0.0, b = 0.0; + double power() const { return 0.5 * (a * a + b * b); } + }; + + inline tone_fit fit_tone_fixed(std::span x, double freq_norm) { + const double w = 2.0 * std::numbers::pi * freq_norm; + double ss = 0.0, sc = 0.0, cc = 0.0, rs = 0.0, rc = 0.0; + for (std::size_t i = 0; i < x.size(); ++i) { + const double s = std::sin(w * static_cast(i)); + const double c = std::cos(w * static_cast(i)); + ss += s * s; + sc += s * c; + cc += c * c; + rs += s * x[i]; + rc += c * x[i]; + } + const double det = ss * cc - sc * sc; + tone_fit f; + f.a = (rs * cc - rc * sc) / det; + f.b = (rc * ss - rs * sc) / det; + return f; + } + + /// Refines a tone's frequency by comparing the fitted phase of the two + /// window halves (fitSineTracked's method, on the double work buffer and + /// without a DC term), returning the refined normalized frequency. + inline double track_tone_freq(std::span x, double freq_norm) { + double f = freq_norm; + const std::size_t half = x.size() / 2; + for (int iter = 0; iter < 4; ++iter) { + const tone_fit a = fit_tone_fixed(x.first(half), f); + const tone_fit b = fit_tone_fixed(x.subspan(half), f); + const double two_pi = 2.0 * std::numbers::pi; + const double predicted = std::atan2(a.b, a.a) + two_pi * f * static_cast(half); + const double dphi = std::remainder(std::atan2(b.b, b.a) - predicted, two_pi); + f += dphi / (two_pi * static_cast(half)); + } + return f; + } + + // ANCHOR: pw_metric + /// Joint least-squares fit of all tones at once (2K unknowns via normal + /// equations), writing per-tone fits and returning the residual out-of-model + /// power. Sequential fit-subtract is NOT enough here: 24 tones on a + /// rectangular window leak into each other far above the -120 dB floors + /// being measured, and Gauss-Seidel over that coupling converges too slowly + /// to be an instrument (measured: it floors near 48 dB on exact synthetic + /// tones; the joint solve reaches the float quantization floor). + inline double joint_fit_residual_power(std::span x, std::span nus, + std::span fits) { + const std::size_t k = nus.size(); + const std::size_t n2 = 2 * k + 1; // +1: a DC column. Subtracting the + // sample mean beforehand is WRONG: a finite window of pure tones has a + // legitimate nonzero mean (partial cycles of the low tones), and + // removing it injects a constant the sine basis cannot absorb — a + // measured -48 dB instrument floor. Modeled jointly, DC costs nothing. + std::vector ata(n2 * n2, 0.0), aty(n2, 0.0), basis(n2), sol(n2); + for (std::size_t i = 0; i < x.size(); ++i) { + for (std::size_t t = 0; t < k; ++t) { + const double w = 2.0 * std::numbers::pi * nus[t] * static_cast(i); + basis[2 * t] = std::sin(w); + basis[2 * t + 1] = std::cos(w); + } + basis[n2 - 1] = 1.0; + for (std::size_t r = 0; r < n2; ++r) { + for (std::size_t q = r; q < n2; ++q) { + ata[r * n2 + q] += basis[r] * basis[q]; + } + aty[r] += basis[r] * x[i]; + } + } + for (std::size_t r = 0; r < n2; ++r) { + for (std::size_t q = 0; q < r; ++q) { + ata[r * n2 + q] = ata[q * n2 + r]; + } + } + solve_dense(ata, aty, sol, n2); + for (std::size_t t = 0; t < k; ++t) { + fits[t] = tone_fit{sol[2 * t], sol[2 * t + 1]}; + } + double resid = 0.0; + for (std::size_t i = 0; i < x.size(); ++i) { + double model = sol[n2 - 1]; // DC: modeled out, in neither bucket + for (std::size_t t = 0; t < k; ++t) { + const double w = 2.0 * std::numbers::pi * nus[t] * static_cast(i); + model += sol[2 * t] * std::sin(w) + sol[2 * t + 1] * std::cos(w); + } + const double r = x[i] - model; + resid += r * r; + } + return resid / static_cast(x.size()); + } + + /// Program-weighted SNR: total fitted tone power over the power of what is + /// left after subtracting every tone — aliases, servo FM, noise, all of it. + /// + /// The comb is generated at physical Hz, so the converted tones sit at + /// freqHz/fsOut regardless of any clock offset. At these SNR levels a fit + /// frequency must be exact to ~1e-9 relative, so a converter's sub-ppm + /// settling residual is estimated from the data: joint-fit at the nominal + /// ratio, re-track each tone against (residual + that tone), pool the + /// implied ratios amplitude-weighted — every tone rides the SAME physical + /// clock ratio, so pooling averages the tracking noise down — then + /// joint-fit once more at the pooled ratio. + inline double program_weighted_snr_db(std::span tail, const tone_comb& comb, double /*fsIn*/, + double fs_out) { + std::vector work(tail.begin(), tail.end()); + const std::size_t k = comb.freq_hz.size(); + std::vector nus(k); + for (std::size_t t = 0; t < k; ++t) { + nus[t] = comb.freq_hz[t] / fs_out; + } + std::vector fits(k); + joint_fit_residual_power(work, nus, fits); + + // Ratio refinement against the joint residual, two rounds. Pooling + // weight is (amplitude * frequency)^2 — Fisher weighting: a tone's + // phase drift over the window scales with its frequency, so the high + // tones carry nearly all the ratio information even though pink + // weighting makes them the quietest. (Amplitude-only weighting leaves + // the ratio unresolved and floors the whole instrument at -48 dB.) + std::vector lone(work.size()); + double resid_power = 0.0; + for (int round = 0; round < 2; ++round) { + std::vector resid(work); + for (std::size_t i = 0; i < work.size(); ++i) { + double model = 0.0; + for (std::size_t t = 0; t < k; ++t) { + const double w = 2.0 * std::numbers::pi * nus[t] * static_cast(i); + model += fits[t].a * std::sin(w) + fits[t].b * std::cos(w); + } + resid[i] -= model; + } + double rho_num = 0.0, rho_den = 0.0; + for (std::size_t t = 0; t < k; ++t) { + const double w = 2.0 * std::numbers::pi * nus[t]; + for (std::size_t i = 0; i < lone.size(); ++i) { + lone[i] = resid[i] + fits[t].a * std::sin(w * static_cast(i)) + + fits[t].b * std::cos(w * static_cast(i)); + } + const double refined = track_tone_freq(lone, nus[t]); + const double wt = comb.amplitude[t] * nus[t]; + rho_num += wt * wt * (refined / nus[t]); + rho_den += wt * wt; + } + const double rho = rho_num / rho_den; + for (std::size_t t = 0; t < k; ++t) { + nus[t] *= rho; + } + resid_power = joint_fit_residual_power(work, nus, fits); + } + double signal = 0.0; + for (const auto& f : fits) { + signal += f.power(); + } + return 10.0 * std::log10(signal / resid_power); + } + // ANCHOR_END: pw_metric + +} // namespace tap::dsp::analysis diff --git a/include/tap/dsp/analysis/sine_analysis.h b/include/tap/dsp/analysis/sine_analysis.h new file mode 100644 index 0000000..63d8de0 --- /dev/null +++ b/include/tap/dsp/analysis/sine_analysis.h @@ -0,0 +1,116 @@ +/// @file sine_analysis.h +/// @brief Least-squares sine fitting for THD+N-style quality measurements. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap's test harness (tests/support/), promoted here +// as a shared measurement instrument: fit a sine of known (or tracked) +// frequency by least squares, subtract it, and score the residual. Used by +// converter quality suites (SampleRateTap, RatioTap) and available to any +// primitive that needs a single-tone SNR number. +#pragma once + +#include +#include +#include +#include + +namespace tap::dsp::analysis { + + struct sine_fit { + double amplitude = 0.0; + double phase = 0.0; + double dc = 0.0; + double residual_rms = 0.0; + double freq_norm = 0.0; + }; + + /// Fits x[i] ~ a*sin(w i) + b*cos(w i) + c by least squares (3x3 normal + /// equations) at the known normalized frequency, then measures the residual. + inline sine_fit fit_sine(std::span x, double freq_norm) { + const double w = 2.0 * std::numbers::pi * freq_norm; + double m[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; + double rhs[3] = {0, 0, 0}; + for (std::size_t i = 0; i < x.size(); ++i) { + const double s = std::sin(w * static_cast(i)); + const double c = std::cos(w * static_cast(i)); + const double basis[3] = {s, c, 1.0}; + for (int r = 0; r < 3; ++r) { + for (int q = 0; q < 3; ++q) { + m[r][q] += basis[r] * basis[q]; + } + rhs[r] += basis[r] * static_cast(x[i]); + } + } + // Gaussian elimination with partial pivoting. + int order[3] = {0, 1, 2}; + for (int col = 0; col < 3; ++col) { + int piv = col; + for (int r = col + 1; r < 3; ++r) { + if (std::abs(m[order[r]][col]) > std::abs(m[order[piv]][col])) { + piv = r; + } + } + std::swap(order[col], order[piv]); + const int p = order[col]; + for (int r = col + 1; r < 3; ++r) { + const int rr = order[r]; + const double f = m[rr][col] / m[p][col]; + for (int q = col; q < 3; ++q) { + m[rr][q] -= f * m[p][q]; + } + rhs[rr] -= f * rhs[p]; + } + } + double sol[3]; + for (int col = 2; col >= 0; --col) { + const int p = order[col]; + double v = rhs[p]; + for (int q = col + 1; q < 3; ++q) { + v -= m[p][q] * sol[q]; + } + sol[col] = v / m[p][col]; + } + sine_fit fit; + fit.amplitude = std::hypot(sol[0], sol[1]); + fit.phase = std::atan2(sol[1], sol[0]); + fit.dc = sol[2]; + double sq = 0.0; + for (std::size_t i = 0; i < x.size(); ++i) { + const double s = std::sin(w * static_cast(i)); + const double c = std::cos(w * static_cast(i)); + const double r = static_cast(x[i]) - (sol[0] * s + sol[1] * c + sol[2]); + sq += r * r; + } + fit.residual_rms = std::sqrt(sq / static_cast(x.size())); + fit.freq_norm = freq_norm; + return fit; + } + + /// Like fitSine, but refines the frequency first (a few iterations comparing + /// the fitted phase of the two window halves). A converter's rate estimate + /// converges asymptotically, so the tail of a run can sit a fraction of a ppm + /// off the nominal ratio; a rigid fixed-frequency fit would book that + /// (inaudible) offset as residual. Tracking the fundamental is standard + /// THD-analyzer practice. + inline sine_fit fit_sine_tracked(std::span x, double freq_norm_guess) { + double f = freq_norm_guess; + const std::size_t half = x.size() / 2; + for (int iter = 0; iter < 4; ++iter) { + const sine_fit a = fit_sine(x.first(half), f); + const sine_fit b = fit_sine(x.subspan(half), f); + // b.phase is relative to the second half's start; predict it from a. + const double two_pi = 2.0 * std::numbers::pi; + const double predicted = a.phase + two_pi * f * static_cast(half); + const double dphi = std::remainder(b.phase - predicted, two_pi); + f += dphi / (two_pi * static_cast(half)); + } + return fit_sine(x, f); + } + + /// Signal-to-(residual) ratio in dB for a fitted sine. + inline double snr_db(const sine_fit& f) { + return 10.0 * std::log10((f.amplitude * f.amplitude * 0.5) / (f.residual_rms * f.residual_rms)); + } + +} // namespace tap::dsp::analysis diff --git a/include/tap/dsp/fir_kernels.h b/include/tap/dsp/fir_kernels.h new file mode 100644 index 0000000..70069b0 --- /dev/null +++ b/include/tap/dsp/fir_kernels.h @@ -0,0 +1,168 @@ +/// @file fir_kernels.h +/// @brief FIR dot-product kernels: planar, SMLALD dual-MAC, and channel-parallel. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap (include/srt/polyphase_filter.h), where these +// kernels are the ASRC's hot loop; promoted here so RatioTap's fixed-ratio +// converter runs the identical, already-measured code. The performance claims +// cited below were measured in SampleRateTap's optimization campaign and are +// regression-gated there by instruction-count CI on Cortex-M33/M55 and +// Hexagon (see that repo's docs/PERFORMANCE.md); a change here lands in every +// consumer on its next submodule bump, so treat the measured comments as +// contracts. +#pragma once + +#include +#include +#include +#include + +#include "tap/dsp/sample_traits.h" + +// No-alias qualifier for the kernel hot loops: without it the compiler +// versions loops over distinct row/history spans behind a runtime aliasing +// check (verified with -fopt-info-vec; SampleRateTap docs/PERFORMANCE.md, +// hypothesis 2). +#if defined(_MSC_VER) +#define TAP_DSP_RESTRICT __restrict +#else +#define TAP_DSP_RESTRICT __restrict__ +#endif + +// ANCHOR: opt_smlald_gate +// Dual 16x16 MAC (SMLALD) for the Q15 dot product on Arm cores that have +// the DSP extension but no Helium — the Cortex-M33/M4/M7 class (e.g. +// Raspberry Pi Pico 2). Gated off when MVE is present: on M55 the compiler +// already auto-vectorizes the scalar loop with Helium and the intrinsic +// path would replace vectors with dual-MACs (SampleRateTap +// docs/PERFORMANCE.md, hypothesis 4). Bit-exactness: each 16x16 product is +// exact in int32 and the int64 accumulation is associative, so pairing +// changes no output bit. +#if defined(__ARM_FEATURE_DSP) && !defined(__ARM_FEATURE_MVE) +#include +#define TAP_DSP_Q15_SMLALD 1 +#else +#define TAP_DSP_Q15_SMLALD 0 +#endif +// ANCHOR_END: opt_smlald_gate + +// Channel-parallel dot product for high channel counts (SampleRateTap +// docs/PERFORMANCE.md, hypothesis C6): history stored frame-major so the +// per-tap inner loop runs across channels — contiguous loads, one +// accumulator lane per channel, coefficient broadcast. Bit-exact because +// each channel's accumulation order over taps is unchanged (lanes are +// channels, not taps), which is what lets the FLOAT path vectorize at all: +// its strict per-channel double accumulation forbids tap-axis SIMD +// (hypothesis 5), but the channel axis is free. Float-only by measurement: +// fixed-point planar dots already auto-vectorize over taps on hosts +// (integer reduction is exactly reassociable) and measured ~1.5x FASTER +// than the channel-parallel form. Host-only: the embedded targets keep +// their proven planar codegen (Helium on M55, SMLALD on M33-class, +// Hexagon's measured scalar floor — hypotheses C4/C5). +#if !defined(__ARM_FEATURE_MVE) && !defined(__ARM_FEATURE_DSP) && !defined(__hexagon__) +#define TAP_DSP_CHANNEL_PARALLEL 1 +#else +#define TAP_DSP_CHANNEL_PARALLEL 0 +#endif +// Minimum channel count for the frame-major path (overridable for A/B +// measurements; a blend-share planar path stays better at low counts). +#ifndef TAP_DSP_CP_MIN_CHANNELS +#define TAP_DSP_CP_MIN_CHANNELS 4 +#endif + +namespace tap::dsp { + + // ANCHOR: rs_dot_row + /// Dot product of a coefficient row against a history window, in the + /// sample type's accumulator domain (see sample_traits.h): tap-order + /// accumulation, single rounding in finalize. + template + inline S dot_row(const typename sample_traits::coeff* TAP_DSP_RESTRICT row, const S* TAP_DSP_RESTRICT hist, + std::size_t taps) noexcept { + using tr = sample_traits; +#if TAP_DSP_Q15_SMLALD + if constexpr (std::is_same_v) { + std::int64_t acc = 0; + std::size_t t = 0; + for (; t + 1 < taps; t += 2) { + // memcpy keeps the 16-bit pair loads alignment-safe; both + // compile to a single 32-bit load (little-endian packing + // matches SMLALD's lo/hi lanes). + std::uint32_t h; + std::uint32_t r; + std::memcpy(&h, hist + t, sizeof h); + std::memcpy(&r, row + t, sizeof r); + acc = __smlald(static_cast(h), static_cast(r), acc); + } + for (; t < taps; ++t) // odd-tap tail + acc = tr::mac(acc, hist[t], row[t]); + return tr::finalize(acc); + } +#endif + typename tr::accum acc{}; + for (std::size_t t = 0; t < taps; ++t) { + acc = tr::mac(acc, hist[t], row[t]); + } + return tr::finalize(acc); + } + // ANCHOR_END: rs_dot_row + + // ANCHOR: opt_dot_tile + /// One K-channel tile of the channel-parallel dot (hypothesis C6): K + /// accumulators live in a constexpr-size local array — registers, not + /// memory — while the tap loop walks the frame-major window with stride + /// `stride` samples per frame. K is the register-blocking factor; a naive + /// channels-inner loop with accumulators in memory measures ~2.8x SLOWER + /// than planar (each mac round-trips its accumulator through the stack). + template + inline void dot_tile_frame_major(const typename sample_traits::coeff* TAP_DSP_RESTRICT row, + const S* TAP_DSP_RESTRICT x, std::size_t taps, std::size_t stride, + S* TAP_DSP_RESTRICT out) noexcept { + using tr = sample_traits; + typename tr::accum acc[K]{}; + for (std::size_t t = 0; t < taps; ++t) { + const auto coeff = row[t]; + const S* TAP_DSP_RESTRICT frame = x + t * stride; + for (std::size_t k = 0; k < K; ++k) { + acc[k] = tr::mac(acc[k], frame[k], coeff); + } + } + for (std::size_t k = 0; k < K; ++k) { + out[k] = tr::finalize(acc[k]); + } + } + // ANCHOR_END: opt_dot_tile + + // ANCHOR: rs_dot_rows_frame_major + // ANCHOR: opt_dot_rows + /// Channel-parallel dot products over a frame-major history block: all + /// channels' outputs for one frame in register-blocked tiles of 8/4/2/1. + /// Per channel the accumulation order over taps equals dot_row's, so the + /// outputs are bit-exact vs the planar path for every sample type — float + /// included, since each channel's double accumulator still sums the taps + /// in the same order (lanes are channels, not taps). + template + inline void dot_rows_frame_major(const typename sample_traits::coeff* TAP_DSP_RESTRICT row, + const S* TAP_DSP_RESTRICT x, std::size_t taps, std::size_t channels, + S* TAP_DSP_RESTRICT out) noexcept { + std::size_t c = 0; + for (; c + 8 <= channels; c += 8) { + dot_tile_frame_major(row, x + c, taps, channels, out + c); + } + if (c + 4 <= channels) { + dot_tile_frame_major(row, x + c, taps, channels, out + c); + c += 4; + } + if (c + 2 <= channels) { + dot_tile_frame_major(row, x + c, taps, channels, out + c); + c += 2; + } + if (c < channels) { + dot_tile_frame_major(row, x + c, taps, channels, out + c); + } + } + // ANCHOR_END: rs_dot_rows_frame_major + // ANCHOR_END: opt_dot_rows + +} // namespace tap::dsp diff --git a/include/tap/dsp/kaiser.h b/include/tap/dsp/kaiser.h new file mode 100644 index 0000000..99a3fd7 --- /dev/null +++ b/include/tap/dsp/kaiser.h @@ -0,0 +1,375 @@ +/// @file kaiser.h +/// @brief Kaiser-window FIR prototype design for polyphase filter banks. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap (include/srt/detail/kaiser.h), where it designs +// the ASRC's interpolation bank; promoted here so RatioTap's fixed-ratio +// tables and any future FIR design share one implementation. All methods are +// published literature: Kaiser window design (Kaiser 1974), band-limited +// interpolation (J. O. Smith, CCRMA), the harris length estimate, and the +// zeros-at-k*fs droop-compensated variant worked out in the June 2026 +// music-dsp thread (verified in SampleRateTap's asrc_rbj_analysis notebook). +// +// ANCHOR: kai_design_note +/// Design note — runtime vs constexpr: prototype tables run 12K-33K taps and +/// each tap needs sin/sqrt plus a ~50-term Bessel I0 series. Constexpr +/// evaluation is interpreted (roughly 1e3-1e4x slower than native), would need +/// hand-rolled constexpr transcendentals before C++26, and would cost tens of +/// seconds to minutes of compile time in every including translation unit. +/// Runtime design takes well under 10 ms, runs once in a constructor, and is +/// off the audio path, so all design math here is plain runtime double +/// precision. +// ANCHOR_END: kai_design_note +#pragma once + +#include +#include +#include +#include +#include + +namespace tap::dsp { + + // ANCHOR: kai_besseli0 + /// Modified Bessel function of the first kind, order zero, by power series. + /// Converges for all practical Kaiser betas (|x| < ~40); terms are added until + /// they no longer contribute at double precision. + inline double bessel_i0(double x) noexcept { + const double half_x = 0.5 * x; + double term = 1.0; + double sum = 1.0; + for (int k = 1; k < 1000; ++k) { + const double r = half_x / static_cast(k); + term *= r * r; + sum += term; + if (term < 1e-21 * sum) { + break; + } + } + return sum; + } + // ANCHOR_END: kai_besseli0 + + // ANCHOR: kai_beta + /// Kaiser window shape parameter for a given stopband attenuation in dB + /// (Kaiser's published empirical fit). + inline double kaiser_beta(double atten_db) noexcept { + if (atten_db > 50.0) { + return 0.1102 * (atten_db - 8.7); + } + if (atten_db > 21.0) { + return 0.5842 * std::pow(atten_db - 21.0, 0.4) + 0.07886 * (atten_db - 21.0); + } + return 0.0; + } + // ANCHOR_END: kai_beta + + // ANCHOR: kai_estimate + /// Kaiser/harris FIR length estimate, expressed per polyphase branch. + /// + /// \param attenDb target stopband attenuation in dB + /// \param transWidthNorm transition width normalized to the *input* sample rate + /// (e.g. 8 kHz transition at 48 kHz -> 8000/48000) + /// \return estimated taps per polyphase phase: N = (A - 8) / (2.285 * 2*pi * df) + inline std::size_t estimate_taps(double atten_db, double trans_width_norm) noexcept { + // Clamp pathological inputs (attenDb < 8, non-positive width): the raw + // formula goes negative/infinite there and casting that to size_t is UB. + if (!(trans_width_norm > 0.0)) { + return 4; + } + const double n = (atten_db - 8.0) / (2.285 * 2.0 * std::numbers::pi * trans_width_norm); + return n > 4.0 ? static_cast(std::ceil(n)) : 4; + } + // ANCHOR_END: kai_estimate + + // ANCHOR: kai_sinc + /// sin(pi x)/(pi x) with the removable singularity handled. + inline double sinc(double x) noexcept { + if (std::abs(x) < 1e-12) { + return 1.0; + } + const double px = std::numbers::pi * x; + return std::sin(px) / px; + } + // ANCHOR_END: kai_sinc + + // ANCHOR: kai_prototype + /// Designs the Kaiser-windowed sinc prototype lowpass for an L-phase + /// interpolation bank. + /// + /// \param h output, length L*T (the full oversampled prototype) + /// \param numPhases L; the prototype is sampled on a grid of 1/L input samples + /// \param cutoffNorm cutoff normalized to the input Nyquist, i.e. 2*fc/fs in + /// (0, 1]; for a near-unity interpolator centered between a + /// 20 kHz passband and 28 kHz stopband at 48 kHz this is + /// (20000+28000)/48000 = 1.0 (cutoff at input Nyquist) + /// \param beta Kaiser shape parameter (see kaiserBeta) + /// + /// The result is normalized so that sum(h) == L, giving each polyphase branch a + /// DC gain of ~1 (deviation bounded by the stopband leakage). + inline void design_prototype(std::span h, std::size_t num_phases, double cutoff_norm, + double beta) noexcept { + const std::size_t n = h.size(); + const double center = 0.5 * static_cast(n - 1); + const double i0_beta = bessel_i0(beta); + double sum = 0.0; + for (std::size_t i = 0; i < n; ++i) { + const double t = (static_cast(i) - center) / static_cast(num_phases); + const double u = (static_cast(i) - center) / center; // window argument, [-1, 1] + const double w = bessel_i0(beta * std::sqrt(std::max(0.0, 1.0 - u * u))) / i0_beta; + h[i] = cutoff_norm * sinc(cutoff_norm * t) * w; + sum += h[i]; + } + const double gain = static_cast(num_phases) / sum; + for (auto& v : h) { + v *= gain; + } + } + // ANCHOR_END: kai_prototype + + /// Solves the dense n x n system m * out = rhs in place (Gaussian elimination + /// with partial pivoting; row-major m). Small systems only — the compensated + /// design below solves at most 15 unknowns, and the analysis instruments' + /// joint tone fits stay under ~50. + inline void solve_dense(std::span m, std::span rhs, std::span out, std::size_t n) noexcept { + std::vector order(n); + for (std::size_t i = 0; i < n; ++i) { + order[i] = i; + } + for (std::size_t col = 0; col < n; ++col) { + std::size_t piv = col; + for (std::size_t r = col + 1; r < n; ++r) { + if (std::abs(m[order[r] * n + col]) > std::abs(m[order[piv] * n + col])) { + piv = r; + } + } + std::swap(order[col], order[piv]); + const std::size_t p = order[col]; + for (std::size_t r = col + 1; r < n; ++r) { + const std::size_t rr = order[r]; + const double f = m[rr * n + col] / m[p * n + col]; + for (std::size_t q = col; q < n; ++q) { + m[rr * n + q] -= f * m[p * n + q]; + } + rhs[rr] -= f * rhs[p]; + } + } + for (std::size_t col = n; col-- > 0;) { + const std::size_t p = order[col]; + double v = rhs[p]; + for (std::size_t q = col + 1; q < n; ++q) { + v -= m[p * n + q] * out[q]; + } + out[col] = v / m[p * n + col]; + } + } + + // ANCHOR: pw_comp_design + /// Designs a prototype with transmission zeros at every integer multiple of + /// the sample rate, passband droop pre-compensated (the music-dsp thread's + /// suggestion, done inside the passband spec — see SampleRateTap's epilogue + /// chapter and its asrc_rbj_analysis notebook). + /// + /// Construction: the zeros come from convolving with a one-input-sample rect + /// (multiplies the response by sinc(f/fs), zero at every k*fs — exactly where + /// the images of low-frequency program energy sit). The rect's passband droop + /// (-2.64 dB at 20 kHz) is cancelled by tilting the design target by + /// 1/sinc(f/fs), expressed as a short cosine series in f — which in time is a + /// weighted sum of the brickwall kernel at small integer shifts, so the whole + /// design stays closed-form: no FFT, no dependency. One correction pass + /// (measure the built passband by direct DFT, fold the deviation back into + /// the tilt) holds every preset's ripple within +/-0.005 dB, a >=2x margin + /// on the spec. (A second pass buys ~1 dB of margin for another kernel + /// build plus probe sweep — real money on soft-double targets, where every + /// design flop is a libcall; the pass count and probe count below are + /// sized by the M33 instruction-count ledger in SampleRateTap's + /// docs/PERFORMANCE.md.) + /// + /// \param h output, length L*T for the TOTAL taps per phase T; the + /// sinc design uses T-1 taps and the rect supplies the + /// +1 (composite length L*(T-1)+L-1, one zero of padding) + /// \param numPhases L + /// \param cutoffNorm as designPrototype + /// \param beta Kaiser shape parameter for the T-1-tap base design + /// \param passbandNorm passband edge / sample rate (flatness is corrected and + /// verified up to here) + /// + /// Costs a few ms more than designPrototype (three kernel builds plus ~100 + /// direct-DFT probes); still constructor-only, off the audio path. Allocates + /// workspace; may throw std::bad_alloc. + // ANCHOR_END: pw_comp_design + inline void design_prototype_compensated(std::span h, std::size_t num_phases, double cutoff_norm, + double beta, double passband_norm) { + const std::size_t L = num_phases; + const std::size_t total = h.size() / L; // total taps per phase (with rect) + const std::size_t td = total - 1; // sinc-design taps per phase + const std::size_t n = L * td; // fine-grid design length + const std::size_t nc = n + L - 1; // composite length after rect + // Compensator order: enough cosine terms to hold the passband tilt to + // ~1e-4, capped so the shifted kernels stay well inside short windows. + const std::size_t M = std::min(14, (td - 1) / 5); + + constexpr std::size_t k_grid = 1001; // fit grid over f/fs in [0, 0.5] + constexpr std::size_t k_probe = 24; // passband correction probes + std::vector target(k_grid), a(M + 1), fine(n), probe(k_probe); + for (std::size_t g = 0; g < k_grid; ++g) { + const double f = 0.5 * static_cast(g) / static_cast(k_grid - 1); + const double pf = std::numbers::pi * f; + target[g] = f < 1e-9 ? 1.0 : pf / std::sin(pf); // 1/sinc(f/fs) + } + + const auto fit_cosine_series = [&] { + // Weighted LS of target on cos(2*pi*m*f): exact where flatness is + // specified (the passband, heavy weight), merely tracked above it. + // Basis by Chebyshev recurrence: cos(m x) from the two previous + // orders, one real cosine per grid point. + std::vector nm((M + 1) * (M + 1), 0.0), rhs(M + 1, 0.0), basis(M + 1); + for (std::size_t g = 0; g < k_grid; ++g) { + const double f = 0.5 * static_cast(g) / static_cast(k_grid - 1); + const double w2 = f <= passband_norm + 0.02 ? 1e8 : 1.0; // (weight 1e4)^2 + const double c1 = std::cos(2.0 * std::numbers::pi * f); + basis[0] = 1.0; + if (M >= 1) { + basis[1] = c1; + } + for (std::size_t m = 2; m <= M; ++m) { + basis[m] = 2.0 * c1 * basis[m - 1] - basis[m - 2]; + } + for (std::size_t r = 0; r <= M; ++r) { + for (std::size_t q = 0; q <= M; ++q) { + nm[r * (M + 1) + q] += w2 * basis[r] * basis[q]; + } + rhs[r] += w2 * basis[r] * target[g]; + } + } + solve_dense(nm, rhs, a, M + 1); + }; + + const auto build = [&] { + // Tilted ideal kernel: sum of the brickwall sinc at integer shifts. + // Centered at n/2 (not (n-1)/2): the even-length rect below shifts + // the composite by (L-1)/2, and n/2 + (L-1)/2 == (L*total - 1)/2 — + // the exact center designPrototype uses, so a consuming bank's + // phase/delay convention is identical for both designs. (Getting + // this wrong is a half-fine-sample delay error: ~-72 dB at 1 kHz, + // worse by 6 dB per octave — SampleRateTap's fractional-delay + // accuracy tests catch it.) + // + // Transcendental budget: the naive form of this loop calls libm sin + // once per tap per compensator term (~2M calls across the design; a + // measured +225M constructor instructions on Cortex-M55, and worse + // where doubles are soft). Instead: sin(pi*c*(t -+ m)) expands by + // angle addition over precomputed sin/cos(pi*c*m), and sin/cos of + // the per-tap angle advance by a unit rotator, re-synced with real + // libm calls every 4096 taps to bound drift far below the design's + // own accuracy floor. + const double center = 0.5 * static_cast(n); + const double i0_beta = bessel_i0(beta); + std::vector cs(M + 1), sn(M + 1); + for (std::size_t m = 0; m <= M; ++m) { + cs[m] = std::cos(std::numbers::pi * cutoff_norm * static_cast(m)); + sn[m] = std::sin(std::numbers::pi * cutoff_norm * static_cast(m)); + } + const double step = std::numbers::pi * cutoff_norm / static_cast(L); + const double step_c = std::cos(step), step_s = std::sin(step); + double ang_s = 0.0, ang_c = 1.0; // sin/cos(pi*c*t_i), re-synced below + for (std::size_t i = 0; i < n; ++i) { + const double t = (static_cast(i) - center) / static_cast(L); + if (i % 4096 == 0) { + ang_s = std::sin(std::numbers::pi * cutoff_norm * t); + ang_c = std::cos(std::numbers::pi * cutoff_norm * t); + } + const auto shifted_sinc = [&](double dm, double sin_shift, double cos_shift) { + const double x = cutoff_norm * (t - dm); // dm may be negative + if (std::abs(x) < 1e-12) { + return 1.0; + } + // sin(pi*c*(t - dm)) = sin(pi*c*t)cos(pi*c*dm) - cos(..)sin(..) + return (ang_s * cos_shift - ang_c * sin_shift) / (std::numbers::pi * x); + }; + double v = a[0] * cutoff_norm * shifted_sinc(0.0, 0.0, 1.0); + for (std::size_t m = 1; m <= M; ++m) { + const double dm = static_cast(m); + v += 0.5 * a[m] * cutoff_norm * (shifted_sinc(dm, sn[m], cs[m]) + shifted_sinc(-dm, -sn[m], cs[m])); + } + const double u = (static_cast(i) - center) / center; + fine[i] = v * bessel_i0(beta * std::sqrt(std::max(0.0, 1.0 - u * u))) / i0_beta; + const double next_s = ang_s * step_c + ang_c * step_s; + ang_c = ang_c * step_c - ang_s * step_s; + ang_s = next_s; + } + // ANCHOR: pw_comp_rect + // Rect convolution as a running sum: exact zeros at every k*fs. + double run = 0.0; + for (std::size_t i = 0; i < nc; ++i) { + run += i < n ? fine[i] : 0.0; + if (i >= L) { + run -= fine[i - L]; + } + h[i] = run / static_cast(L); + } + for (std::size_t i = nc; i < h.size(); ++i) { + h[i] = 0.0; + } + // ANCHOR_END: pw_comp_rect + double sum = 0.0; + for (std::size_t i = 0; i < nc; ++i) { + sum += h[i]; + } + const double gain = static_cast(L) / sum; + for (std::size_t i = 0; i < nc; ++i) { + h[i] *= gain; + } + }; + + for (int pass = 0; pass < 1; ++pass) { + fit_cosine_series(); + build(); + // Probe the built passband by direct DFT (cos projection about the + // composite's symmetry center, (L*total - 1)/2 == nc/2) and fold the + // deviation into the tilt. One rotator per probe frequency: two libm + // calls each instead of one per tap (rotator drift over ~2^14 steps + // is ~1e-12, five orders below the ripple being measured). + const double center = 0.5 * static_cast(nc); + for (std::size_t j = 0; j < k_probe; ++j) { + const double f = passband_norm * static_cast(j + 1) / k_probe; + const double th = 2.0 * std::numbers::pi * f / static_cast(L); + const double th_c = std::cos(th), th_s = std::sin(th); + double rc = std::cos(th * -center), rs = std::sin(th * -center); + double acc = 0.0; + for (std::size_t i = 0; i < nc; ++i) { + acc += h[i] * rc; + const double nrc = rc * th_c - rs * th_s; + rs = rs * th_c + rc * th_s; + rc = nrc; + } + probe[j] = std::abs(acc) / static_cast(L); + } + for (std::size_t g = 0; g < k_grid; ++g) { + const double f = 0.5 * static_cast(g) / static_cast(k_grid - 1); + if (f > passband_norm) { + continue; + } + // probe[j] sits at f = passbandNorm*(j+1)/kProbe, i.e. x = j+1 + const double x = f / passband_norm * k_probe - 1.0; + double d; + if (x <= 0.0) { + d = probe[0]; + } + else if (x >= static_cast(k_probe - 1)) { + d = probe[k_probe - 1]; + } + else { + const auto j = static_cast(x); + const double fr = x - static_cast(j); + d = probe[j] * (1.0 - fr) + probe[j + 1] * fr; + } + target[g] /= std::max(d, 0.5); + } + } + fit_cosine_series(); + build(); + } + +} // namespace tap::dsp diff --git a/include/tap/dsp/quantize.h b/include/tap/dsp/quantize.h new file mode 100644 index 0000000..ff46f28 --- /dev/null +++ b/include/tap/dsp/quantize.h @@ -0,0 +1,83 @@ +/// @file quantize.h +/// @brief Row-sum-preserving coefficient quantization for polyphase tables. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap, where this correction ran inline in the +// polyphase bank constructor; promoted here because every polyphase table +// with unity per-phase DC gain wants it — RatioTap's fixed-ratio tables +// exactly as much as the ASRC's interpolation bank. +#pragma once + +#include +#include +#include +#include +#include + +#include "tap/dsp/sample_traits.h" + +namespace tap::dsp { + + // ANCHOR: pw_row_sum + /// Quantizes one polyphase branch's coefficients (double, in units where + /// 1.0 is unity gain) to the sample type's coefficient format, preserving + /// the row's DC sum exactly. + /// + /// "Note that for DC, this should get you infinite S/N ratio... for every + /// phase or fractional delay, the FIR coefficients must add to 1." + /// -- R. Bristow-Johnson, music-dsp. In double, a well-designed bank + /// makes every branch's DC sum identical to machine epsilon (zeros at + /// k*fs ARE branch-DC uniformity, stated in frequency), but independent + /// per-tap rounding re-breaks it by several LSB. This distributes each + /// row's total rounding residual to the taps that were rounded furthest + /// from it (largest-remainder method): every row then sums to + /// llround(exact_sum * k_coeff_scale), so a table built row-by-row holds + /// DC gain within one coefficient LSB across all phases. + /// + /// For floating-point coefficient types this is a plain conversion (no + /// correction runs; k_coeff_scale is 1). Design-time code: allocates a + /// scratch vector for integer formats, so keep it off the audio path. + /// + /// \pre dst.size() == src.size() + template + inline void quantize_row_preserving_sum(std::span src, + std::span::coeff> dst) { + using tr = sample_traits; + using coeff = typename tr::coeff; + const std::size_t n = src.size(); + if constexpr (std::is_floating_point_v) { + for (std::size_t t = 0; t < n; ++t) { + dst[t] = tr::make_coeff(src[t]); + } + } + else { + std::vector remainder(n); + double exact_sum = 0.0; + std::int64_t quant_sum = 0; + for (std::size_t t = 0; t < n; ++t) { + const double scaled = src[t] * tr::k_coeff_scale; + const coeff q = tr::make_coeff(src[t]); + dst[t] = q; + remainder[t] = scaled - static_cast(q); + exact_sum += scaled; + quant_sum += static_cast(q); + } + std::int64_t residual = static_cast(std::llround(exact_sum)) - quant_sum; + while (residual != 0) { + const double sgn = residual > 0 ? 1.0 : -1.0; + std::size_t best = 0; + for (std::size_t u = 1; u < n; ++u) { + if (sgn * remainder[u] > sgn * remainder[best]) { + best = u; + } + } + dst[best] = static_cast(dst[best] + (residual > 0 ? 1 : -1)); + remainder[best] -= sgn; + residual -= residual > 0 ? 1 : -1; + } + } + } + // ANCHOR_END: pw_row_sum + +} // namespace tap::dsp diff --git a/include/tap/dsp/sample_traits.h b/include/tap/dsp/sample_traits.h new file mode 100644 index 0000000..7ad396e --- /dev/null +++ b/include/tap/dsp/sample_traits.h @@ -0,0 +1,202 @@ +/// @file sample_traits.h +/// @brief Sample-format customization point for FIR datapaths: float, Q15, Q31. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap (include/srt/sample_traits.h), where it +// parameterizes the ASRC datapath; promoted here as the family's sample-format +// substrate, shared with RatioTap's fixed-ratio hot loop and available to any +// primitive growing an embedded fixed-point profile. This header is the +// *format core* only — how each sample type stores coefficients, accumulates +// dot products, and rounds/saturates back to samples. Engine-specific +// extensions (e.g. SampleRateTap's inter-phase coefficient blending) layer on +// top by deriving from these specializations and refining the concept. +// +// Three sample types are provided: +// +// - float : float I/O and coefficients, double accumulation +// - std::int16_t : Q15 samples, Q1.14 coefficients, int64 accumulation, +// saturating output +// - std::int32_t : Q31 samples, Q1.30 coefficients, int64 accumulation, +// saturating output +// +// The fixed-point profiles exist for targets where double (and sometimes any +// float) is unaffordable: M33-class cores measured ~19x the instruction count +// running SampleRateTap's float datapath vs Q15 (soft-double accumulation). +// Expected deployments include Bluetooth-adjacent converters (RatioTap) and +// M33/M55-class eurorack/pedal targets (TapTools primitives, opt-in per +// primitive — each adoption is its own documented Q-format design). +// +// Deliberately raw sample types + traits, not wrapper classes with operator +// overloads: the Q-format ladder (headroom bits, accumulator pre-shifts, +// where the single rounding happens) is a per-datapath design decision that +// must stay visible at the use site and pinnable by tests — see the Q31 +// mac() below for an example a `q31` class could not express honestly. Raw +// int16_t/int32_t also match what codecs and C ABIs actually deliver, and +// keep memcpy/SIMD access patterns (e.g. the SMLALD pair loads in +// fir_kernels.h) legal. +#pragma once + +#include +#include +#include +#include + +namespace tap::dsp { + + namespace detail { + + // ANCHOR: st_roundsat + /// Round-and-saturate a double to a signed integer coefficient/sample type. + template + constexpr I round_sat(double v) noexcept { + constexpr double lo = static_cast(std::numeric_limits::min()); + constexpr double hi = static_cast(std::numeric_limits::max()); + const double r = v < 0.0 ? v - 0.5 : v + 0.5; // round half away from zero + if (r <= lo) { + return std::numeric_limits::min(); + } + if (r >= hi) { + return std::numeric_limits::max(); + } + return static_cast(r); + } + // ANCHOR_END: st_roundsat + + /// Saturate a 64-bit accumulator result to a narrower signed integer. + template + constexpr I clamp_sat(std::int64_t v) noexcept { + constexpr auto lo = static_cast(std::numeric_limits::min()); + constexpr auto hi = static_cast(std::numeric_limits::max()); + return static_cast(v < lo ? lo : (v > hi ? hi : v)); + } + + } // namespace detail + + /// Primary template intentionally undefined; specialize per sample type. + template + struct sample_traits; + + /// Float datapath: float samples and coefficients, double accumulation. + /// The double accumulator keeps the dot-product noise floor far below a + /// 120 dB transparency target; float coefficient storage quantizes the + /// filter at roughly -150 dB, negligible against the same target. + template <> + struct sample_traits { + using coeff = float; ///< stored filter coefficient type + using accum = double; ///< dot-product accumulator type + + /// Convert a double-precision designed coefficient to storage form. + static coeff make_coeff(double c) noexcept { return static_cast(c); } + /// Coefficient units per 1.0 (used by row-sum-preserving quantization; + /// unity for floating storage, where no correction runs). + static constexpr double k_coeff_scale = 1.0; + + /// acc + x * c, in the accumulator domain. + static accum mac(accum acc, float x, coeff c) noexcept { + return acc + static_cast(x) * static_cast(c); + } + + /// Convert the accumulator to an output sample (saturates for fixed point). + static float finalize(accum acc) noexcept { return static_cast(acc); } + + /// The zero/silence sample value. + static float silence() noexcept { return 0.0f; } + }; + + // ANCHOR: st_q15_core + /// Q15 fixed-point datapath (samples are int16_t in Q0.15). + /// + /// Coefficients are stored in Q1.14: a unity-DC prototype's peak tap + /// reaches ~1.0, which does not fit Q0.15, so one headroom bit is traded + /// for one precision bit. Products are Q0.15 x Q1.14 = Q29 and are summed + /// exactly in int64 (48-80 taps add ~6-7 bits — no overflow, no + /// intermediate rounding). The single rounding happens in finalize(): + /// Q29 -> Q15 with round-half-up and saturation. Coefficient quantization + /// (Q14, ~-86 dB) and output quantization (Q15) set the noise floor — both + /// at the format's own limit, so a converter built on this is + /// Q15-transparent. + template <> + struct sample_traits { + using coeff = std::int16_t; + using accum = std::int64_t; + + // ANCHOR: st_q15_coeff + static coeff make_coeff(double c) noexcept { + return detail::round_sat(c * 16384.0); // Q1.14 + } + static constexpr double k_coeff_scale = 16384.0; // Q1.14 units per 1.0 + // ANCHOR_END: st_q15_coeff + + // ANCHOR: st_q15_mac + static accum mac(accum acc, std::int16_t x, coeff c) noexcept { + return acc + static_cast(static_cast(x) * static_cast(c)); + } + // ANCHOR_END: st_q15_mac + + // ANCHOR: st_q15_finalize + static std::int16_t finalize(accum acc) noexcept { + // Round-half-up, not half-even: the bias is a fraction of one + // sub-LSB rounding step, far below the Q15 noise floor. + return detail::clamp_sat((acc + (1 << 13)) >> 14); // Q29 -> Q15 + } + // ANCHOR_END: st_q15_finalize + + static std::int16_t silence() noexcept { return 0; } + }; + // ANCHOR_END: st_q15_core + + // ANCHOR: st_q31_core + /// Q31 fixed-point datapath (samples are int32_t in Q0.31). + /// + /// Coefficients are stored in Q1.30 (one headroom bit for the ~1.0 peak + /// tap). A full-precision product would be Q0.31 x Q1.30 = 62 bits, which + /// overflows int64 once ~48 of them are summed, so each product is + /// pre-shifted down 16 bits (Q45) before accumulation; the discarded bits + /// sit 14 bits below the final Q31 LSB, far beneath the format's noise + /// floor. finalize() rounds Q45 -> Q31 with saturation. + template <> + struct sample_traits { + using coeff = std::int32_t; + using accum = std::int64_t; + + static coeff make_coeff(double c) noexcept { + return detail::round_sat(c * 1073741824.0); // Q1.30 + } + static constexpr double k_coeff_scale = 1073741824.0; // Q1.30 units per 1.0 + + // ANCHOR: st_q31_mac + static accum mac(accum acc, std::int32_t x, coeff c) noexcept { + return acc + ((static_cast(x) * c) >> 16); // Q61 -> Q45 + } + // ANCHOR_END: st_q31_mac + + static std::int32_t finalize(accum acc) noexcept { + return detail::clamp_sat((acc + (1 << 13)) >> 14); // Q45 -> Q31 + } + + static std::int32_t silence() noexcept { return 0; } + }; + // ANCHOR_END: st_q31_core + + // ANCHOR: st_core_concept + /// Satisfied by any type with a complete format-core sample_traits + /// specialization — everything the FIR kernels (fir_kernels.h) and table + /// builders require. Consumers with richer datapaths refine this concept + /// over their own trait extensions (SampleRateTap's sample_type adds the + /// inter-phase blend contract on top). + template + concept sample_type = + requires(T x, double d, typename sample_traits::accum a, typename sample_traits::coeff c) { + { sample_traits::make_coeff(d) } -> std::same_as::coeff>; + { sample_traits::mac(a, x, c) } -> std::same_as::accum>; + { sample_traits::finalize(a) } -> std::same_as; + { sample_traits::silence() } -> std::same_as; + }; + + static_assert(sample_type); + static_assert(sample_type); + static_assert(sample_type); + // ANCHOR_END: st_core_concept + +} // namespace tap::dsp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 75f134c..e3abd01 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -20,10 +20,15 @@ set(INSTALL_GTEST OFF CACHE BOOL "" FORCE) FetchContent_MakeAvailable(googletest) add_executable(tap_dsp_tests + test_analysis.cpp test_fft.cpp test_fft_backend.cpp + test_fir_kernels.cpp + test_kaiser.cpp test_psola.cpp test_pvoc.cpp + test_quantize.cpp + test_sample_traits.cpp test_yin.cpp) target_link_libraries(tap_dsp_tests PRIVATE tap::dsp diff --git a/tests/test_analysis.cpp b/tests/test_analysis.cpp new file mode 100644 index 0000000..235e1e8 --- /dev/null +++ b/tests/test_analysis.cpp @@ -0,0 +1,79 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Smoke battery for the measurement instruments. These headers earn their +// keep in consumer suites (SampleRateTap's quality tests, RatioTap's alias +// acceptance); here we pin that the instruments themselves reach their +// documented floors on exact synthetic signals — an instrument that cannot +// measure a clean signal cleanly would silently weaken every consumer's +// quality gate. + +#include +#include +#include + +#include + +#include "tap/dsp/analysis/multitone_analysis.h" +#include "tap/dsp/analysis/sine_analysis.h" + +namespace { + + namespace an = tap::dsp::analysis; + + TEST(SineAnalysis, FitRecoversExactTone) { + const double nu = 997.0 / 48000.0; + std::vector x(16384); + for (std::size_t i = 0; i < x.size(); ++i) { + x[i] = static_cast(0.5 * std::sin(2.0 * std::numbers::pi * nu * static_cast(i) + 0.3)); + } + const an::sine_fit fit = an::fit_sine(x, nu); + EXPECT_NEAR(fit.amplitude, 0.5, 1e-6); + EXPECT_NEAR(fit.dc, 0.0, 1e-6); + // Residual is the float storage quantization only: > 100 dB SNR. + EXPECT_GT(an::snr_db(fit), 100.0); + } + + TEST(SineAnalysis, TrackedFitAbsorbsSmallFrequencyOffset) { + // The tone sits 5 ppm off the guess — the situation the tracker exists + // for (a converter's rate estimate settling asymptotically). A rigid + // fit at the guess books the offset as residual; the tracked fit must + // recover the true frequency and the quantization-limited floor. + const double nu_true = (997.0 / 48000.0) * (1.0 + 5e-6); + std::vector x(65536); + for (std::size_t i = 0; i < x.size(); ++i) { + x[i] = static_cast(0.5 * std::sin(2.0 * std::numbers::pi * nu_true * static_cast(i))); + } + const an::sine_fit fit = an::fit_sine_tracked(x, 997.0 / 48000.0); + EXPECT_NEAR(fit.freq_norm, nu_true, 1e-10); + EXPECT_NEAR(fit.amplitude, 0.5, 1e-6); + EXPECT_GT(an::snr_db(fit), 100.0); + } + + TEST(MultitoneAnalysis, CombIsBoundedByPeakSum) { + const an::tone_comb comb = an::tone_comb::pink(24, 40.0, 20000.0, 0.5); + ASSERT_EQ(comb.freq_hz.size(), 24u); + double amp_sum = 0.0; + for (const double a : comb.amplitude) { + amp_sum += a; + } + EXPECT_NEAR(amp_sum, 0.5, 1e-12); + for (std::uint64_t i = 0; i < 48000; ++i) { + ASSERT_LE(std::abs(comb.sample_at(i, 48000.0)), 0.5 + 1e-9); + } + } + + TEST(MultitoneAnalysis, JointFitReachesQuantizationFloorOnExactTones) { + // The documented instrument floor claim: sequential fit-subtract + // floors near 48 dB on exact synthetic tones; the joint solve reaches + // the float storage quantization floor. Pin the latter. + const an::tone_comb comb = an::tone_comb::pink(12, 40.0, 20000.0, 0.5); + std::vector tail(32768); + for (std::size_t i = 0; i < tail.size(); ++i) { + tail[i] = static_cast(comb.sample_at(i, 48000.0)); + } + const double snr = an::program_weighted_snr_db(tail, comb, 48000.0, 48000.0); + EXPECT_GT(snr, 90.0); + } + +} // namespace diff --git a/tests/test_fir_kernels.cpp b/tests/test_fir_kernels.cpp new file mode 100644 index 0000000..cfc5cc5 --- /dev/null +++ b/tests/test_fir_kernels.cpp @@ -0,0 +1,130 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the FIR dot kernels. The load-bearing promise — +// extracted from SampleRateTap's multichannel suite, where it gates the +// frame-major fast path — is bit-exactness: the channel-parallel kernel must +// produce the identical bits as the planar dot for every sample type, because +// consumers switch between the layouts by channel count and target. + +#include +#include + +#include + +#include "tap/dsp/fir_kernels.h" + +namespace { + + using tap::dsp::dot_row; + using tap::dsp::dot_rows_frame_major; + using tap::dsp::sample_traits; + + // Deterministic pseudo-random generator (xorshift), mapped per sample type + // to well-inside-full-scale values so no test depends on saturation. + inline std::uint32_t next(std::uint32_t& s) { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + return s; + } + + template + S sample_from(std::uint32_t r); + template <> + float sample_from(std::uint32_t r) { + return (static_cast(r % 65536) - 32768.0f) / 65536.0f; + } + template <> + std::int16_t sample_from(std::uint32_t r) { + return static_cast(static_cast(r % 32768) - 16384); + } + template <> + std::int32_t sample_from(std::uint32_t r) { + return static_cast(r) / 4; + } + + template + typename sample_traits::coeff coeff_from(std::uint32_t r) { + // Coefficients in roughly [-0.5, 0.5) of the format's unity. + return sample_traits::make_coeff((static_cast(r % 4096) - 2048.0) / 4096.0); + } + + template + class fir_kernels_test : public ::testing::Test {}; + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(fir_kernels_test, sample_types, ); + + // dot_row must equal the reference accumulation: mac per tap in order, + // one finalize. (On SMLALD targets this also pins the dual-MAC pairing; + // pairing is bit-free because each 16x16 product is exact in int32.) + TYPED_TEST(fir_kernels_test, DotRowMatchesReferenceAccumulation) { + using sample = TypeParam; + using tr = sample_traits; + constexpr std::size_t k_taps = 48; + std::uint32_t seed = 0x12345678u; + std::vector hist(k_taps); + std::vector row(k_taps); + for (std::size_t t = 0; t < k_taps; ++t) { + hist[t] = sample_from(next(seed)); + row[t] = coeff_from(next(seed)); + } + typename tr::accum acc{}; + for (std::size_t t = 0; t < k_taps; ++t) { + acc = tr::mac(acc, hist[t], row[t]); + } + EXPECT_EQ(dot_row(row.data(), hist.data(), k_taps), tr::finalize(acc)); + } + + // The frame-major channel-parallel kernel is bit-exact against the planar + // dot for every channel count that exercises the 8/4/2/1 tiling, float + // included (lanes are channels, not taps: each channel's double + // accumulation order is unchanged). + TYPED_TEST(fir_kernels_test, ChannelParallelMatchesPlanarBitExact) { + using sample = TypeParam; + using tr = sample_traits; + constexpr std::size_t k_taps = 48; + for (std::size_t channels : {1u, 2u, 3u, 4u, 7u, 8u, 11u, 12u, 16u}) { + std::uint32_t seed = 0x9e3779b9u + static_cast(channels); + std::vector x(k_taps * channels); // frame-major + std::vector row(k_taps); + for (auto& v : x) { + v = sample_from(next(seed)); + } + for (auto& c : row) { + c = coeff_from(next(seed)); + } + std::vector out(channels); + dot_rows_frame_major(row.data(), x.data(), k_taps, channels, out.data()); + for (std::size_t c = 0; c < channels; ++c) { + std::vector planar(k_taps); + for (std::size_t t = 0; t < k_taps; ++t) { + planar[t] = x[t * channels + c]; + } + EXPECT_EQ(out[c], dot_row(row.data(), planar.data(), k_taps)) + << "channels=" << channels << " c=" << c; + } + } + } + + // Odd tap counts exercise dot_row's scalar tail on SMLALD targets; on + // hosts this just pins the same contract at a second geometry. + TYPED_TEST(fir_kernels_test, OddTapCountMatchesReference) { + using sample = TypeParam; + using tr = sample_traits; + constexpr std::size_t k_taps = 33; + std::uint32_t seed = 0xdeadbeefu; + std::vector hist(k_taps); + std::vector row(k_taps); + for (std::size_t t = 0; t < k_taps; ++t) { + hist[t] = sample_from(next(seed)); + row[t] = coeff_from(next(seed)); + } + typename tr::accum acc{}; + for (std::size_t t = 0; t < k_taps; ++t) { + acc = tr::mac(acc, hist[t], row[t]); + } + EXPECT_EQ(dot_row(row.data(), hist.data(), k_taps), tr::finalize(acc)); + } + +} // namespace diff --git a/tests/test_kaiser.cpp b/tests/test_kaiser.cpp new file mode 100644 index 0000000..813cd20 --- /dev/null +++ b/tests/test_kaiser.cpp @@ -0,0 +1,196 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the Kaiser prototype designer, ported from +// SampleRateTap's test_kaiser.cpp. The spec structs below carry that +// library's shipping preset numbers (fast/balanced/transparent/economy) so +// the designs it depends on stay pinned at their source. + +#include +#include +#include +#include +#include + +#include + +#include "tap/dsp/kaiser.h" + +namespace { + + using namespace tap::dsp; + + TEST(Kaiser, BesselI0ReferenceValues) { + EXPECT_DOUBLE_EQ(bessel_i0(0.0), 1.0); + EXPECT_NEAR(bessel_i0(1.0), 1.2660658777520084, 1e-12); + EXPECT_NEAR(bessel_i0(5.0), 27.239871823604442, 1e-9); + EXPECT_NEAR(bessel_i0(12.0), 18948.925349296309, 1e-6); + } + + TEST(Kaiser, BetaReferenceValues) { + EXPECT_NEAR(kaiser_beta(120.0), 0.1102 * (120.0 - 8.7), 1e-12); + EXPECT_NEAR(kaiser_beta(40.0), 0.5842 * std::pow(19.0, 0.4) + 0.07886 * 19.0, 1e-12); + EXPECT_DOUBLE_EQ(kaiser_beta(15.0), 0.0); + } + + TEST(Kaiser, TapEstimateMatchesHarrisFormula) { + // 120 dB over a 20->28 kHz transition at 48 kHz: ~47 taps per phase. + const std::size_t taps = estimate_taps(120.0, 8000.0 / 48000.0); + EXPECT_GE(taps, 45u); + EXPECT_LE(taps, 49u); + } + + // A prototype specification in SampleRateTap's filter_spec vocabulary. + struct proto_spec { + std::size_t num_phases; + std::size_t taps_per_phase; + double passband_hz; + double stopband_hz; + double stopband_atten_db; + bool image_zeros; + + proto_spec scaled_to(double sample_rate_hz) const { + constexpr double k_design_rate_hz = 48000.0; + proto_spec s = *this; + s.passband_hz *= sample_rate_hz / k_design_rate_hz; + s.stopband_hz *= sample_rate_hz / k_design_rate_hz; + return s; + } + }; + + // SampleRateTap's shipping presets (filter_spec::fast/balanced/transparent/ + // economy), pinned here at the substrate so a design change is caught where + // the code now lives. + constexpr proto_spec k_fast{128, 32, 18000.0, 30000.0, 96.0, false}; + constexpr proto_spec k_balanced{256, 48, 20000.0, 28000.0, 120.0, true}; + constexpr proto_spec k_transparent{512, 80, 20000.0, 26000.0, 140.0, true}; + constexpr proto_spec k_economy{512, 32, 18000.0, 30000.0, 96.0, true}; + + // Direct DFT magnitude of the double-precision prototype, normalized so the + // passband sits at 0 dB. f is in Hz; the prototype rate is L * fs. + double response_db(const std::vector& h, std::size_t num_phases, double fs, double f) { + const double proto_rate = static_cast(num_phases) * fs; + std::complex acc{0.0, 0.0}; + for (std::size_t m = 0; m < h.size(); ++m) { + const double ang = -2.0 * std::numbers::pi * f * static_cast(m) / proto_rate; + acc += h[m] * std::polar(1.0, ang); + } + return 20.0 * std::log10(std::abs(acc) / static_cast(num_phases)); + } + + void check_prototype_meets_spec(const proto_spec& spec, double fs) { + const std::size_t phases = std::bit_ceil(spec.num_phases); + const std::size_t n = phases * spec.taps_per_phase; + std::vector h(n); + const double cutoff_norm = (spec.passband_hz + spec.stopband_hz) / fs; + if (spec.image_zeros) { + design_prototype_compensated(h, phases, cutoff_norm, kaiser_beta(spec.stopband_atten_db), + spec.passband_hz / fs); + } + else { + design_prototype(h, phases, cutoff_norm, kaiser_beta(spec.stopband_atten_db)); + } + + // Passband: flat within +/-0.01 dB up to the passband edge. For the + // compensated designs this is the claim the droop pre-compensation + // exists to defend (the raw rect would sag -2.64 dB at 20 kHz). + for (double f = 0.0; f <= spec.passband_hz; f += 500.0) { + EXPECT_NEAR(response_db(h, spec.num_phases, fs, f), 0.0, 0.01) << "passband deviation at " << f << " Hz"; + } + + // Stopband: at least the rated attenuation (1 dB grace) from the stopband + // edge out to well past the first few images. + for (double f = spec.stopband_hz; f <= 3.0 * fs; f += 250.0) { + EXPECT_LT(response_db(h, spec.num_phases, fs, f), -(spec.stopband_atten_db - 1.0)) + << "stopband leakage at " << f << " Hz"; + } + + // Transmission zeros at every k*fs: exact in exact arithmetic, so demand + // far below the rated stopband (double rounding measures ~-300 dB). + if (spec.image_zeros) { + for (int k = 1; k <= 3; ++k) { + EXPECT_LT(response_db(h, spec.num_phases, fs, static_cast(k) * fs), -150.0) + << "missing transmission zero at " << k << "*fs"; + } + } + } + + TEST(Kaiser, FastPrototypeMeetsSpec) { + check_prototype_meets_spec(k_fast, 48000.0); + } + + TEST(Kaiser, BalancedPrototypeMeetsSpec) { + check_prototype_meets_spec(k_balanced, 48000.0); + } + + TEST(Kaiser, TransparentPrototypeMeetsSpec) { + check_prototype_meets_spec(k_transparent, 48000.0); + } + + TEST(Kaiser, EconomyPrototypeMeetsSpec) { + check_prototype_meets_spec(k_economy, 48000.0); + } + + // The compensated presets must also hold their specs at scaled rates (the + // 16 kHz deployment path): normalized design, same numbers. + // The k*fs transmission zeros ARE branch-DC uniformity, stated in the + // frequency domain: with exact zeros, every polyphase branch's coefficient + // sum is identical (measured spread 1.8e-15 -- machine epsilon -- vs 4.7e-6 + // for the plain fast design, whose spread is its stopband leakage at fs). + TEST(Kaiser, CompensatedBranchSumsAreUniform) { + const proto_spec spec = k_balanced; + const std::size_t phases = std::bit_ceil(spec.num_phases); + std::vector h(phases * spec.taps_per_phase); + design_prototype_compensated(h, phases, (spec.passband_hz + spec.stopband_hz) / 48000.0, + kaiser_beta(spec.stopband_atten_db), spec.passband_hz / 48000.0); + double lo = 1e9; + double hi = -1e9; + for (std::size_t p = 0; p < phases; ++p) { + double sum = 0.0; + for (std::size_t t = 0; t < spec.taps_per_phase; ++t) { + sum += h[t * phases + p]; + } + lo = std::min(lo, sum); + hi = std::max(hi, sum); + } + EXPECT_LT(hi - lo, 1e-12); + EXPECT_NEAR(lo, 1.0, 1e-9); + } + + TEST(Kaiser, CompensatedSpecsHoldAt16k) { + check_prototype_meets_spec(k_balanced.scaled_to(16000.0), 16000.0); + check_prototype_meets_spec(k_economy.scaled_to(16000.0), 16000.0); + } + + // Non-power-of-two phase counts are first-class here (RatioTap designs at + // L = 147 and 160): the designer must meet spec without the bit_ceil + // rounding its historical consumer applied. 70 dB / 19->22.05 kHz is + // RatioTap's economy-profile shape for the 48->44.1 direction. + TEST(Kaiser, RationalPhaseCountMeetsSpec) { + const std::size_t phases = 147; + const double fs = 48000.0, pass = 19000.0, stop = 22050.0, atten = 70.0; + const std::size_t taps = estimate_taps(atten, (stop - pass) / fs); + std::vector h(phases * taps); + design_prototype(h, phases, (pass + stop) / fs, kaiser_beta(atten)); + for (double f = 0.0; f <= pass; f += 500.0) { + EXPECT_NEAR(response_db(h, phases, fs, f), 0.0, 0.05) << "passband deviation at " << f << " Hz"; + } + for (double f = stop; f <= 3.0 * fs; f += 250.0) { + EXPECT_LT(response_db(h, phases, fs, f), -(atten - 1.0)) << "stopband leakage at " << f << " Hz"; + } + } + + TEST(Kaiser, SolveDenseSolvesKnownSystem) { + // 3x3 that requires pivoting (zero leading entry); the solver works in + // place, so verify by residual against saved copies. + const std::vector m0{0.0, 2.0, 1.0, /**/ 1.0, 1.0, 1.0, /**/ 2.0, 0.0, 3.0}; + const std::vector rhs0{4.0, 6.0, 5.0}; + std::vector m = m0, rhs = rhs0, x(3); + solve_dense(m, rhs, x, 3); + for (int r = 0; r < 3; ++r) { + const double v = m0[3 * r] * x[0] + m0[3 * r + 1] * x[1] + m0[3 * r + 2] * x[2]; + EXPECT_NEAR(v, rhs0[static_cast(r)], 1e-12) << "row " << r; + } + } + +} // namespace diff --git a/tests/test_quantize.cpp b/tests/test_quantize.cpp new file mode 100644 index 0000000..df3a9d3 --- /dev/null +++ b/tests/test_quantize.cpp @@ -0,0 +1,101 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for row-sum-preserving quantization, ported from the +// row-sum checks in SampleRateTap's test_fixed_point.cpp (there stated +// against the assembled polyphase bank; here against the utility itself, +// over rows of a genuinely designed prototype). + +#include +#include +#include + +#include + +#include "tap/dsp/kaiser.h" +#include "tap/dsp/quantize.h" + +namespace { + + using tap::dsp::quantize_row_preserving_sum; + using tap::dsp::sample_traits; + + // A realistic source: the RatioTap-shaped design (L = 160, plain Kaiser) + // whose branches this utility will quantize in production. + std::vector designed_prototype(std::size_t phases, std::size_t taps) { + std::vector h(phases * taps); + tap::dsp::design_prototype(h, phases, (19000.0 + 22050.0) / 48000.0, tap::dsp::kaiser_beta(70.0)); + return h; + } + + template + void check_rows_sum_exact() { + constexpr std::size_t k_phases = 160; + constexpr std::size_t k_taps = 32; + const auto proto = designed_prototype(k_phases, k_taps); + using coeff = typename sample_traits::coeff; + const double scale = sample_traits::k_coeff_scale; + + std::vector row(k_taps); + std::vector q(k_taps); + for (std::size_t p = 0; p < k_phases; ++p) { + double exact = 0.0; + for (std::size_t t = 0; t < k_taps; ++t) { + row[t] = proto[t * k_phases + p]; + exact += row[t] * scale; + } + quantize_row_preserving_sum(row, q); + std::int64_t sum = 0; + for (std::size_t t = 0; t < k_taps; ++t) { + sum += q[t]; + // Each tap stays within the rounding step plus at most the + // residual corrections: never grossly redistributed. + EXPECT_LE(std::abs(static_cast(q[t]) - row[t] * scale), 1.5) << "phase " << p << " tap " << t; + } + ASSERT_EQ(sum, std::llround(exact)) << "phase " << p; + } + } + + TEST(Quantize, RowSumsAreExactQ15) { + check_rows_sum_exact(); + } + + TEST(Quantize, RowSumsAreExactQ31) { + check_rows_sum_exact(); + } + + TEST(Quantize, FloatIsPlainConversion) { + const std::vector row{0.25, -0.125, 1.0, -0.9999, 0.0}; + std::vector q(row.size()); + quantize_row_preserving_sum(row, q); + for (std::size_t t = 0; t < row.size(); ++t) { + EXPECT_FLOAT_EQ(q[t], static_cast(row[t])); + } + } + + // The consequence the correction exists for: a unity-DC row quantizes to a + // row whose coefficient sum is exactly the format's unity, so DC gain + // through any phase deviates by at most one output LSB. + TEST(Quantize, UnityDcRowSumsToFormatUnity) { + constexpr std::size_t k_taps = 48; + // A smooth row normalized to sum exactly 1.0 in double. + std::vector row(k_taps); + double sum = 0.0; + for (std::size_t t = 0; t < k_taps; ++t) { + const double u = (static_cast(t) - 23.5) / 24.0; + row[t] = std::exp(-4.0 * u * u); + sum += row[t]; + } + for (auto& v : row) { + v /= sum; + } + std::vector q(k_taps); + quantize_row_preserving_sum(row, q); + std::int64_t qsum = 0; + for (const auto c : q) { + qsum += c; + } + EXPECT_EQ(qsum, 16384); // Q1.14 unity, exactly + } + +} // namespace diff --git a/tests/test_sample_traits.cpp b/tests/test_sample_traits.cpp new file mode 100644 index 0000000..57aff55 --- /dev/null +++ b/tests/test_sample_traits.cpp @@ -0,0 +1,111 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the sample-format traits, ported and extended from +// SampleRateTap's test_fixed_point.cpp. Every number here is a documented +// contract point (Q formats, rounding mode, saturation, accumulator +// pre-shift); changing one is a breaking change for every consumer. + +#include +#include + +#include + +#include "tap/dsp/sample_traits.h" + +namespace { + + using q15 = tap::dsp::sample_traits; + using q31 = tap::dsp::sample_traits; + using f32 = tap::dsp::sample_traits; + + TEST(SampleTraits, CoefficientConversionRoundsAndSaturates) { + EXPECT_EQ(q15::make_coeff(0.0), 0); + EXPECT_EQ(q15::make_coeff(1.0), 16384); // Q1.14 + EXPECT_EQ(q15::make_coeff(-1.0), -16384); + EXPECT_EQ(q15::make_coeff(10.0), 32767); // saturates + EXPECT_EQ(q15::make_coeff(-10.0), -32768); + EXPECT_EQ(q31::make_coeff(1.0), 1073741824); // Q1.30 + EXPECT_EQ(q31::make_coeff(10.0), 2147483647); // saturates + EXPECT_FLOAT_EQ(f32::make_coeff(0.25), 0.25f); // unity scale, plain cast + } + + TEST(SampleTraits, RoundSatRoundsHalfAwayFromZero) { + using tap::dsp::detail::round_sat; + EXPECT_EQ(round_sat(0.5), 1); + EXPECT_EQ(round_sat(-0.5), -1); + EXPECT_EQ(round_sat(0.49), 0); + EXPECT_EQ(round_sat(-0.49), 0); + EXPECT_EQ(round_sat(1e9), std::numeric_limits::max()); + EXPECT_EQ(round_sat(-1e9), std::numeric_limits::min()); + } + + TEST(SampleTraits, FinalizeSaturates) { + // Far beyond full scale in the accumulator domain must clamp, not wrap. + EXPECT_EQ(q15::finalize(std::int64_t{1} << 40), 32767); + EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 40)), -32768); + EXPECT_EQ(q31::finalize(std::int64_t{1} << 60), 2147483647); + EXPECT_EQ(q31::finalize(-(std::int64_t{1} << 60)), -2147483648LL); + } + + TEST(SampleTraits, FinalizeRoundsHalfUp) { + // Q29 -> Q15: the rounding constant is 2^13; exactly half rounds up, + // one below half rounds down, and the same holds on the negative side + // (round-half-up, not half-away: -0.5 LSB lands on 0). + EXPECT_EQ(q15::finalize((std::int64_t{1} << 13)), 1); + EXPECT_EQ(q15::finalize((std::int64_t{1} << 13) - 1), 0); + EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 13)), 0); + EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 13) - 1), -1); + // Q45 -> Q31 uses the identical constant. + EXPECT_EQ(q31::finalize((std::int64_t{1} << 13)), 1); + EXPECT_EQ(q31::finalize((std::int64_t{1} << 13) - 1), 0); + } + + TEST(SampleTraits, Q15MacIsExactInInt64) { + // Q0.15 x Q1.14 products are exact in int32 and summed without + // intermediate rounding: full-scale sample times unity coefficient. + const auto acc = q15::mac(0, std::int16_t{32767}, q15::make_coeff(1.0)); + EXPECT_EQ(acc, std::int64_t{32767} * 16384); + // finalize returns the sample: one rounding, at the end. + EXPECT_EQ(q15::finalize(acc), 32767); + } + + TEST(SampleTraits, Q31MacPreShiftsSixteenBits) { + // Q0.31 x Q1.30 = Q61 would overflow int64 after ~48 accumulations, so + // each product drops 16 bits before the add (Q45). Contract: exactly + // 16, no more (precision), no fewer (headroom). + const auto acc = q31::mac(0, std::int32_t{1} << 30, std::int32_t{1} << 30); + EXPECT_EQ(acc, std::int64_t{1} << 44); // (2^60) >> 16 + // Headroom check: ~80 full-scale accumulations stay inside int64. + std::int64_t big = 0; + for (int i = 0; i < 80; ++i) { + big = q31::mac(big, std::numeric_limits::max(), std::numeric_limits::max()); + } + EXPECT_GT(big, 0); // no wrap + EXPECT_EQ(q31::finalize(big), std::numeric_limits::max()); + } + + TEST(SampleTraits, FloatMacAccumulatesInDouble) { + static_assert(std::is_same_v); + // A double accumulator must hold a contribution a float one would drop: + // 1.0 + 2^-30 survives in double, vanishes in float. + const double acc = f32::mac(f32::mac(0.0, 1.0f, 1.0f), 0x1p-15f, 0x1p-15f); + EXPECT_GT(acc, 1.0); + EXPECT_DOUBLE_EQ(acc, 1.0 + 0x1p-30); + } + + TEST(SampleTraits, SilenceIsZero) { + EXPECT_EQ(q15::silence(), 0); + EXPECT_EQ(q31::silence(), 0); + EXPECT_EQ(f32::silence(), 0.0f); + } + + // The core concept is the kernels' requirement set; all three shipping + // formats satisfy it (also statically asserted in the header). + static_assert(tap::dsp::sample_type); + static_assert(tap::dsp::sample_type); + static_assert(tap::dsp::sample_type); + static_assert(!tap::dsp::sample_type); // no specialization on purpose: + // double is the golden-model/accumulator domain, not an I/O sample format. + +} // namespace From 728e842f32d43d5efcc02d1c491342cf00df4433 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 12:17:19 +0000 Subject: [PATCH 2/2] Fix a sign-conversion error AppleClang catches and GCC does not The solve_dense verification loop indexed the saved matrix with an int: clang's -Wconversion implies -Wsign-conversion (GCC's does not in C++), so the macos-vdsp job failed under -Werror while linux-ooura passed. Index with std::size_t. Verified locally with a clang -Werror build of the full test suite, which also covers the TUs the failed CI run never reached. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Ldeq57sBySx2nFTQsGcQB6 --- tests/test_kaiser.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_kaiser.cpp b/tests/test_kaiser.cpp index 813cd20..ee866dd 100644 --- a/tests/test_kaiser.cpp +++ b/tests/test_kaiser.cpp @@ -187,9 +187,9 @@ namespace { const std::vector rhs0{4.0, 6.0, 5.0}; std::vector m = m0, rhs = rhs0, x(3); solve_dense(m, rhs, x, 3); - for (int r = 0; r < 3; ++r) { + for (std::size_t r = 0; r < 3; ++r) { const double v = m0[3 * r] * x[0] + m0[3 * r + 1] * x[1] + m0[3 * r + 2] * x[2]; - EXPECT_NEAR(v, rhs0[static_cast(r)], 1e-12) << "row " << r; + EXPECT_NEAR(v, rhs0[r], 1e-12) << "row " << r; } }