function [K, k] = cqkern(sr, bpo, f_min, f_max)
% [K, k, W] = cqkern(sr, bpo = 12, f_min = midifreq(12), f_max = sr/2)
%
% Calculate the spectral (K) and temporal (k) kernels for an efficient
% constant-Q transform [1].
%
% 1. An efficient algorithm for the calculation of a constant Q transform
%    Judith C. Brown and Miller S. Puckette
%    J. Acoust. Soc. Am. 92 (5), November 1992

error(nargchk(1, 4, nargin));
if nargin < 2 || isempty(bpo)
    % semitones
    bpo = 12;
end
if nargin < 3 || isempty(f_min)
    % C0, MIDI note 12
    f_min = 16.352;
end
if nargin < 4 || isempty(f_max)
    % Nyquist
    f_max = sr/2;
end

% Constants
r = 2^(1/bpo);
Q = 1 / (r - 1);
B = ceil(log(f_max/f_min) / log(r)); % number of bins
o_min = 2 * pi * f_min / sr; % o for omega (w for window)
o = o_min * r.^(0:B-1)'; % vector of omegas
N = round(2*pi*Q ./ o); % vector of Nkcqs
N0 = N(1); % Longest N

k = zeros(B, N0); % temporal kernel
for kcq = 1:B
    % centered window
    w = postpad(hamming(N(kcq)), N0)';
    k(kcq,:) = (w .* exp(-j*o(kcq)*(0:N0-1))) / N(kcq);
end

K = fft(k',N0)';

% Make it sparse.
%K(abs(K) < 0.15) = 0;
%K = sparse(K);
