function s = wv_audio(x, fs, scale)
% WV_AUDIO  A signal packed for the page's audio player.
%
% The samples are divided by `scale` and sent as 16-bit integers, which keeps
% the JSON about a third the size of full-precision numbers and loses nothing
% that a loudspeaker could reproduce.  Pass the same scale to every clip in a
% comparison, so that "louder" on the page means louder in the computation;
% with no scale the clip is normalised to its own peak.  A short fade at each
% end keeps the player from clicking; it is applied to every clip alike.
%
% Browsers are only obliged to play sample rates from 8 kHz up, and a coarse
% simulation's own rate can be lower, so such a clip is first raised to a
% whole multiple of its rate by band-limited interpolation (zero-padding its
% spectrum), which adds nothing and changes nothing below its Nyquist
% frequency.
x = x(:)';
if fs < 8000
    k = ceil(8000 / fs);
    x = wv_upsample(x, k);
    fs = fs * k;
end
if nargin < 3 || scale <= 0
    scale = max(abs(x));
    if scale == 0
        scale = 1;
    end
end
n = numel(x);
nf = min(floor(n / 4), round(0.004 * fs));
if nf > 1
    r = 0.5 - 0.5 * cos(pi * (0:nf-1) / nf);
    x(1:nf) = x(1:nf) .* r;
    x(n-nf+1:n) = x(n-nf+1:n) .* fliplr(r);
end
y = round(29000 * x / scale);
y(y > 32767) = 32767;
y(y < -32767) = -32767;
s = struct();
s.fs = fs;
s.pcm = int16(y);
end

function y = wv_upsample(x, k)
n = numel(x);
X = fft(x);
Y = zeros(1, k * n);
h = floor(n / 2);
Y(1:h + 1) = X(1:h + 1);
Y(k * n - (n - h - 2):k * n) = X(h + 2:n);
if mod(n, 2) == 0
    % the Nyquist bin is shared between the two halves
    Y(h + 1) = X(h + 1) / 2;
    Y(k * n - h + 1) = X(h + 1) / 2;
end
y = real(ifft(Y)) * k;
end