clear;
L=80;
N=400;
d=4;
irun=1e2;
tBK=1:15;
M=N/d;

MSE_BOMP=zeros(1,length(tBK));
MSE_BP=zeros(1,length(tBK));
MSE_OMP=zeros(1,length(tBK));
MSE_DDBOMP=zeros(1,length(tBK));
SR_BOMP=zeros(1,length(tBK));
SR_BP=zeros(1,length(tBK));
SR_OMP=zeros(1,length(tBK));
SR_DDBOMP=zeros(1,length(tBK));
for BK=tBK
    mse_bomp=zeros(1,length(irun));
    mse_bp=zeros(1,length(irun));
    mse_omp=zeros(1,length(irun));
    mse_ddbomp=zeros(1,length(irun));
    sr_bomp=zeros(1,length(irun));
    sr_bp=zeros(1,length(irun));
    sr_omp=zeros(1,length(irun));    
    sr_ddbomp=zeros(1,length(irun));
    for in0=1:irun
        xindex=d.*randperm(M,BK);
        x=zeros(N,1);
        for in1=0:d-1
            x(xindex-in1)=random('norm',0,1,length(xindex),1);
        end
        D=random('norm',0,1,L,N);
        for in1=1:N
            D(:,in1)=D(:,in1)./norm(D(:,in1));
        end
        y=D*x;
        
        rx_BOMP=BOMP(y,D,d);
        mse_bomp(in0)=sqrt(norm(rx_BOMP-x).^2/N);
        if norm(rx_BOMP-x)/norm(x)<1e-2
            sr_bomp(in0)=1;
        end

        dd=4;
        rx_DDBOMP=BOMP(y,D,d+dd);
        mse_ddbomp(in0)=sqrt(norm(rx_DDBOMP-x).^2/N);
        if norm(rx_DDBOMP-x)/norm(x)<1e-2
            sr_ddbomp(in0)=1;
        end
        
        rx_BP=BP(y,D);
        mse_bp(in0)=sqrt(norm(rx_BP-x).^2/N);
        if norm(rx_BP-x)/norm(x)<1e-2
            sr_bp(in0)=1;
        end
        
        rx_OMP=OMP(y,D);
        mse_omp(in0)=sqrt(norm(rx_OMP-x).^2/N);
        if norm(rx_OMP-x)/norm(x)<1e-2
            sr_omp(in0)=1;
        end
    end
    MSE_BOMP(BK)=mean(mse_bomp);
    MSE_BP(BK)=mean(mse_bp);
    MSE_OMP(BK)=mean(mse_omp);   
    MSE_DDBOMP(BK)=mean(mse_ddbomp);
    SR_BOMP(BK)=sum(sr_bomp)/irun;
    SR_BP(BK)=sum(sr_bp)/irun;
    SR_OMP(BK)=sum(sr_omp)/irun;   
    SR_DDBOMP(BK)=sum(sr_ddbomp)/irun;
%     figure(BK);
%     subplot(5,1,1);
%     plot(1:N,x.');
%     xlabel('x');
%     subplot(5,1,2);
%     plot(1:N,rx_BOMP.');
%     xlabel('rx BOMP');
%     subplot(5,1,3);
%     plot(1:N,rx_BP.');
%     xlabel('rx BP');
%     subplot(5,1,4);
%     plot(1:N,rx_OMP.');
%     xlabel('rx OMP');
%     subplot(5,1,5);
%     plot(1:N,rx_DDBOMP.');
%     xlabel('rx DDBOMP');
end

figure(BK+1);
plot(tBK,MSE_BOMP,tBK,MSE_BP,tBK,MSE_OMP,tBK,MSE_DDBOMP);
legend('BOMP','BP','OMP','DDBOMP');
xlabel('Sparsity K');
ylabel('MSE');
figure(BK+2);
plot(tBK,SR_BOMP,tBK,SR_BP,tBK,SR_OMP,tBK,SR_DDBOMP);
legend('BOMP','BP','OMP','BOMP(dd=4)');
xlabel('Sparsity K');
ylabel('Success Rate');