ExactBER(varargin)

function ExactBER(varargin)
%
% Generates MatLab-code expressions describing the exact theoretical
% uncoded ("raw") bit error ratio vs. Eb/N0 and Es/N0 for the transmission
% of symmetrical binary signals using Gray-coded 2^2n-QAM modulation over
% an AWGN channel with a coherent detection receiver.
%
% LastMod 2008-04-25    Paul A.M. Bun - Alcatel-Lucent Stuttgart, Germany
% Created 2008-02-04    Paul A.M. Bun - Alcatel-Lucent Stuttgart, Germany
%
%
% Use:
% ===
%
% ExactBER(n)
% ExactBER(n, 'details')
%
% n = 1: QPSK
% n = 2: QAM16
% n = 3: QAM64
% n = 4: QAM256
% n = 5: QAM1024
% n = 6: QAM4096
% n = 7: QAM16384
% n = 8: QAM65536
%     etc.
%
%
% Example:
% =======
%
% >> ExactBER(2)
%
% Theoretical AWGN uncoded BER for QAM16:
% 
% EsN0 = 4 * EbN0
%
% F = sqrt(EsN0/10)
%
% BER = (3/8) * erfc(F) + (1/4) * erfc(3*F) - (1/8) * erfc(5*F)
%
%

ShowDetailsFlag = false;

if nargin == 0 || numel(varargin{1}) > 1 || varargin{1} < 1 || ...
    mod(varargin{1}, 1)
  disp('Please type "help ExactBER" for instructions.');
  return;
else
  NrOfClasses = varargin{1};
end

if nargin == 2
  if ischar(varargin{2}) && lower(varargin{2}(1)) == 'd'
    ShowDetailsFlag = true;
  else
    disp('Unknown option.');
    disp('Please type "help ExactBER" for instructions.');
    return;
  end
end

if NrOfClasses > 1
  ModulationString = ['QAM' num2str(bitshift(1, 2 * NrOfClasses))];
else
  ModulationString = 'QPSK';
end

disp(' ');
disp(['Theoretical AWGN uncoded BER for ' ModulationString ':']);
disp(' ');
disp(['EsN0 = ' num2str(2 * NrOfClasses) ' * EbN0']);
disp(' ');

Power = 0;
for iX = 1:bitshift(1, NrOfClasses - 1)
  for iY = 1:bitshift(1, NrOfClasses - 1)
    Power = Power + (2 * iX - 1)^2 + (2 * iX - 1)^2;
  end
end
Power = Power / (bitshift(1, NrOfClasses - 1))^2;
disp(['F = sqrt(EsN0/' num2str(Power) ')']);

BasicPattern = [1; 0; 0; 1];
NrOfPositions = bitshift(1, NrOfClasses);
BitConstellation = zeros(NrOfClasses, NrOfPositions);
Terms = zeros(NrOfClasses, NrOfPositions - 1, 2);
for iClass = 1:NrOfClasses
  if ShowDetailsFlag
    disp(' ');
  end
  SpreadingFactor =  bitshift(1, NrOfClasses - iClass);
  NrOfBitChanges = 0;
  BitChange = zeros(1, NrOfPositions);
  for iPosition = 1:NrOfPositions
    BasicPatternIndex = 1 + ...
      mod(floor((iPosition - 1)/ SpreadingFactor), 4);
    BitConstellation(iClass, iPosition) = BasicPattern(BasicPatternIndex);
    if iPosition > 1
      if BitConstellation(iClass, iPosition) ~= ...
          BitConstellation(iClass, iPosition - 1)
        NrOfBitChanges = NrOfBitChanges + 1;
        BitChange(NrOfBitChanges) = 2 * (iPosition - 1) - NrOfPositions;
      end
    end
    BitChange = BitChange(1:NrOfBitChanges);
  end
  for iBit = 0:1
    BitDimension = iBit + 1;
    for iPosition = 1:NrOfPositions
      if BitConstellation(iClass, iPosition) == iBit
        Position = 2 * iPosition - 1 - NrOfPositions;
        BitChangePosition = find(Position < BitChange, 1);
        if isempty(BitChangePosition)
          BitChangePosition = NrOfBitChanges + 1;
        end
        SignValue = -1;
        for iBitChange = (BitChangePosition - 1):-1:1
          SignValue = -SignValue;
          TermPosition = ...
            bitshift(abs(Position - BitChange(iBitChange)) + 1, -1);
          Terms(iClass, TermPosition, BitDimension) = ...
            Terms(iClass, TermPosition, BitDimension) + SignValue;
        end
        SignValue = -1;
        for iBitChange = BitChangePosition:NrOfBitChanges
          SignValue = -SignValue;
          TermPosition = ...
            bitshift(abs(Position - BitChange(iBitChange)) + 1, -1);
          Terms(iClass, TermPosition, BitDimension) = ...
            Terms(iClass, TermPosition, BitDimension) + SignValue;
        end
      end
    end
  end
  if ShowDetailsFlag
    if sum(abs(Terms(iClass, :, 1) - Terms(iClass, :, 2))) > 0
      for iBit = 0:1
        disp(['Class ' num2str(iClass) ', bit = ' num2str(iBit) ':']);
        WriteLine(NrOfPositions, Terms(iClass, :, iBit + 1), 2);
      end
    end
    if NrOfClasses > 1
      disp(['Class ' num2str(iClass) ':']);
      WriteLine(2 * NrOfPositions, sum(Terms(iClass, :, :), 3), 2);
    end
  end
end
if ShowDetailsFlag
  if NrOfClasses > 1
    disp(' ');
  end
end
if ShowDetailsFlag
  disp(['Overall ' ModulationString ':']);
  WriteLine(2 * NrOfPositions * NrOfClasses, ...
    sum(sum(Terms(:, :, :), 3), 1), 2);
else
  disp(' ');
  WriteLine(2 * NrOfPositions * NrOfClasses, ...
    sum(sum(Terms(:, :, :), 3), 1), 0);
end
disp(' ');
return;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function WriteLine(Denominator, TermVector, Indentation)

Line = [repmat(' ', 1, Indentation) 'BER = '];
FirstTerm = true;
for iTerm = 1:length(TermVector)
  Term = TermVector(iTerm);
  if Term ~= 0
    if ~FirstTerm
      if Term < 0
        Line = [Line ' - '];
      else
        Line = [Line ' + '];
      end
    else
      if Term < 0
        Line = [Line '-'];
      end
    end
    FirstTerm = false;
    if abs(Term) == Denominator
      Line = [Line 'erfc('];
    else
      TermGCD = gcd(abs(Term), Denominator);
      Line = [Line '(' num2str(abs(Term)/TermGCD) '/' ...
        num2str(Denominator/TermGCD) ') * erfc('];
    end
    if iTerm > 1
      Line = [Line num2str(2 * iTerm - 1) '*'];
    end
    Line = [Line 'F)'];
  end
end
disp(Line);
return;