function f = wv_peakfreq(a, dt, flo, fhi)
% WV_PEAKFREQ  The frequency of the strongest partial of a between flo and fhi.
%
% This is how every frequency in the tutorial is measured, from a recorded
% signal and nothing else: taper the record with wv_window, find the largest
% peak of its zero-padded DFT inside the band, and then polish the peak with
% Newton's method on |X(f)|^2, using the exact derivatives of the windowed
% discrete-time Fourier transform.  For a clean partial well separated from
% its neighbours the result is good to about 1e-10 relative; the method used
% on it never needs to be.
%
% a is a real record sampled every dt seconds; the band is in hertz.
a = a(:);
n = numel(a);
w = wv_window(n);
b = w .* (a - mean(a));
nfft = 2 ^ nextpow2(4 * n);
B = abs(fft(b, nfft));
df = 1 / (nfft * dt);
k0 = max(2, floor(flo / df) + 1);
k1 = min(floor(nfft / 2), ceil(fhi / df) + 1);
[~, i] = max(B(k0:k1));
f = (k0 + i - 2) * df;
tt = (0:n-1)' * dt;
for it = 1:30
    e = exp(-2i * pi * f * tt);
    X = sum(b .* e);
    X1 = sum(b .* e .* (-2i * pi * tt));
    X2 = sum(b .* e .* (-4 * pi ^ 2 * tt .^ 2));
    g = 2 * real(conj(X) * X1);
    h = 2 * real(conj(X1) * X1 + conj(X) * X2);
    if h >= 0
        break;
    end
    step = -g / h;
    if abs(step) > df
        step = sign(step) * df;
    end
    f = f + step;
    if abs(step) < 1e-13 * f
        break;
    end
end
end