function [b,w] = descent (h)

% Gradient descent
% Given a (noisy) barcode signal h
% Returns cleaned up 0-1 barcode b and width w.
% Assume a = 1.

% Maximum number of iterations.
itMAX = 10000;

% Initial guesses for constants.
[temp,N] = size(h);
a = 1;
w = 0.1;
alpha = 0.1 * N^2;
epsilon = 0.000001;

% Make sure h has even length.
[temp,N] = size(h);
if mod(N,2) == 1
  h = [h 0];
end;
[temp,N] = size(h);

% Initialize b.
b = ceil (h ./ max(h) - 0.5);
%b = rand([temp,N]);
x = [1:N] ./ N - 0.5;

subplot(3,1,1);
plot(h,'g');
title('Original signal h');

for it = 1:itMAX
  % Calculate convolution.
  g = a.*exp(-(x/w).^2);
  bg = conv(b,g);
  bg = bg(N/2+1:3*N/2);
  subplot(3,1,2);
  plot(1:length(bg),bg,1:length(bg),h);
  title('b*g vs. h');
  drawnow;
  
  % Calculate dJdb to minimize b.
  b_pad = [0 b 0];
  b_pad_x = (b_pad(2:N+2) - b_pad(1:N+1)) .* N;
  inside = b_pad_x ./ sqrt(epsilon + b_pad_x .^2);
  b_diff = (inside(2:N+1) - inside(1:N)) .* N;
  bg2 = -2 .* conv(g,h-bg);
  bg2 = bg2(N/2+1:3*N/2);
  dJdb = bg2 - alpha .* b_diff;
  
  % Calculate dJdw.
  Dg = 2*a / w^3 .* (x.^2) .* exp(-(x/w).^2);
  Dgb = conv(Dg,b);
  Dgb = Dgb(N/2+1:3*N/2);
  term = -2 .* (h - bg) .* Dgb;
  dJdw = sum(term);
  %w = w - lambda * dJdw;
  
  % Determine step size lambda (B-F p. 571).
  [temp,N] = size(dJdb);
  z = [dJdb dJdw]/ norm([dJdb dJdw]);   % Normalized.
  %z = [dJdb dJdw];    % Un-normalized
  lambda1 = 0;
  lambda3 = 100;       % Set bigger?
  g1 = J_fun3 (h,b,w,a,alpha,epsilon);
  g3 = J_fun3 (h,b - lambda3.*z(1:N),w - lambda3*z(N+1),a,alpha,epsilon);
  while g3 >= g1 & lambda3 >= epsilon/2
    lambda3 = lambda3 / 2;
	g3 = J_fun3 (h,b - lambda3.*z(1:N),w - lambda3*z(N+1),a,alpha,epsilon);
  end;
  lambda2 = lambda3 / 2;
  g2 = J_fun3 (h,b - lambda2.*z(1:N),w - lambda2*z(N+1),a,alpha,epsilon);
  h1 = (g2 - g1) / lambda2;
  h2 = (g3 - g2) / (lambda3 - lambda2);
  h3 = (h2 - h1) / lambda3;
  lambda0 = 0.5 * (lambda2 - h1 / h3);
  g0 = J_fun3 (h,b-lambda0.*z(1:N),w - lambda0*z(N+1),a,alpha,epsilon);
  if g0 < g3
    lambda = lambda0;
  else
    lambda = lambda3;
  end;
  
  %lambda = 0.01;    % Fixed lambda.
  b = b - lambda .* z(1:N);
  w = w - lambda * z(N+1);
  
  % Project b.
  [temp,N] = size(b);
  for i = 1:N
    if b(i) < 0
	  b(i) = 0;
	elseif b(i) > 1
	  b(i) = 1;
	end;
  end;
  
  J = J_fun3 (h,b,w,a,alpha,epsilon);
  
  subplot(3,1,3);
  stairs(b);
  axis([0 length(b) -0.1 1.1]);
  title(strcat('It#',num2str(it),' a=',num2str(a),' w=',num2str(w)));
  xlabel(strcat('J = ',num2str(J)));
  pause (0.01);
  
end;