Basic Matching Pursuit Channel Estimation

Np=24;
Nd=128;
L=19;
TapsL=[1 10 12 19 20];
S=length(TapsL);

Index=zeros(S,1);%Tap Locations   
    c=zeros(L+1,1);
    cHatMP=zeros(L+1,1);
    hHatMP=zeros(S,1);%Taps CS Estimates
    
    h=(1/sqrt(2.*(S))).*(randn(S,1)+1i.*randn(S,1));
    c(TapsL)=h;
   
   
    
    %Channel Estimation
    np=(1/sqrt(2*SNR(ii))).*(randn(Np+L,1)+1i.*randn(Np+L,1));
    n=(1/sqrt(2*SNR(ii))).*(randn(Nd+L,1)+1i.*randn(Nd+L,1));
    
    xpCP=[xp(end-L+1:end);xp];
    
    yp=filter(c,1,xpCP)+np;
    
    yp(1:L)=[];
    
    cHatLS=inv(transpose(conj(A))*A)*transpose(conj(A))*yp;
    CHatLS=fft(cHatLS,Nd);
    WZFLS=1./CHatLS;
    WMMSELS=conj(CHatLS)./(abs(CHatLS).^2+1/SNR(ii));
    
    for ss=1:S
        max=0;
        for pp=1:L+1
            if ss>1
                set=0;
                for mm=1:ss-1
                    if Index(mm)==pp
                        set=1;
                    end
                end
                if set==0
                    Proj=(abs(transpose(conj(A(:,pp)))*yp)^2)/norm(A(:,pp))^2;
                    if Proj>max
                       max=Proj;
                       maxIndex=pp;
                    end 
                end
            else
                Proj=(abs(transpose(conj(A(:,pp)))*yp)^2)/norm(A(:,pp))^2;
                if Proj>max
                   max=Proj;
                   maxIndex=pp;
                end
            end                   
        end
        hHatMP(ss)=transpose(conj(A(:,maxIndex)))*yp/norm(A(:,maxIndex))^2;
        Index(ss)=maxIndex;
        yp=yp-((transpose(conj(A(:,maxIndex)))*yp*A(:,maxIndex))/norm(A(:,maxIndex))^2);
    end 
    
     
    cHatMP(Index)=hHatMP;