% function [x,phi,phitilde,psi,psitilde] = biphivals(h0,h1,f0,f1,i)
% Generate biorthogonal scaling functions and their associated
% wavelets using the given filter coefficients
%               Kevin Amaratunga
%               4 August, 1993
%
%    h0, h1, f0, f1 = wavelet filters (from BIORFILT).
%    i = discretization parameter.  The number of points per integer
%        step is 2^i.  Thus, setting i = 0 gives the scaling function
%        and wavelet values at integer points.
%
function [x,phi,phitilde,psi,psitilde] = biphivals(h0,h1,f0,f1,i)
[N,dum] = size(f0);
if N == 1
  f0 = f0';
  f1 = f1';
  h0 = h0';
  h1 = h1';
  N = dum;
end
[tmp,dum] = size(h0);
if i < 0
  error('biphivals: i must be non-negative')
end
[m,n] = size(f0);
[tmp,dum] = size(h0);
if m ~= tmp
  error('biphivals: filters f0 and h0 must be the same length')
end
%
% Make sure the lowpass filters sum up to 2.
%
fac = 2/sum(h0);
h0 = flipud(h0) * fac;
h1 = flipud(h1) * fac;
f0 = f0 * fac;
f1 = f1 * fac;

cf0 = [f0; zeros(m,1)];
rf0 = [f0(1), zeros(1,m-1)];
tmp = toeplitz(cf0,rf0);
M = zeros(m,m);
M(:) = tmp(1:2:2*m*m-1);
M = M - eye(m);
M(m,:) = ones(1,m);
tmp = [zeros(m-1,1); 1];
phi = M \ tmp;		% Integer values of phi

ch0 = [h0; zeros(m,1)];
rh0 = [h0(1), zeros(1,m-1)];
tmp = toeplitz(ch0,rh0);
M = zeros(m,m);
M(:) = tmp(1:2:2*m*m-1);
M = M - eye(m);
M(m,:) = ones(1,m);
tmp = [zeros(m-1,1); 1];
phitilde = M \ tmp;		% Integer values of phitilde

if i > 0
  for k = 0:i-1
    p = 2^(k+1) * (m-1) + 1;	% No of rows in toeplitz matrix
    q = 2^k * (m-1) + 1;	% No of columns toeplitz matrix
    if (k == 0)
      cf00 = [f0; zeros(p-1-m,1)];
      cf0 = [cf00; 0];
      ch10 = [h1; zeros(p-1-m,1)];
      ch00 = [h0; zeros(p-1-m,1)];
      ch0 = [ch00; 0];
      cf10 = [f1; zeros(p-1-m,1)];
    else
      cf0 = zeros(p-1,1);
      cf0(:) = [1; zeros(2^k-1,1)] * cf00';
      cf0 = [cf0; 0];
      ch0 = zeros(p-1,1);
      ch0(:) = [1; zeros(2^k-1,1)] * ch00';
      ch0 = [ch0; 0];
    end
    rf0 = [cf0(1), zeros(1,q-1)];
    Tf0 = toeplitz(cf0,rf0);
    rh0 = [ch0(1), zeros(1,q-1)];
    Th0 = toeplitz(ch0,rh0);
    if k == i-1
      ch1 = [1; zeros(2^k-1,1)] * ch10';
      ch1 = ch1(:);
      ch1 = [ch1; 0];
      rh1 = [ch1(1), zeros(1,q-1)];
      Th1 = toeplitz(ch1,rh1);
      cf1 = [1; zeros(2^k-1,1)] * cf10';
      cf1 = cf1(:);
      cf1 = [cf1; 0];
      rf1 = [cf1(1), zeros(1,q-1)];
      Tf1 = toeplitz(cf1,rf1);
      psi = Tf1 * phi;
      psitilde = Th1 * phitilde;
    end
    phi = Tf0 * phi;
    phitilde = Th0 * phitilde;
  end
elseif i == 0
  ch10 = [h1; zeros(m-2,1)];
  ch1 = [ch10; 0];
  rh1 = [ch1(1), zeros(1,m-1)];
  Th1 = toeplitz(ch1,rh1);
  cf10 = [f1; zeros(m-2,1)];
  cf1 = [cf10; 0];
  rf1 = [cf1(1), zeros(1,m-1)];
  Tf1 = toeplitz(cf1,rf1);
  psi = Tf1 * phi;
  psi = psi(1:2:2*m-1);
  psitilde = Th1 * phitilde;
  psitilde = psitilde(1:2:2*m-1);
end

[a,b] = size(phi);
x = (0:a-1)' / 2^i;