%Finds P - symbol error probability (symbol error rate)
%for given channel over all chip sequences, all noise events, and all ambient light events.
%It can be speculated that BER=P/M.
%Core algorithm module.
%In this procedure, P is denoted as BAmb.
%
%See main_PAPM_default Test Case for description of input parameters.
function retv=simulatePAPM()
global a
global shotNoisePresented
%global ambientLightPresented
global amSAR
global amInteferenceSummationPoints
global amInterferencePeriodTi
global Amax
global L
global M
global SN
global T
global S
global tapsNumber
global lambda
global scaled_chip_length
%Numerical integration amount:
NIPoints=100
%-------------------------------------------------
%Prevent errors:
%- - - - - - - - - - - - - - - - - - - - - - - - -
tapsSimulationLimit=10000;
if tapsSimulationLimit<=tapsNumber
message='Taps Number limit exceeded.'
return;
end
if L<2
message='Incorrect value: L<2.'
return;
end
%- - - - - - - - - - - - - - - - - - - - - - - - -
%Prevent errors:
%-------------------------------------------------
%Adjust x-scale adopted in MatLab for erfc:
sqrt2m1=1.0/sqrt(2);
%First, find out number of preceding symbols:
sslots=0;
%We'l try to take enough slots to cover ISI tapsNumber:
while sslots*L<tapsNumber
sslots=sslots+1;
end
%Recalculated taps number occupied by preceding symbols:
pastTaps=sslots*L;
%Recalculate tapsNumber including primary symbol:
symTapsNumber=pastTaps+L;
b=[1:symTapsNumber];
bh=[1:symTapsNumber];
%Ambient contribution to primary symbol chips:
V=[1:L];
%Set type of this variable to "float":
overFlowProtector=1.0;
%-------------------------------------
%Generate symbols and chip sequences.
%- - - - - - - - - - - - - - - - - - -
AAmax=Amax; %Number of active levels.
%For PAM Amax+1.
overFlowProtector=1.0*L*AAmax;
if overFlowProtector>2.0E9
sprintf(' L*AAmax>2E9, ~ 32 bit limit.');
L
AAmax
return;
end
symSlotWeight=L*AAmax;
%Find out number of all combinations of preceding symbols:
symbol_events=1;
overFlowProtector=1.0;
for i=1:sslots
overFlowProtector=overFlowProtector*symSlotWeight;
if overFlowProtector>2.0E9
sprintf('symSlotWeight^preSymbolSlots>2E9, ~ 32 bit limit.');
L
AAmax
preSymbolSlots
return;
end
symbol_events=symbol_events*symSlotWeight; %NewForPAPM
end
%Now, symbol_events=L^sslots
%We don't know at the moment which limit to take:
%Take it from a cloud now:
PPMSymbolSimulationLimit=1000000;
%Protect against long calculations:
if PPMSymbolSimulationLimit<=symbol_events
sprintf('Symbol Slots Limit exceeded.');
symbol_events
PPMSymbolSimulationLimit
return;
end
%- - - - - - - - - - - - - - - - - - -
%Generate symbols and chip sequences.
%-------------------------------------
%Shortcuts:
amK=1.0/amSAR; %Parameter K-declared in [7, Wong, ...]
%%%KVK.Fix1. Please note (L-1) removed. Apparently was a mistake.
%weight_ISI_NOISE=1.0/symbol_events/L/(L-1);
weight_ISI_NOISE=1.0/symbol_events/(symSlotWeight*M); %NewForPAPM
%Prepare constants for numerical integration:
GaussNorm=1.0/sqrt(2.0*pi); %NewForPAPM
BAmb=0.0;
amInterferenceStep=amInterferencePeriodTi/amInteferenceSummationPoints;
for iXAm=0:amInteferenceSummationPoints-1
tt=amInterferenceStep*iXAm;
BB=0.0; %"BER under integration sign" by time.
%Prepare ambient contributions to current symbol:
for i=0:L-1
V(i+1)=amK*am_vi(tt+T*i, T);
end
for e=0:symbol_events-1
%-------------------------------------
%Generate symbols and chip sequences.
%- - - - - - - - - - - - - - - - - - -
%
% Symbol is encoded as couple
% (d,A) = (PostionOfPrimaryChip, AmplitudeLevelOfPrimaryChip).
% 0<=d<=L-1, 0<A<=Amax
%
% (d,A) is encoded as number sym:
% sym=d*Amax+(A-1)
%
% Hence, sym has range 0<=sym<symSlotWeight=L*Amax.
% In turn, for sequence of Lps symbols, the full number
% of possible sequences is
%
% symbol_events = Lps ^ symSlotWeight
%
% Notations in program:
% preSymbolSlots=Lps
% primaryChip=d
%
mask=symbol_events;
for slot=0:sslots-1
%sym=rem(mask,L);
sym=rem(mask,symSlotWeight); %NewForPAPM
mask=mask-sym;
mask=mask/L;
%NewForPAPM
amplitude=rem(sym,L);
primaryChip=(sym-amplitude)/L;
amplitude=amplitude+1;
for i=0:L-1
b(slot*L+i+1)=0;
%Bug discovered:
%Was:
%if sym==i
%After bug fix:
if primaryChip==i
%
b(slot*L+i+1)=amplitude; %NewForPAPM
end
end
end
%- - - - - - - - - - - - - - - - - - -
%Generate symbols and chip sequences.
%-------------------------------------
%Cycle through primary chips:
for i=0:L-1
%Fill primary symbol's chips with zeros:
for k=0:L-1
b(sslots*L+k+1)=0;
end
for ia=1:Amax
%Make i-th primary chip non-zero:
%b(sslots*L+i+1)=Amax;
b(sslots*L+i+1)=ia; %NewForPAPM
%Calculate convolved chips for primary symbol:
%from pastTaps+1 to pastTaps+L:
bh=convolve(pastTaps+1,b,bh);
Zi=lambda*bh(pastTaps+i+1);
if amInteferenceSummationPoints>1
Zi=Zi+V(i+1);
end
%NewForPAPM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
Gplus=lambda*(ia+0.5)-Zi;
Gminus=lambda*(ia-0.5)-Zi;
%Setup integration over y.
skipSummation = false;
normedGplus=0.0;
normedGminus=0.0;
%We have no processor power to integrate numerically.
%Deploy tricks ...
if ~shotNoisePresented
if Gplus<0.0 || Gminus>0.0
skipSummation=true;
end
else
normedGplus=Gplus*SN;
normedGminus=Gminus*SN;
%-- Fix3Sept14 -----------------------
%Set limits for extreme levels for PAPM:
if 1==ia
normedGminus=-100.0;
end
if Amax==ia
normedGplus=100.0;
end
%Hence, when Amax==1, then limit of integration for
%y is -oo to +oo.
%-- Fix3Sept14 -----------------------
normedGplus=min(10.0,normedGplus);
normedGminus=max(-10.0,normedGminus);
if normedGplus<-9.0
skipSummation=true;
end
if normedGminus>9.0
skipSummation=true;
end
%Now, we excluded infinite trunks ...
end
if skipSummation
continue;
end
PSuccess=0.0;
%Now, do integrate numerically ...
NIStep= (normedGplus-normedGminus)/NIPoints;
for iNIy=0:NIPoints-1
%Take values in the middle of intervals: add 0.5 to index:
y=(iNIy+0.5)*NIStep+normedGminus;
yWeight=GaussNorm*exp(-y*y*0.5)*NIStep;
%NewForPAPM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Cycle through competing chips:
PRODUCT=1.0; %NewForPAPM
for j=0:L-1
if i==j
continue;
end
Zj=lambda*bh(pastTaps+j+1);
if amInteferenceSummationPoints>1
Zj=Zj+V(j+1);
end
G=Zi-Zj;
SNRf=sqrt2m1; %SNR factor
if ~shotNoisePresented
SNRf=-1;
end
%BB=BB+erfh(G,SNRf);
%do1 make y unnormed:
%arg=y+G*SN
qq=max((1-erfh(y+G*SN,SNRf)),0.0);
PRODUCT=PRODUCT*qq; %NewForPAPM
end %Cycle through competing chips
PSuccess=PSuccess+yWeight*PRODUCT; %NewForPAPM
end %Cycle through y-integration
BB=BB+1.0-PSuccess; %NewForPAPM
end %Cycle through amplitude levels of primary chip.
end % Cycle through primary chips:
% for i=0:L-1
end % for(e
BAmb=BAmb+weight_ISI_NOISE*BB;
end % for iXAm=0 ..
BAmb=BAmb/amInteferenceSummationPoints
retv=BAmb;
end
Copyright (C) 2009 Konstantin Kirillov