% qamdemo.m by drb on 17-Feb-2004
% modified by drb on 14-Feb-2005
% function [y,e,pe]=qamdemo(M,Eb,N0,NumberOfOutputs);
%
% y contains the noisy PSK symbols
% e is the symbol error rate (actual)
% pe is the predicted symbol error rate
%
% M is the number of distinct QAM symbols (log2(M) is an even integer)
% Eb is the desired average energy per bit
% N0 is the desired noise PSD (it's actually N0/2, you just spec N0)
% NumberofOutputs is the the total number of correlator outputs to generate

function [y,e,pe]=qamdemo(M,Eb,N0,NumberOfOutputs);

% Set up constellation points
MM=sqrt(M);			% MM needs to be an integer
A_pam=2*[1:MM]-(MM+1);		% basic MM-PAM pattern
A=zeros(2,M);			% make space for QAM alphabet
for i=1:MM,			% loop to make QAM alphabet
  for j=1:MM,
    A(:,i+(j-1)*MM)=[A_pam(i) ; A_pam(j)];
  end
end
Es_av=mean(sum(A.^2));		% average energy per symbol
Eb_av=Es_av/log2(M);		% average energy per bit
A=A/sqrt(Eb_av)*sqrt(Eb);	% scale alphabet so the avg energy/bit = Eb

% Make symbols
index=ceil(rand(NumberOfOutputs,1)*M);	% index 1 to M
S=A(:,index);				% M-PAM symbols with correct Eb

% Make Gaussian noise with correct variance
n=sqrt(N0/2)*randn(2,NumberOfOutputs);  % uncorrelated noise

y=S+n;					% generate outputs

% Compute symbol error rate
decisionindex=zeros(1,NumberOfOutputs);
for i=1:NumberOfOutputs,
    % find minimum distance constellation point
    [junk,decisionindex(i)]=min(sum((y(:,i)*ones(1,M)-A).^2));
end
e=sum(decisionindex~=index.')/NumberOfOutputs;
pe=2*(1-1/sqrt(M))*erfc(sqrt(3*log2(M)*Eb/2/(M-1)/N0));