function [h,g] = hwlet(K,L)
% Hilbert transform pair of orthogonal wavelet bases
% h, g - scaling filters of length 2*(K+L)
% K - number of zeros at z=-1
% L - degree of fractional delay

n = 0:L-1;
t = 1/2;
d = cumprod([1, (L-n).*(L-n-t)./(n+1)./(n+1+t)]);
s1 = binom(2*K,0:2*K);
s2 = conv(d,d(end:-1:1));
s = conv(s1,s2);
M = K+L;
C = convmtx(s',2*M-1);
C = C(2:2:end,:);
b = zeros(2*M-1,1);
b(M) = 1;
r = (C\b)';
q = sfact(r);
f = conv(q,binom(K,0:K));
h = conv(f,d);
g = conv(f,d(end:-1:1));