home top contents previous up next

%Finds BER for given channel over all chip sequences, all noise events, and all ambient light events.
%Core algorithm module.
function retv=simulateThresholding()
        
       tapsSimulationLimit=31;

       global a       
       global shotNoisePresented
       %global ambientLightPresented
       global amSAR
       global amInteferenceSummationPoints
       global amInterferencePeriodTi
       global Amax
       global OOK_threshold
       global SN
       global b
       global bh
       global S
       global tapsNumber
       global lambda
	   global scaled_chip_length
            
	   if tapsSimulationLimit<=tapsNumber
		   message='Taps Number limit exceeded.'
           return;
       end
       %tapsLimit:
       tL=min(tapsSimulationLimit,tapsNumber);
       unitEventsCount=1;
       eventsCount=1;
       for i=1:tL
           unitEventsCount=2*unitEventsCount;
           eventsCount=eventsCount*(Amax+1);
       end

       unitEventsCount
       eventsCount

       EVENTS_MEASURE_LIMIT=1000000;
 	   if eventsCount>EVENTS_MEASURE_LIMIT
            message='Stat. events limit exceeded'
            EVENTS_MEASURE_LIMIT
            eventsCount
	        return;
 	   end
            
       teta=lambda*OOK_threshold
       amK=1.0/amSAR
   	   amInterferenceStep=amInterferencePeriodTi/amInteferenceSummationPoints;
       BAmb=0.0;
       for iXAm=0:amInteferenceSummationPoints-1

    	    t=amInterferenceStep*iXAm;
            
            B=0.0;  %=BER
            Bup=0.0;
            Bdown=0.0;
            k=tL-1; %Probe slot to count statistics.
            for e=0:unitEventsCount-1
            	%-------------------------------
            	%Generate signal of units:
            	mask=e;
            	weight=1;
          	    for slot=0:tL-1
                    ampl=rem(mask,2); %ampl=bitand(uint32(mask),uint32(1));  %http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/ref/bitshift.html&http://www.mathworks.com/access/helpdesk/help/techdoc/ref/bitand.html
          	    	if ampl>0
                      weight=weight*Amax; 
                    end
          	    	b(k-slot+1)=ampl;
                    mask=mask-ampl;
                    mask=mask/2;  %mask=bitshift(mask,-1) %http://www.mathworks.com/access/helpdesk_r13/help/techdoc/ref/uint8.html
          	    end
            	%-------------------------------
          	    for iW=0:weight-1
          	    	%----------------------------------------------
          	    	%Decompose iW and assign amplitudes to a signal  of units:
          	    	weightS=iW;
                    for j=1:tL
                        i=tL-j;
          	    		if(b(i+1)>0)
          	    			reminder=rem(weightS,Amax);
          	    			b(i+1)=reminder+1;
          	    			weightS=(weightS-reminder)/Amax;
          	    		end
                    end
          	    	%----------------------------------------------
                    bh=convolve(k+1,b,bh);
                    %bh(k+1)=bh(k+1)
                    S=lambda*bh(k+1);
                    Z=S; %no ambient light yet.
                    if amInteferenceSummationPoints>1
                       	Z=Z+amK*am_vi(t, scaled_chip_length*a);
                    end 
                
                    bk=b(k+1);         %shortcut for computer performance 
                    nUp=1.0*(lambda*bk+teta-Z);     %noise Up
           	        nDown=1.0*(Z-(lambda*bk-teta)); %noise Down

                    SNRf=SN; %SNR factor
                    if(~shotNoisePresented)
                      SNRf=-1;
                    end  

                    up=0.0;
                    down=0.0;
                    %We have three principal cases, bk=0, bk=A, bk in the middle.
                    if(0==bk)
                       up=erfh(nUp,SNRf);
                    elseif(Amax==bk)
                       down=erfh(nDown,SNRf);
                    else
                       down=erfh(nDown,SNRf);
                       up=erfh(nUp,SNRf);
           	        end
                    Bup=Bup+up;
                	B=B+up;
                    Bdown=Bdown+down;
                	B=B+down;
                end % iW    
            end %for e
            B=B/eventsCount
            Bup=Bup/eventsCount
            Bdown=Bdown/eventsCount
            BAmb=BAmb+B/amInteferenceSummationPoints;
         end %       for(int iXAm=0; iXAm<amInteferenceSummationPoints; iXAm++){
         BAmb
         retv=BAmb;
        end

Copyright (C) 2009 Konstantin Kirillov