function yhat = Rdwt(dwt,clo,jo,dec);

% reconstruction from a single dwt
% input:
%	dwt	discrete wavelet transform with length a power of 2
%		last half coefs are highest freq difs (mcoefs),
%		second 1/4 are next highest, ... , first 2^jo are
%		the new "avgs"
%	clo	lo-pass filter values as made by MakeFilter.m
%	jo	row # of lowest freq difs = 2^jo
%	dec	order of decimation for dwt to be constructed
% output:
%	yhat	reconstructed vector
%

n = length(dwt);
je = nextpow2(n);
nfil = length(clo);

% Check -- is length of dwt a power of 2?
if (pow2(je)-n) ~= 0;
  disp('Length of vector is not a power of 2 -- cannot reconstruct.');
  return;
end;

% Is jo large enough?
if jo == -1
  jou = ceil(log2(nfil))-1 ;
elseif 2^(jo+1) < nfil
  disp('jo too small for this filter');
else jou = jo;
end;

% Find hi-pass filter
chi = -( (-1).^(1:nfil) ).*clo ;
chi = flipud(chi(:));

oldavg = dwt(1:2^jou);

k = length(dec);

for j = jou:(je-1)
  nj = pow2(j);
  e = dec(k);

  dif = dwt((nj+1): 2*nj);
  dif = [0;  ZeroInsert(dif)];
  avg = [ZeroInsert(oldavg);  0];

  na = length(avg);
  nd = length(dif);

%  na = length(avg)-1 ;
%  nd = length(dif)-1 ;
%  avg = avg((2-e) : (na+1-e));
%  dif = dif((1+e) : (nd+e));

  if nfil <= na ;
     avgpad = [avg((na+1-nfil):na);  avg];
  else
    z = zeros(nfil,1);
    for i=1:nfil
	imod = 1 + rem((nfil*na-nfil+i-1),n);
	z(i) = avg(imod);
    end;
    avgpad = [z; avg];
  end;

  if nfil < nd ; 
    difpad = [dif;  dif(1:nfil)];
  else
    z = zeros(nfil,1);
    for i = 1:nfil
      imod = 1 + rem(i-1,nd);
      z(i) = dif(imod);
    end;
    difpad = [dif; z];
  end;

  nlo = filter(clo,1,avgpad);
  nhi = filter(chi,1,difpad);

  nlo = nlo((nfil+1): (na+nfil));
  nhi = nhi(nfil:(nd+nfil-1));

  oldavg = nlo + nhi ;
  oldavg = downshift(oldavg,e);
  k = k-1 ;
end;

yhat = oldavg;

