function sp = postfilter(s,wavelet,filter_type,boundary_handling)

% POSTFILTER - multiwavelet postfiltering step
%
%        sp = postfilter(s,'wavelet','filter_type','boundary_handling')
%
% Supported wavelets and types are
%
%  'dghm':  'ap1' and 'ap2'
%           (orthogonal postfilters of Hardin-Roach,
%            preserving approximation orders 1 and 2, respectively)
%
%  'cl2', 'cl3', 'stt', 'hc': 'ap1'
%           (orthogonal postfilter preserving approximation order 1)
%
% Default for the filter type is 'ap1'.
%
% Supported boundary handling is only 'periodic' at this point, so this
% parameter is optional at the moment.

% Copyright (c) 2004 by Fritz Keinert (keinert@iastate.edu),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

if (nargin < 4)
    boundary_handling = 'periodic';
end

if (nargin < 3)
    filter_type = 'ap1';
end

% get postfilter coefficients
switch wavelet
 case 'dghm'
  m = 2;
  r = 2;
  switch filter_type
   case 'ap1'
    F = [-1+sqrt(2),1+sqrt(2);
	 1+sqrt(2),1-sqrt(2)]/sqrt(6);
   
   case 'ap2'
    F = [0.04967860804828,0.99158171438258,0.11942337067748,-0.00598315472909;
	 -0.00598315472909,-0.11942337067748,0.9915817143825,-0.04967860804828];
    
   otherwise
    error('unsupported postfiltering type for this wavelet');
  end
  
 case {'bat', 'daubbal'}
  error('This multiwavelet is balanced; it does not need a postfilter');

 case {'cl2', 'cl3', 'stt', 'hc'}
  m = 2;
  r = 2;
  switch filter_type
   case 'ap1'
    F = [1,1;1,-1] / sqrt(2);
   otherwise
    error('unsupported prefiltering type for this wavelet');
  end

 otherwise
  error('postfiltering is not supported yet for this wavelet');
end

filter_length = size(F,2)/size(F,1);

% check signal length, make sure it is a column vector
s = s(:);
n = length(s);
nsteps = round(n/r);
if (nsteps*r ~= n)
    error(['signal length is must be a multiple of ',num2str(r)]);
end
sp = zeros(size(s));

% set up boundary handling
switch lower(boundary_handling(1))
 case 'p'
  s = [s(end-r*(filter_length - 1)+1:end);s];
       
 otherwise
  error('unsupported boundary handling for this postfilter');
end

% do the postfiltering
for i = 1:nsteps
    sp((i-1)*r+1:i*r) = F * s((i-1)*r+1:(i+filter_length-1)*r);
end