LB=500000; % Symbols transmitted
%L=LB*SF; % chips transmitted

    Z=zeros(1,LB);
    bk=zeros(1,LB);
    bk_trans1=zeros(1,LB);
    bk_trans2=zeros(1,LB);
    ray0=zeros(1,LB);   % tx0 -> rx0
    ray1=zeros(1,LB);   % tx1 -> rx0
    ray2=zeros(1,LB);   % tx0 -> rx1
    ray3=zeros(1,LB);   % tx1 -> rx1
    
    bk(1,:)=randint(1,LB,2)*2-1;
        
    noisex=randn(1,LB/2+128);
    noisey=randn(1,LB/2+128);
    fmn=0.01;
    h=fir1(255,fmn);
    noisex1=filter(h,1,noisex);
    noisey1=filter(h,1,noisey);
    ray0(1,1:2:end)=noisex1(1+128:end)+sqrt(-1)*noisey1(1+128:end);
    ray0(1,2:2:end)=ray0(1,1:2:end);
    ray0(1,:)=ray0(1,:)./std(ray0(1,:));
        
    noisex=randn(1,LB/2+128);
    noisey=randn(1,LB/2+128);
    fmn=0.01;
    h=fir1(255,fmn);
    noisex1=filter(h,1,noisex);
    noisey1=filter(h,1,noisey);
    ray1(1,1:2:end)=noisex1(1+128:end)+sqrt(-1)*noisey1(1+128:end);
    ray1(1,2:2:end)=ray1(1,1:2:end);
    ray1(1,:)=ray1(1,:)./std(ray1(1,:));
    
    noisex=randn(1,LB/2+128);
    noisey=randn(1,LB/2+128);
    fmn=0.01;
    h=fir1(255,fmn);
    noisex1=filter(h,1,noisex);
    noisey1=filter(h,1,noisey);
    ray2(1,1:2:end)=noisex1(1+128:end)+sqrt(-1)*noisey1(1+128:end);
    ray2(1,2:2:end)=ray2(1,1:2:end);
    ray2(1,:)=ray2(1,:)./std(ray2(1,:));
    
    noisex=randn(1,LB/2+128);
    noisey=randn(1,LB/2+128);
    fmn=0.01;
    h=fir1(255,fmn);
    noisex1=filter(h,1,noisex);
    noisey1=filter(h,1,noisey);
    ray3(1,1:2:end)=noisex1(1+128:end)+sqrt(-1)*noisey1(1+128:end);
    ray3(1,2:2:end)=ray3(1,1:2:end);
    ray3(1,:)=ray3(1,:)./std(ray3(1,:));
    
noise0=(randn(1,LB)+sqrt(-1)*randn(1,LB))/sqrt(1);
noise1=(randn(1,LB)+sqrt(-1)*randn(1,LB))/sqrt(1);
noise2=(randn(1,LB)+sqrt(-1)*randn(1,LB))/sqrt(1);
noise3=(randn(1,LB)+sqrt(-1)*randn(1,LB))/sqrt(1);

for k=1:30; %%%%%%%% Calculate the different BER while SNR increaseing 
    
    SNRdB(k)=k-1;
    SNR(k)=10^(-SNRdB(k)/20);
    for p=1:2:LB
        r0=ray0(1,p)*bk(1,p)+ray1(1,p)*bk(1,p+1)...
                                        +sqrt(2)*SNR(k)*noise0(1,p)/sqrt(2);
        r1=-ray0(1,p)*conj(bk(1,p+1))+ray1(1,p)*conj(bk(1,p))...
                                        +sqrt(2)*SNR(k)*noise1(1,p)/sqrt(2);
        r2=ray2(1,p)*bk(1,p)+ray3(1,p)*bk(1,p+1)...
                                        +sqrt(2)*SNR(k)*noise2(1,p)/sqrt(2);
        r3=-ray2(1,p)*conj(bk(1,p+1))+ray3(1,p)*conj(bk(1,p))...
                                        +sqrt(2)*SNR(k)*noise3(1,p)/sqrt(2);
                                    
        s0=conj(ray0(1,p))*r0+ray1(1,p)*conj(r1)...
            +conj(ray2(1,p))*r2+ray3(1,p)*conj(r3);
        
        s1=conj(ray1(1,p))*r0-ray0(1,p)*conj(r1)...
            +conj(ray3(1,p))*r2-ray2(1,p)*conj(r3);
        
        Z(1,p)=s0;
        Z(1,p+1)=s1;
    end

    bk2=zeros(1,LB);
    bk2(1,Z(1,:)>0)=1;
    bk2(1,Z(1,:)==0)=1;
    bk2(1,Z(1,:)<0)=-1;

    delayB=0;
    diffB_BPSK= (bk(1,1:end-delayB) - bk2(1,delayB+1:end))/2;
    BER_BPSK(1,k) = sum(abs(diffB_BPSK))/(length(bk)-delayB);
    disp(sprintf('bit error probability (BPSK) = %f',BER_BPSK(1,k)));
end;
%disp(sprintf('bit error probability (BPSK) = %f',BER_BPSK(U,12)));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BER of BPSK theorey
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:length(SNRdB);
    Q_BPSK(k)=0.5*erfc(1/10^(-SNRdB(k)/20));
end