# Copyright © 2006  Francesco Potortì
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software Foundation,
## Inc., 51 Franklin Street, 5th Floor, Boston, MA 02110-1301, USA.

## -*- texinfo -*-
## @deftypefn {Function File} {} wifiper (@var{d})
## @deftypefnx {Function File} {} wifiper (@var{d}, @var{len}, @var{s}, @var{plo}, @var{arq}, @var{h}, @var{teta}, @var{f}, @var{epsr})
##
## IEEE 802.11 packet error rate
##
## Use a 2-ray propagation model in AWGN, good for unobstructed rural area.
##
##  @itemize @w
##  @item @var{d} is the distance in metres as a scalar or matrix greater than 0.
##  @item @var{len} is the length of the transmitted data in bytes as a scalar or
## matrix, not including the LLC and SNAP headers; default is 1500,
## accepted values are 0..2304.
##  @item @var{s} is the speed in Mb/s; default is 11, value is one of 1,
## 2, 5.5, 6, 9, 11, 12, 18, 24, 36, 48, 54.
##  @item @var{plo} is the path loss offset in dB with respect to a
## station pair having a frame error rate of 8e-2 at 200 m for frames
## containing 1024 bytes of data at 11 Mb/s.  Negative values indicate
## less powerful transmitter, less sensitive receiver, lower antenna
## gain, more noise or interference; default is 0.
##  @item @var{arq} is the maximum number of times a packet is transmitted;
## default is 8, value must be not less than 1.
##  @item @var{h} is  the height in metres from the ground of Tx and Rx
## stations; default is [1 1], values must be positive.
##  @item @var{teta} is the polarisation angle of the signal.  Values must
## range from 0 (horizontal polarisation) to @code{pi}/2 (vertical
## polarisation); default is vertical polarisation.
##  @item @var{f} is the channel frequency in MHz; default is 2442 MHz, which
## is channel 7 of IEEE 80211b, reasonable values are those specified
## for 802.11a and 802.11b, no checks are done.
##  @item @var{epsr} is the relative permittivity of the ground surface, assumed
## to be non-conductive, non-ferromagnetic; default is 15, reasonable
## values are between 4 and 25 for ground, 80 for water, no checks are
## done.
## @end itemize
##
## To account for slow fading having Rice factor @var{ricek} in dB and
## bandwidth @var{bw} in Hz, use a variable @var{sfplo} initialised to
## @math{10^(plo/20)}. Each time you call @var{wifiper}, extract a random
## value @var{x} from a Rice distribution having mean @math{10^(plo/20)} and
## Rice factor @math{10^(ricek/20)}, set @math{a=2*pi*bw*timed}, where @var{timed}
## is the time difference between the current call and the previous
## call, then do @math{sfplo=a*x+(1-a)*sfplo}, and use @math{20*log10(sfplo)} as
## the @var{plo} argument of @var{wifiper}.  For the standard case of
## slow fading having Rice factor equal to 20 dB and bandwidth of 1 Hz,
## you can do @math{x=10^(plo/20)*(1+randn/10)}.
##
## @strong{Reference}
## Paolo Barsocchi, Gabriele Oligeri, Francesco Potortì, "Frame error
## model in rural Wi-Fi networks</I>", proceedings of the International
## Symposium on Modeling and Optimization (Wiopt), WiNMee/WiTMeMo
## workshop, IEEE, Limassol (CY), April 2007, No. 4146.
## @end deftypefn
##
## @seealso{wifitmt}

## Author: Francesco Potortì <Potorti@isti.cnr.it>
## $Revision: 1.8 $
## Usage: p = wifiper(d, [len, s, plo, arq, h, teta, f, epsr])
## Description: Compute IEEE 802.11 packet error rate

