function [fr, db] = wv_spectrum(x, fs, fmax, npts)
% WV_SPECTRUM  Magnitude spectrum of x in decibels, for plotting.
%
% Tapered with wv_window and zero-padded, then reduced to npts points up to
% fmax by keeping the largest value in each bin of the plot, so that narrow
% peaks survive the thinning.  0 dB is the largest peak.
x = x(:);
n = numel(x);
w = wv_window(n);
nfft = 2 ^ nextpow2(2 * n);
X = abs(fft(w .* x, nfft));
X = X(1:nfft/2+1);
X = X / max(max(X), 1e-300);
f = (0:nfft/2)' * fs / nfft;
fr = linspace(0, fmax, npts);
db = -200 * ones(1, npts);
bw = fr(2) - fr(1);
idx = floor(f / bw + 0.5) + 1;
% real() because log10 is typed as possibly complex, which would keep this
% loop off numbl's compiled path
v = real(20 * log10(max(X, 1e-300)));
for k = 1:numel(f)
    j = idx(k);
    if j >= 1 && j <= npts && v(k) > db(j)
        db(j) = v(k);
    end
end
db(db < -160) = -160;
end