MakeFilter(WaveType,m,n,disp)

function [h,g,rh,rg,d,rd] = MakeFilter(WaveType,m,n,disp)

% all-purpose wavelet generator based on freeware toolboxes
% delays are computed to synchronize all cascades
% Analysis filters h and g satisfy h[n] = h1[-n] and g[n] = g1[-n]
% where h1 and g1 define the scaling equations on the analysis
% scaling function and wavelet
% d is a left delay which makes the analysis filters causal.
% rd is a right delay which makes the synthesis filters causal.
% d+rd are computed in order to synchronize the cascades.
% For analysis filters, the z-transform of h and g are
% h(z) = z^d sum h[k+1]z^(-k) and g(z) = z^d sum g[k+1]z^(-k) 
% For reconstruction filters, the z-transform of rh and rg are
% rh(z) = z^rd sum rh[k+1]z^(-k) and rg(z) = z^rd sum rg[k+1]z^(-k)

if nargin <= 3
	disp = 0;
end
% disp(lib);


switch WaveType
case {'Daubechies'}
	load(['Daub' num2str(2*m)]);
	l = length(h);
	% for orthogonal filters with an even length, wavelets can be shifted right
	% by a length of l/2-1, which amounts to a renumbering of the basis and coefficients.
	% If the support of rh is [0,l-1], then the support of the shifted rg is still [0,l-1],
	% and the supports of h and g are [1-l,0], which accounts for the delay d = l-1.
	if disp
		doDisp(h,g),
	end
	d = l-1;
	rd = 0;
case {'Symlets'}
	load(['Syml' num2str(2*m)]);
	l = length(h);
	% for orthogonal filters with an even length, wavelets can be shifted right
	% by a length of l/2-1, which amounts to a renumbering of the basis and coefficients.
	% If the support of rh is [0,l-1], then the support of the shifted rg is still [0,l-1],
	% and the supports of h and g are [1-l,0], which accounts for the delay d = l-1.
	if disp
		doDisp(h,g);
	end
	d = l-1;
	rd = 0;
case {'Splines'}
	[h,g,rh,rg] = splineFilterBanks(m,n);
	if disp
		doDisp(h,g,rh,rg);
	end
	if mod(m,2) == 0
		rl = m/2; % support of rh = [-rl,rl], hence supp(g) = [-rl-1,rl+1]
		l = n+rl-1; % support of h = [-l,l], hence supp(rg) = [1-l,1+l]
		d = l; % shift h to make supp = [0,2l]
		rd = l-1; % shift rg to make supp = [0,2l]
		g = [zeros(1,l-rl-1) g]; % then supp(g) = [l-rl-1,l+rl-1]
		rh =[zeros(1,l-rl-1) rh]; % then supp(rh) = [l-rl-1,l+rl-1]
	else
		rl = (m-1)/2; % support of rh = [-rl,rl+1], hence supp(g) = [-rl-1,rl]
		l = n+rl-1; % support of h = [-l-1,l], hence supp(rg) = [-l,1+l]
		d = l+1; % shift h to make supp = [0,2l+1]
		rd = l; % shift rg to make supp = [0,2l+1]
		g = [zeros(1,l-rl) g];% then supp(g) = [l-rl,rl+l+1]
		rh =[zeros(1,l-rl) rh]; % then supp(rh) = [l-rl,rl+l+1]
	end
	% cone of influence: observe that the delays d are rd have a length equal to the half
	% width of the cone of influence! operating within these delays should be enough to 
	% take care regularization/regularity analysis
otherwise
	error([WaveType ' wavelet not supported']);
end


function [s,w] = getWave(h,g)
[s,w] = simuWavelet(h,g,6);

function doDisp(h,g,rh,rg)
if nargin <= 3
	[s,w] = getWave(h,g);
	subplot(2,1,1);
	plot(s);
	title('Scaling function');
	subplot(2,1,2);
	plot(w);
	title('Wavelet');
else
	subplot(2,1,1);
	[s,w] = getWave(h,g);
	subplot(2,2,1);
	plot(s);
	title('Analysis scaling function');
	subplot(2,2,3);
	plot(w);
	title('Analysis wavelet'); 
	[rs,rw] = getWave(rh,rg);
	subplot(2,2,2);
	plot(rs);
	title('Synthesis scaling function');
	subplot(2,2,4);
	plot(rw);
	title('Synthesis wavelet');
end