function X = SCOMP(F,G,b,s,band)
%OMP Summary of this function goes here
% Detailed explanation goes here

[m,N] = size(F); % size(F), size(G) should be equal
itermax = s;
%%

r = b;
X = zeros(2*N,1);

S = [];

Region = 1:N;

for iter=1:itermax
%step1
mu = abs(F'*r) + abs(G'*r);
[q,lambda] = max(mu(Region));

NewIndex = Region(lambda);
Region = setdiff(Region,NewIndex-band:NewIndex+band);

%step2
S = [S,NewIndex]; % display(Lda); %Lda = union(Lda,lambda);

F_S = F(:,S);
G_S = G(:,S);
A_S = [F_S G_S];

%step3, pseudo inverse
x = (A_S'*A_S)\(A_S'*b);

%step4, residual
r = b - A_S*x;

%dump
res = norm(r,2);
%display(sprintf('iter=%d, res=%g', iter, res))

X([S,N+S]) = x;
if(res<1e-14) % || length(S)==s)
return;
end
end