function per = wifiper(d, len, s, plo, arq, h, teta, f, epsr)

  ## Constants
  c = 3e8;			# light speed [m/s]
  f7 = 2442;			# 802.11 frequency (channel 2) [MHz]
  llcsnap = 8;			# length of LLC plus SNAP headers [B]
  fcs = 4;			# length of 3-address MAC header plus FCS [B]
  maxdata = 2312;		# 802.11 max data length [B]
  ethdata = 1500;		# an Ethernet maximum size packet length [B]

  ## Nominal speeds [Mb/s], preamble and data overheads [B], rate gains for
  ## preamble and data relative to 11 Mb/s [dB].
  ss=[ 1     2     5.5  11    6   9   12   18   24   26    48    54  ];
  lp=[ 6     6     6     6    3   3    3    3    3    3     3     3  ];
  ld=[30    30    30    30   32  32   32   32   32   32    32    32  ];
  gp=[ 7.9   7.9   7.9   7.9  5   5    5    5    5    5     5     5  ];
  gd=[ 7.9   4.9   3.0   0    5   3.5  1.9 -0.6 -3.8 -7.1 -11.5 -12.8];
  ld -= 6;			# use 3-address MAC header
  ld += llcsnap + fcs;		# add FCS, LLC and SNAP to MAC header length
  if (false)			# set to true to use short preambles
    gp((ss==2)|(ss==5.5)|(ss==11)) -= 3; # short preamble Eb/No is 3 dB less
  endif
  sno = length(ss);		# the number of available speeds


  ## Check input parameters
  if (nargchk(1,9,nargin))
    usage("wifiper (d[, len, s, plo, arq, h, teta, f, epsr])");
  endif
  if (any(d) < 0)		# check distance
    error("d must be positive");
  endif
  if (nargin < 2)		# length excluding LLC/SNAP [B]
    len = ethdata;		# default data length [B]
  elseif (len == 0.42)		# an undocumented magic number
    len = 1024-ld(ss==11);	# data len for sensitivity benchmark in 802.11
  elseif (any(len < 0) || any(len > maxdata-llcsnap)) # check length bounds
    error("len must be in 0..%d", maxdata-llcsnap);
  endif
  if (nargin < 3)		# nominal speed [Mb/s]
    si = find(ss == 11);	# default value is 11 Mb/s
  elseif (!any(s == ss))	# check admissible values for speed
    error(["s must be one of" repmat(" %d",1,sno)], ss);
  else
    si = find(ss == s);		# speed index
  endif
  if (nargin < 4)		# link budget offset [dB]
    plo = 0;			# default value
  endif
  if (nargin < 5)		# max number of transmissions of a frame
    arq = 8;			# default value is up to 7 retransmissions
  elseif (arq < 1)		# check bounds for arq
    error("arq must not be less than 1");
  endif
  if (nargin < 6)		# heights of stations from the ground [m]
    h = [1 1];			# default value is Tx and Rx 1m above the ground
  elseif (length(h) != 2 || any(h <= 0))	# check h
    error("the two elements of h must be greater than 0");
  else
    h = reshape(h, 1, 2);	# make it a row vector
  endif
  if (nargin < 7)		# polarisation angle [rad]
    teta = pi/2;		# default value is vertical polarisation
  elseif (teta < 0 || teta > pi/2)	# check bounds for teta
    error("teta must be in 0..pi/2");
  endif
  if (nargin < 8)		# polarisation angle [rad]
    f = f7;			# default value is channel 7 of IEEE 802.11b
  endif
  if (nargin < 9)		# polarisation angle [rad]
    epsr = 15;			# default value, good for ground
  endif


  ## Compute equivalent QPSK Eb/No at the receiver
  freq = f * 1e6;		# frequency [Hz]
  lambda = c/freq;		# wave length [m]
  dstart = 1;			# reference point [m]
  K = 61.54235;			# constant for 8% 1024B 11Mb/s FER at 200m

  delta = sqrt(sum(h)^2+d.^2) - sqrt(diff(h)^2+d.^2); # path difference
  cosinc = d ./ sqrt(sum(h)^2+d.^2); # cos incidence angle
  sininc = sum(h) ./ sqrt(sum(h)^2+d.^2); # sin incidence angle
  tmp1 = sqrt(epsr-cosinc.^2);
  gammahor = (epsr*sininc-tmp1)./(epsr*sininc+tmp1); # horizontal component of gamma
  gammaver = (sininc-tmp1)./(sininc+tmp1); # vertical component of gamma
  tmp2 = exp(-j*2*pi*delta/lambda) ./ (d+delta);
  Lhor = cos(teta) * (1./d + gammahor .* tmp2);
  Lver = sin(teta) * (1./d + gammaver .* tmp2);
  L = dstart * sqrt(abs(Lhor).^2 + abs(Lver).^2);
  Ldb = 20*log10(L) + plo;	# path loss
  R = Ldb + K;			# equivalent QPSK Eb/No at the receiver


  ## Compute FER (frame error rate)
  preamblnoerrp = 1;		# assuming perfect synchronisation
  plcpnoerrp = (1 - erfc(10.^((R+gp(si))/20))/2).^(8*lp(si));
  pdunoerrp = (1 - erfc(10.^((R+gd(si))/20))/2).^(8*(ld(si)+len));
  fer = 1 - preamblnoerrp .* plcpnoerrp .* pdunoerrp;


  ## Compute PER (packet error rate)
  ## WARNING: this works only if FER is deterministic, else FER should
  ##          be recomputed arq times
  per = fer.^arq;

endfunction
