1function [y, peakPa] = maraca_output(P, p, oversample)
2% MARACA_OUTPUT Anti-aliased decimation and output conditioning.
3% Kaiser-windowed-sinc FIR lowpass at 0.45 * fs_lo, decimate by oversample,
4% pre-roll, tail fade, then either normalize to peak_dbfs or divide by
5% full_scale_pa. Identical in structure to the woodblock output stage.
7peakPa = max(abs(p));
8R = oversample;
9T = length(p);
11if R > 1
12 taps = 12 * R + 1;
13 M2 = (taps - 1) / 2;
14 fcNorm = 0.45 / R; % cycles per input sample
15 m = (0:taps - 1) - M2;
16 h = sin(2 * pi * fcNorm * m) ./ (pi * m);
17 h(M2 + 1) = 2 * fcNorm; % sinc value at m = 0
18 r = m / M2;
19 beta = 8;
20 h = h .* (mr_i0(beta * sqrt(max(0, 1 - r.^2))) / mr_i0(beta));
21 h = h / sum(h); % unity DC gain
22 pf = conv(p, h);
23 outLen = floor(T / R);
24 y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen); % centered pick
25else
26 y = p;
27end
29% Pre-roll, then truncate to the nominal duration
30padN = floor(P.pre_delay_ms * 1e-3 * P.sample_rate);
31duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
32total = floor(duration * P.sample_rate);
33if padN > 0
34 y = [zeros(1, padN) y];
35end
36y = y(1:min(length(y), total));
38% Squared fade on the tail
39fade = min(256, floor(length(y) / 8));
40if fade > 1
41 y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
42end
44if P.normalize ~= 0
45 pk = max(abs(y));
46 if pk > 0
47 y = y / pk * 10^(P.peak_dbfs / 20);
48 end
49else
50 y = y / P.full_scale_pa;
51end
52end
54function v = mr_i0(x)
55% Modified Bessel function of the first kind, order zero (series expansion).
56v = ones(size(x));
57term = ones(size(x));
58hx = x / 2;
59for k = 1:40
60 term = term .* (hx / k).^2;
61 v = v + term;
62end
63end