%file qpsk_berest.m
function [peideal,pesystem]=qpsk_berest(xx,yy,ebn0db,eb,tb,nbw)
%ebn0db is an array of Eb/N0 values in db (specified at the receiver
%input);tb is the bit duration and nbw is the noise BW
%xx is the reference (ideal)input; yy is the distorted output;
%
[n1 n2]=size(xx); nx=n1*n2;
[n3 n4]=size(yy); ny=n3*n4;
[n5 n6]=size(ebn0db); neb=n5*n6;
%
%for comparision purposes,set the noise BW of the ideal
%receive(integrate and dump) to be equal to rs/2.
%
nbwideal=1/(2*tb*2);
for m=1:neb
    peideal(m)=0.0;pesystem(m)=0.0;%initialize
    %
    %find n0 and the variance of the noise
    %
    string1=['Eb/N0=',num2str(ebn0db(m))];
    disp(string1)                       %track execution
    ebn0(m)=10^(ebn0db(m)/10); %dB to linear
    n0=eb/ebn0(m);                       %noise power
    sigma=sqrt(n0*nbwideal*2);  %variance
    sigma1=sqrt(n0*nbwideal*2); %variance of ideal
    %
    % mutiply the input constellation/signal by a scale facotr so that
    % input to receive filter have the same ave power
    %a=sqrt(2*eb/(2*tb))
    %
    b=sqrt(2*eb/tb)/sqrt(sum(abs(xx).^2)/nx);
    for n=1:nx
        theta=angle(xx(n));
        if (theta<0)
            theta=theta+2*pi;
        end
        %
        %rotate x and y to the first quadrant and compute BER
        %
        xxx(n)=b*xx(n)*exp(-i*(theta-(pi/4)));
        yyy(n)=yy(n)*exp(-i*(theta-(pi/4)));
        d1=real(xxx(n));d2=imag(xxx(n)); %reference
        d3=real(yyy(n)); d4=imag(yyy(n)); %system
        pe1=q(d1/sigma1)+q(d2/sigma1); %reference
        pe2=q(d3/sigma)+q(d4/sigma); %system
        peideal(m)=peideal(m)+pe1;  %SER of reference
        pesystem(m)=pesystem(m)+pe2;
    end
end
peideal=(1/2)*peideal./nx;    %convert BER
pesystem=(1/2)*pesystem/nx;