function [x, cost] = tvd_ncvx_ver1(y, lam, pen, Nit)
% [x, cost] = tvd_ncvx_ver1(y, lam, pen, Nit)
% Total variation denoising with non-convex penalty.
%
% INPUT
%   y - noisy signal
%   lam - regularization parameter
%   pen - penalty ('log', 'atan', 'L1')
%   phi - penalty function
%   wfun - function wfun(x) = x / phi'(x)
%   Nit - number of iterations
%
% OUTPUT
%   x - denoised signal
%   cost - cost function history

% Ivan Selesnick, selesi@nyu.edu, 2012
% Revised June 2014
% Revised Januar 2017 (simplify D matrix)

% The algorithm uses majorization-minimization
% and a fast solver for banded linear systems.

y = y(:);                                             % Ensure column vector
cost = zeros(1, Nit);                                 % Cost function history
N = length(y);

a = 1 / (4 * lam);

[phi, dphi, wfun, ds] = penalty_funs(pen, a);

I = speye(N);
D = I(2:N, :) - I(1:N-1, :);
DDT = D * D';

x = y;                                                % Initialization
Dx = D*x;
Dy = D*y;

for k = 1:Nit
    F = sparse(1:N-1, 1:N-1, wfun(Dx)/lam) + DDT;       % F : Sparse banded matrix    
    x = y - D'*(F\Dy);                                  % Solve banded linear system
    Dx = D*x;
    
    cost(k) = 0.5*sum(abs(x-y).^2) + lam*sum(phi(Dx));  % Save cost function history
end
