function [y, peakPa] = woodblock_output(P, p, oversample)
% WOODBLOCK_OUTPUT Anti-aliased decimation and output conditioning.
% Kaiser-windowed-sinc FIR lowpass at 0.45 * fs_lo, decimate by oversample,
% pre-roll, tail fade, then either normalize to peak_dbfs or divide by
% full_scale_pa (woodblock.md section 8).

peakPa = max(abs(p));
R = oversample;
T = length(p);

if R > 1
  taps = 12 * R + 1;
  M2 = (taps - 1) / 2;
  fcNorm = 0.45 / R;                     % cycles per input sample
  m = (0:taps - 1) - M2;
  h = sin(2 * pi * fcNorm * m) ./ (pi * m);
  h(M2 + 1) = 2 * fcNorm;                % sinc value at m = 0
  r = m / M2;
  beta = 8;
  h = h .* (wb_i0(beta * sqrt(max(0, 1 - r.^2))) / wb_i0(beta));
  h = h / sum(h);                        % unity DC gain
  pf = conv(p, h);                       % pf(i) = sum_k p(k) h(i-k+1)
  outLen = floor(T / R);
  y = pf((0:outLen - 1) * R + M2 + 1);   % centered polyphase pick
else
  y = p;
end

% Pre-roll, then truncate to the nominal duration
padN = floor(P.pre_delay_ms * 1e-3 * P.sample_rate);
total = floor(P.duration_s * P.sample_rate);
if padN > 0
  y = [zeros(1, padN) y];
end
y = y(1:min(length(y), total));

% Squared fade on the tail
fade = min(256, floor(length(y) / 8));
if fade > 1
  y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
end

if P.normalize ~= 0
  pk = max(abs(y));
  if pk > 0
    y = y / pk * 10^(P.peak_dbfs / 20);
  end
else
  y = y / P.full_scale_pa;
end
end

function v = wb_i0(x)
% Modified Bessel function of the first kind, order zero (series expansion).
v = ones(size(x));
term = ones(size(x));
hx = x / 2;
for k = 1:40
  term = term .* (hx / k).^2;
  v = v + term;
end
end