function [ TF_l1, TF_l2, epsilon ] = sparsetfr( s, larg, epsilon )
% sparsetfr.m 
%
% Written by: Patrick Flandrin and Pierre Borgnat, CNRS-ENS de Lyon
% Email: {flandrin,pborgnat}@ens-lyon.fr
% Created: April 2015
%
% adapted from: 
%
% tveq_phantom_example.m
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005

% Computes and plots sparse time-frequency representations as described in:
% P. Flandrin and P. Borgnat, 
% "Time-Frequency Energy Distributions Meet Compressed Sensing,"
% IEEE Trans. on Signal Proc., Vol. 58, No. 6, pp. 2974-2982, 2010.
% Digital Object Identifier: 10.1109/TSP.2010.2044839  
%
% 1. input:
%    - s       : signal              : 
%    - larg    : mask size (must be even)
%    - epsilon : trade-off parameter (default = 5e-2)
% 2. output:
%    - TF_l2   : min l_2 reconstruction (backprojection)
%    - TF_l1   : min l_1 solution with l_2 ambiguity function 
%      approximation fixed by epsilon
%    - epsilon

% calls m-scripts from:
%
% - Time-Frequency Toolbox (tftb.nongnu.org)
% - l1-magic Toolbox (http://users.ece.gatech.edu/~justin/l1magic/)

if nargin ==2
    epsilon = 5e-2 ;
end

% normalization to unit energy
s = s/sqrt(sum(abs(s).^2)) ;

n = length(s) ;
N = n ;
    
% Wigner-Ville distribution
X = tfrwv(s) ;  
x = X(:) ;

% square mask in the central region of the ambiguity function
M = zeros(N) ;
M(1+n/2-larg/2:1+n/2+larg/2,1+n/2-larg/2:1+n/2+larg/2) = ones(larg+1,larg+1) ;

% upper half plane mask (not including origin)
Mh = zeros(N) ;
Mh = M;
Mh(N/2+2:N,:) = 0;
Mh(N/2+1,N/2+1:N) = 0;

M = ifftshift(M);
Mh = ifftshift(Mh);
mi = find(M);
mhi = find(Mh);

OMEGA = mhi;
A = @(z) A_fhp(z, OMEGA);
At = @(z) At_fhp(z, OMEGA, n);

% initialization
TF_l1 = zeros(n) ;
TF_l2 = zeros(n) ;

% measurements
y = A(x);

% min l_2 reconstruction (backprojection) -- initial guess
xbp = At(y);
TF_l2 = reshape(xbp,n,n);

% large scale l_1 minimization with l_2 minimization  [from l1qc_example]
xp2 = l1qc_logbarrier(xbp, A, At, y, epsilon, 1e-3, 25, 1e-8, 500);
TF_l1 = reshape(xp2, n, n);

end

