function pd=dec(snrindb,PN1,PN2,PN3,PN4)
r11=1;r12=sum(PN1.*PN2)/31;r13=sum(PN1.*PN3)/31;r14=sum(PN1.*PN4)/31;
  r22=1;r21=sum(PN1.*PN2)/31;r23=sum(PN2.*PN3)/31;r24=sum(PN2.*PN4)/31;
  r33=1;r31=sum(PN1.*PN3)/31;r32=sum(PN2.*PN3)/31;r34=sum(PN3.*PN4)/31;
  r44=1;r41=sum(PN1.*PN4)/31;r42=sum(PN2.*PN4)/31;r43=sum(PN3.*PN4)/31;
  R=[r11 r12 r13 r14;r21 r22 r23 r24;r31 r32 r33 r34;r41 r42 r43 r44];
  snr=10^(snrindb/10);
%噪聲方差假設為1，信號的幅度也設為1
sgma=1;
%比特能量
Eb=2*sgma^2*snr;
%碼片能量
%E_chip=Eb/Lc;
N=10000;
numoferr=0;
%產生10000個數據比特
for i=1:N;
    %產生第一個用戶數據
    temp=rand;
    if (temp<0.5);
        data1=-1;
    else
        data1=1;
    end;
    %產生第二個用戶數據
    temp=rand;
    if (temp<0.5);
        data2=-1;
    else
        data2=1;
    end;
    %產生第三個用戶數據
    temp=rand;
    if (temp<0.5),
        data3=-1;
    else
        data3=1;
    end;
    %產生第四個用戶數據
    temp=rand;
    if (temp<0.5);
        data4=-1;
    else
        data4=1;
    end;
    b=[data1 data2 data3 data4]';
  v=[sqrt(Eb) sqrt(Eb) sqrt(Eb) sqrt(Eb)];
  w=diag(v);
    y=R*w*b+sgma*randn(1,4)';
    y1=inv(R*w)*y;
    for j=1:4;
        if y1(j)<0;
            b1(j)=-1;
        else b1(j)=1;
            if b1(j)~=b(j);
                numoferr=numoferr+1;
            end;
        end;
    end;
end;
pd=numoferr/(4*N);