% function [Ps,Pb] = pb_pam_awgn (EbN0db,M,method,map);
% ------------------------------------------------------------------------------
% EINGABE:
%   EbN0db: Signal/Noise pro Infobit [db] 
%           (Skalar, Spalten- oder Zeilenvektor)
%        M: Stufigkeit der Modulation (default=2)
%           (Skalar)
%   method: Loesungsverfahren, das fuer M-PAM angewendet werden soll
%           'approx' (default)   
%             Approximation nach [Pro96,S.268,(5-2-45)]
%           'int'
%             Integration und Gewichtung der "Abschnitswahrscheinlichkeiten"
%           (string)
%      map: Art des Mappings
%           'gray' (default)   
%             Gray-Codierung
%           'nat'
%             natuerliches Mapping
%           (string)
%
% AUSGABE:
%   Pb: Bitfehlerwahrscheinlichkeit (Spaltenvektor)
%   Ps: Symbolfehlerwahrscheinlichkeit (Spaltenvektor)
%
% QUELLE: 
%    [Pro96,S.268,(5-2-45)]
%
% AUTOR: Juergen Rinas,  08.12.1998
% Erweitert um Symbolfehlerraten von Volker Kuehn, 04.01.01
% ------------------------------------------------------------------------------

function [Ps,Pb] = pb_pam_awgn (EbN0db,M,method,map)

if (sscanf(version,'%d')<5) 
  disp(['# Achtung: ',mfilename,' wurde unter Matlab 5.1 entwickelt.']); 
end;
if (nargin<2) M=2; end;             % default: M=2
if (nargin<3) method='approx'; end; % default: method='approx'
if (nargin<4) map='gray'; end;      % default: map='gray'
K=log2(M);
if (K~=fix(K))
  disp(['# Achtung(',mfilename,') M ist keine Potenz von 2.']); 
end;

EbN0=10.^(EbN0db(:)./10);
EsN0=K*EbN0;

% exakte Loesung fuer Ps
Ps=(M-1)/M*erfc(sqrt(EsN0*3/(M^2-1)));

if ((upper(method(1))=='A') & (upper(map(1))=='G'))
  % Naeherung fuer grosse EbNo
  Pb=1/K*Ps;
else
  if ((M==2) & (upper(map(1))=='G'))
    Pb=1/2*erfc(sqrt(EbN0));
  elseif ((M==4) & (upper(map(1))=='G'))
    Pb=3/8*erfc(sqrt(2/5*EbN0))...
      +1/4*erfc(sqrt(18/5*EbN0))...
      -1/8*erfc(sqrt(10*EbN0));
  elseif ((M==8) & (upper(map(1))=='G'))
    Pb=7/24*erfc(sqrt(1/7*EbN0))...
      +1/4*erfc(sqrt(9/7*EbN0))...
      -1/24*erfc(sqrt(25/7*EbN0))...
      +1/24*erfc(sqrt(81/7*EbN0))...
      -1/24*erfc(sqrt(169/7*EbN0));
  else
    if (upper(map(1))=='N')
        codetable=nattable(K);
    else
        codetable=graytable(K); 
    end
    
    % mittlere Leistung des PAM-Signals [Pro95,S.267,(5-2-40)]
    P0=1/6*(M^2-1)*(0.5^2)*2;
    
    % Berechnung der Gewichtung
    P=zeros(M-1,1);
    for APriori=1:M
      AP=codetable(APriori,:);
      for RXdiff=1:M-1
        for RX=[APriori-RXdiff APriori+RXdiff]
          Pidx=abs(APriori-RX);
          RX=max(RX,1);
          RX=min(RX,M);
          P(Pidx)=P(Pidx)+ sum(AP~=codetable(RX,:));
        end;
      end;
    end;
    
    P=P/M/K;
    
    % Umsortierung der Gewichtung, Zusammenfassung gleicher erfc-Terme
    Pw=zeros(length(P),1);
    for i=1:length(P)-1
      Pw(i)=Pw(i)+P(i);
      Pw(i+1)=Pw(i+1)-P(i);
    end;
    i=length(P);
    Pw(i)=Pw(i)+P(i);
    
    % Berechnung und Gewichtung der Wahrscheinlichkeiten
    Pb=zeros(length(EsN0),1);
    for i=1:length(Pw)
      Pb=Pb+1/2*Pw(i)*erfc(sqrt((i-0.5)^2*EsN0/P0));
    end;
  end;
end;  
  
% EOF --------------------------------------------------------------------------
