% Approximation of a functions by the scaling function and its translates. 

p = input('Order of wavelet, p (defaults to 1) = ');
if isempty(p)
  p = 1;
end
h = daub(2*p);
%h = [-1 0 9 16 9 0 -1 0]'/16;

N = length(h);
[x1,phi] = phivals(h,4);

nx = 2^4;
L = 12;
n = L * nx;
x = (0:n)/n*6;
sf = [phi; zeros((L-N+1)*2^4,1)];

for num = 1:5

  y = zeros(size(sf));
  tmp = ((-N+2)*nx:(-N+2+32)*nx-1)'/n;
  v = x'/6;
  if num == 1
    ck = [zeros(N-2,1); h; zeros(L,1)];
  elseif num == 2
    f = ones(n+1,1);
  elseif num == 3
    f = v;
    ck = scalecoeffs(tmp,32,h,0,4);
  elseif num == 4
    f = 4*v.*v-4*v+1;
    ck = scalecoeffs(4*tmp.*tmp-4*tmp+1,32,h,0,4);
  elseif num == 5
    f = -6*v.*v.*v+9*v.*v-3*v;
    ck = scalecoeffs(-6*tmp.*tmp.*tmp+9*tmp.*tmp-3*tmp,32,h,0,4);
  end

  clf
  minval = 0;
  maxval = 0;
  for k = -N+2:L-1
    if num == 2
      g = eoshift(sf,k*nx);
    else
      g = ck(k+N-1)*eoshift(sf,k*nx);
    end
    hold on
    plot(x,g,':')
    hold off
    y = y + g;
    minval = min(minval,min(g));   
    maxval = max(maxval,max(g));    
  end
  hold on
  plot(x,y)
  hold off

  minval = min(0,min(y));
  maxval = max(y);
  s = maxval - minval;
  minval = minval - 0.2 * s;
  maxval = maxval + 0.2 * s;
  axis([min(x) max(x) minval maxval])
  xlabel('x')
  ylabel('f(x)')
  if num == 1
    title('Representation of a scaling function by its translates')
    v = axis;
    v(2) = 6;
    axis(v);
  elseif num == 2
    title('Representation of a constant function by translates of a scaling function')
  elseif num == 3
    title('Representation of a linear function by translates of a scaling function')
  elseif num == 4
    title('Representation of a quadratic by translates of a scaling function')
  elseif num == 5
    title('Representation of a cubic by translates of a scaling function')
  end

  if num > 1
    pause
    plot(x(1:n),f(1:n)-y(1:n))
    xlabel('x')
    ylabel('f(x)-f_{approx}(x)')
    title('Approximation error')
  end
  pause

end