home top contents previous up next

%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_PPM_default Test Case for description of input parameters.
function retv=simulatePPM()

       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
       
       %-------------------------------------------------
       %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:
       %-------------------------------------------------

       SNR2=SN/sqrt(2.0);  %Jump is built up "with" two noise events.
       %Adjust x-scale adopted in MatLab for erfc:
       SNR2=SNR2/sqrt(2.0);
       
       %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];

       %Find out number of all combination of preceding symbols:
       symbol_events=1;
	   for i=1:sslots 
		   symbol_events=symbol_events*L;
	   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
     
       %Shortcuts:
       amK=1.0/amSAR;  %Parameter K-declared in [7, Wong, ...]

       weight_ISI_NOISE=1.0/symbol_events/L/(L-1);

       
       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.
            	%- - - - - - - - - - - - - - - - - - -
            	mask=symbol_events;
          	    for slot=0:sslots-1
          	    	sym=rem(mask,L);
          	        mask=mask-sym;
                    mask=mask/L;
                    for i=0:L-1
                    	b(slot*L+i+1)=0;
                        if sym==i 
                           b(slot*L+i+1)=Amax;
                        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  
                    %Make i-th primary chip non-zero:
          	    	b(sslots*L+i+1)=Amax;
        
                    
                    %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   
                    
                    %Cycle through competing chips:
                    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=SNR2; %SNR factor
                        if ~shotNoisePresented
                           SNRf=-1;
                        end  
                        BB=BB+erfh(G,SNRf);
              	    end
                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