function w = invWavTransf(v)

% Calculate the inverse wavelet transform, using Daubechies fourth-order
% coefficients, of a vector of length 2^12, with a 1 at indexes determines
% by the vector v of integers from 1 to 2^12.

% create my input vector
nbpts = 2^12;
w = zeros([1,nbpts]);
for i=1:length(v),
	w(v(i)) = 1;
	w(v(i)) = 1;
end

% get the number of steps in the decomposition
nsteps = log(nbpts)/log(2) - 1; 

% inverse wavelet transform with daubechies fourth order
% initialize daubechies coefficients
c0 = (1+sqrt(3))/(4*sqrt(2));
c1 = (3+sqrt(3))/(4*sqrt(2));
c2 = (3-sqrt(3))/(4*sqrt(2));
c3 = (1-sqrt(3))/(4*sqrt(2));

% double convolution and shuffling
% Scaling function: X = c2*xi + c1*xi+1 + c0*xi+2 + c3*xi+3
% Wavelet function: W = c3*xi - c0*xi+1 + c1*xi+2 - c2*xi+3
% w contains the output

for j=1:nsteps,
   % shuffling X and W's
   xstep = w(1:2^j);
   wstep = w(2^j + 1:2^(j+1));
   w(1:2:2^(j+1)) = xstep;
   w(2:2:2^(j+1)) = wstep;

   % first step convolutions
   xstep(1) = c0*w(1) + c3*w(2) + c2*w(nbpts-1) + c1*w(nbpts);
   xstep(2) = c1*w(1) - c2*w(2) + c3*w(nbpts-1) - c0*w(nbpts);

   % we loop on the changed part of the signal for this step
   for i=3:2:2^(j+1),
   	xstep(i)   = c2*w(i-2) + c1*w(i-1) + c0*w(i) + c3*w(i+1);
   	xstep(i+1) = c3*w(i-2) - c0*w(i-1) + c1*w(i) - c2*w(i+1);
	end

   % update w with new step
   w(1:2^(j+1)) = xstep;
end

% draw inverse wavelet transform
plot(w);

