function [x,f,cost] = lpfcsd(y, d, fc, lam0, lam1, Nit, mu)
% [x, f, cost] = lpfcsd(y, d, fc, lam0, lam1, Nit, mu)
% Simultaneous low-pass filtering and compound sparsity denoising 
%
% INPUT
%   y - noisy data
%   d - degree of filter is 2d (use d = 1, 2, or 3)
%   fc - cut-off frequency (normalized frequency, 0 < fc < 0.5)
%   lam0, lam1 - regularization parameters for x and diff(x)
%   Nit - number of iterations
%   mu - ADMM parameter
%
% OUTPUT
%   x - TV component
%   f - LPF component
%   cost - cost function history

% Ivan Selesnick
% Polytechnic Institute of New York University
% New York, USA
% March 2012
% Revised June 2012
% Revised Jan 2013

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

[A, B] = ABfilt(d, fc, N);
Id = @(x) x(d+1:N-d);
H = @(x) A\(B*x);                      % H: high-pass filter
G = mu*(A*A') + B*B';                  % G: banded matrix [sparse]
bn = nan + zeros(d, 1);     	       % bn : nan's to extend f to length N

v = zeros(N, 1);                       % initializations
d = zeros(N, 1);
b = (1/mu) * B'*((A*A')\(B*y));

for k = 1:Nit
    g = b + v - d;
    x = g - B' * (G \ (B*g));          % banded system solve (G)
    v = tvd(x + d, N, lam1/mu);        % TV denoising
    v = soft(v, lam0/mu);        
    v = v(:);
    d = d + x - v;
    cost(k) = lam0 * sum(abs(x)) + lam1 * sum(abs(diff(x))) + 0.5 * sum(abs(H(x-y)).^2);
end

f = y - x - [bn; H(y-x); bn];          % f : low-pass component

