function [y, peakAccel] = cymbal2_output(P, acc)
% CYMBAL2_OUTPUT Output conditioning: DC blocker on the raw pickup
% acceleration, Kaiser-windowed-sinc anti-alias decimation from the
% simulation rate to the output rate, tail fade, and normalization. Same
% structure as the sibling instruments' output stages.

peakAccel = max(abs(acc));
R = P.oversample;
T = length(acc);

% One-pole DC blocker (the pickup is an acceleration, but the contact
% leaves a small quasi-static offset while the washer relaxes)
fc = 15.0;
a1 = exp(-2 * pi * fc / (P.sample_rate * R));
p = zeros(1, T);
xp = 0; yp = 0;
for k = 1:T
  yv = acc(k) - xp + a1 * yp;
  xp = acc(k);
  yp = yv;
  p(k) = yv;
end

if R > 1
  taps = 12 * R + 1;
  M2 = (taps - 1) / 2;
  fcNorm = 0.45 / R;
  m = (0:taps - 1) - M2;
  hh = sin(2 * pi * fcNorm * m) ./ (pi * m);
  hh(M2 + 1) = 2 * fcNorm;
  r = m / M2;
  beta = 8;
  hh = hh .* (c2_i0(beta * sqrt(max(0, 1 - r.^2))) / c2_i0(beta));
  hh = hh / sum(hh);
  pf = conv(p, hh);
  outLen = floor(T / R);
  y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen);
else
  y = p;
end

% Squared fade on the tail
fade = min(2048, 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_accel;
end
end

function v = c2_i0(x)
% Modified Bessel function of the first kind, order zero (series)
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