1function [y, peakPa] = woodblock_output(P, p, oversample)
2% WOODBLOCK_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 (woodblock.md section 8).
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 .* (wb_i0(beta * sqrt(max(0, 1 - r.^2))) / wb_i0(beta));
21 h = h / sum(h); % unity DC gain
22 pf = conv(p, h); % pf(i) = sum_k p(k) h(i-k+1)
23 outLen = floor(T / R);
24 y = pf((0:outLen - 1) * R + M2 + 1); % centered polyphase 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);
31total = floor(P.duration_s * P.sample_rate);
32if padN > 0
33 y = [zeros(1, padN) y];
34end
35y = y(1:min(length(y), total));
37% Squared fade on the tail
38fade = min(256, floor(length(y) / 8));
39if fade > 1
40 y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
41end
43if P.normalize ~= 0
44 pk = max(abs(y));
45 if pk > 0
46 y = y / pk * 10^(P.peak_dbfs / 20);
47 end
48else
49 y = y / P.full_scale_pa;
50end
51end
53function v = wb_i0(x)
54% Modified Bessel function of the first kind, order zero (series expansion).
55v = ones(size(x));
56term = ones(size(x));
57hx = x / 2;
58for k = 1:40
59 term = term .* (hx / k).^2;
60 v = v + term;
61end
62end