pe4.m

以cos（2*pi*k*t/N）信号空间

function [p]=Pe4(N,snr_in_dB)
M=4;
k=sqrt(M);

Rb=1;
Tb=1/Rb;        %比特周期
T=k*Tb;         %符号周期  
Es=T/2;         %符号能量归一划到1
Eb=Es/k;        %比特能量
fs=5000*Rb;
ts=1/fs;
t=0:ts:T;

of=1/M;
fc0=0;
fc1=fc0+of;
fc2=fc0+2*of;
fc3=fc0+3*of;

s0=cos(2*pi*t*fc0);
s1=cos(2*pi*t*fc1);
s2=cos(2*pi*t*fc2);
s3=cos(2*pi*t*fc3);

Q0=cos(2*pi*fc0*t);
Q1=cos(2*pi*fc1*t);
Q2=cos(2*pi*fc2*t);
Q3=cos(2*pi*fc3*t);

SNR=exp(snr_in_dB*log(10)/10);
sgma=sqrt(Es^2/(4*SNR));
for i=1:N,
    temp=rand;
    if(temp<0.25),
        dsource1(i)=0;
        dsource2(i)=0;
    elseif(temp<0.5),
        dsource1(i)=0;
        dsource2(i)=1;
    elseif(temp<0.75)
        dsource1(i)=1;
        dsource2(i)=0;
    else
        dsource1(i)=1;
        dsource2(i)=1;
    end
end;
numoferr=0;
for i=1:N,
    % matched filter outputs
    if ((dsource1(i)==0)&(dsource2(i)==0)),
        r=s0+gngauss(sgma);
    elseif ((dsource1(i)==0)&(dsource2(i)==1)),
        r=s1+gngauss(sgma);
    elseif ((dsource1(i)==1)&(dsource2(i)==0)),
        r=s2+gngauss(sgma);
    else
        r=s3+gngauss(sgma);
    end;
    r0=r*(Q0');
    r1=r*(Q1');
    r2=r*(Q2');
    r3=r*(Q3');
max_r=max([r0 r1 r2 r3]);
    if (r0==max_r),
        decis1=0;
        decis2=0;
    elseif (r1==max_r),
        decis1=0;
        decis2=1;
    elseif (r2==max_r),
        decis1=1;
        decis2=0;
    else
        decis1=1;
        decis2=1;
    end;
 if (decis1~=dsource1(i)),
        numoferr=numoferr+1;
    end;
    if (decis2~=dsource2(i)),
        numoferr=numoferr+1;
    end;
end;
 
p=numoferr/(2*N);