% This function recover the signal support (CTF) from the sample sequences.
% Expanding the inputs in factor q is performed if q>1
%ExSamples, ExTimeAxis
function [Support  Success] = RecoverSupport(Signal, MWC, Samples, Time_Axis, handles)

% get relevant parameters
fmax = getfield(Signal.Structure, 'fmax');
N=getfield(Signal.Structure, 'N');
fnyq = 2*fmax;
M = getfield(MWC, 'M');
m = getfield(MWC, 'm');
fp = fnyq/M;
Ts = getfield(MWC, 'Ts');
fs = 1/Ts;
L = M;
if (mod(L,2)==0)
    L=L+1;
end   % L is the number of specturm slices, namely the length of z[n]
L0 = (L-1)/2;


% compute the sensing matrix (same as done in the sampling stage), we
% repeat the computations so that the reconstruction function can smoothly
% interface with samples from the hardware (not numerical simulations).

S = getfield(MWC, 'SignPatterns');
DFTM = fft(eye(M));
F = [ conj(fliplr(DFTM(:,2:(L0+1))))    DFTM(:,1:(L0+1)) ];
theta = exp(-j*2*pi/L);
np = 1:L0;
nn = (-L0):1:-1;
dn = [   (1-theta.^nn)./(j*2*pi*nn)      1/L    (1-theta.^np)./(j*2*pi*np)  ];
D = diag(dn);
A = conj(S*F*D);

C = A;
y = Samples;

% Run CTF
Q = y * y';

% decompose Q to find frame V
NumDomEigVals= FindNonZeroValues(eig(Q),5e-8);
[V,d] = eig_r(Q,min(NumDomEigVals,2*N));
v = V*diag(sqrt(d));
[RecSupp] = RunOMP_forMB(v, C,N,0);
Support = sort(unique(RecSupp));

% Decide on success
if (is_contained(handles.S,Support)  && (rank(C(:,Support)) == length(Support)))
    Success= 1;
else
    Success = 0;
end
