function [rec] = medrec(Y,A,Thi,beta,ite,epsilon)
%This function implements a robust reconstruction algorithm for compressive
%sensing:
%
%INPUTs:
%Y --> measuremnt vector
%A --> Measurement matrix or Holographic dictionary 
%Thi --> Initial Threshsol
%Beta --> Tuning parameter 
%ite --> Number of iterations
%Epsilon -->Target residual energy
%
%OUTPUTS
%rec --> Reconstruction sparse signal

    [M,N]=size(A);
    Aa =abs(A);
        Ai = 1./A;
        Th = Thi;
        rec = zeros(N,1);
        Y2 = zeros(M,1);
        Yn = norm(Y);
        R=Y;
        km=0;
        while ((norm(Y2-Y,2) > epsilon) && (km <=ite) )
            for k = 1:N
                ynewa = (Y - Y2 + A(:,k)*rec(k));
                [reca] = wmediana(ynewa.*Ai(:,k),Aa(:,k),0);
                skt = sum((abs(ynewa))- abs(ynewa-A(:,k)*reca));
                rectemp=rec(k);
                rec(k) = (skt>Th)*reca;
                Y2=Y2+(rec(k)-rectemp)*A(:,k);
            end
            Th = Th*beta;
            km=km+1;
        end
        
function [y]=wmediana(X,W,a)
% Weighted median with negative weights
% The firt input argument is the data, the second argument is the weigths
% a is a flag that indicates the type of filtering operation
% a = 0  lowpass filtering operation
% a = 1  bandpass or highpass filtering operation.

if (nargin==2) a = 0;end;

% Assigning signs of weights to samples

x=reshape(X,1,size(X,1)*size(X,2));
w=reshape(W,1,size(W,1)*size(W ,2));
sw=sign(w);     % signs of weights
sw1=sw;
w=abs(w);       % non-negative weights
x=x.*(sw+~sw);  % assigning signs to samples


[z,ind]=sort(x,'descend');
%z=fliplr(z);
%ind=fliplr(ind);
tau = 0.5*sum(w);
i=0;
sw=0;
while sw<tau;
        i=i+1;
        sw=sw+w(ind(i));
end;


if (a==0)
       y = z(i);
else
    y=0.5*z(i)+0.5*z(i+1);
end;