Source Code
Infrared Signal Coding Schemes
home top contents previous up next

	import java.awt.Color;
	import java.awt.Dimension;
    import java.awt.Image;

	abstract public class MultiPathConvolution extends Content{
		public MultiPathConvolution(Dimension dim, Image img, String problemName){
			super(dim,img,problemName);
	    }

		
 	    ////////////////////////////////////////////////////
		//  AFFECTS CHILD CLASSES ONLY WHEN NOT OVERRIDDEN.
		//==================================================
 	    protected void spawnParametersAfterUserApplied(){
			 spanParsInMPLevel();
			 displySpawnedParsInMPLevel();
		     simulationBlocked=false;
			 simulateThresholding();
		}
 	    ////////////////////////////////////////////////////

		
		
		//=============================
		// Sub Model
		//-----------------------------
        int subModel; //Index of following subModels:
        String[] subModelTitle={ "OOK and PAM", "DAPPM", "PPM", "PAPM"};
        int shotNoisePresented=1;
		int amInteferenceSummationPoints=2;
		//=============================
		
		//=============================
		// Default parameters
        // See description at line:
        // User Input Prompts: 
		//-----------------------------
		double ceilingHeight=3.5;
		int Amax=3;
	    int L=1;
		double Rb=100000000;
		double SNR=7; 
		double amSAR=1.0;
		double amInterferencePeriodTi=25.0E-6;
		//=============================

		//=============================
		//Derivative parameters:
        // See description at line:
        // To display this strings in graph background: 
		//-----------------------------
		double a;
	    //SNR not in dB form:
	    double SN;	    
		double T;
	    double avLength;
		double aphabetCount;
		double M;
		double bitsPerChip;
		double scaled_chip_length;
	    int tapsNumber;
	    double[] beta;  //Discretized h.
	    //min non-zero Intensity/average Intensity:
	    double lambda;
	    
	    //Submodel specific:
	    double OOK_threshold=0.5;
		//=============================

	    //Simulation control:
		protected boolean simulationBlocked;
	    
	    //Intermediate Results:
	    //double B;  //BER
 	    double BAmb=0.0; //BER over Ti - Ambient interference period.	    
 	    
 	    
		/////////////////////////////////////////////////////////////////
		// Assign Parameters
		//================================================================

 	    
		//================================================================
 	    //Can be used in child classes directly.
		//If not used in child classes, then do not affect them.
 	    //----------------------------------------------------------------
		final protected void initPreconstructorInMPLevel(){
			modelDescription=new String[200];
			modelDescription[0]="Program:  Stub for MultiPath Indoor Dispersion.";

			tapsSimulationLimit=31;

			amPrepare();
			
			dmnStartF=0.0;
			dmnRangeF=8.0;
			dmnRangeFDown=1.0;

			dmnStartX=0.0;
			dmnRangeX=3.0;
			grPoints=400;
			grStartF=600;

			functionCOUNT=7;
	    	functionColor=new Color[]{
	    			new Color(255,0,0),
	    			new Color(220,100,0),
	    			new Color(0,0,200),
	    			new Color(150,0,150),
	    			new Color(0,150,0),
	    			new Color(0,255,0),
	    			new Color(100,100,100,50),
	    			new Color(200,0,0,50),
	    			new Color(222,0,200,50)
	    	};
	        functionTitle=new String[]{
	        	"h",
	        	"beta",
	        	"Origin",
	        	"Loss",
	        	"Origin",
	        	"Gain",
	        	"FirstLossThreshold",
	        	"Ambient V",
	        	"Ambient V over chip"
	        };
		}
		
		//To be explicitly used in child classes:
	    final protected void spanParsInMPLevel(){
	    	//con("Spawning parameters InMPLevel...");
			a=2.*ceilingHeight/300000000.0; 
			setSubModel();
			bitsPerChip=M/avLength;
			T=bitsPerChip/Rb;
			scaled_chip_length=T/a;
			
			
	        //-------------------------------------------------------
	        ///stimation of size of sequence beta:
	        //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
	        double accuracyEps=1.0e-3;
	        double hThresholdTs= Math.exp(  Math.log(1.0/accuracyEps)*(1.0/6.0) ) - 1;
	        if(hThresholdTs<1.0)hThresholdTs=1.0;
	        double wTapsNumber=hThresholdTs/scaled_chip_length;
	        tapsNumber = ((int)Math.floor(wTapsNumber)) + 1; //1 is taken for safety.
	        double wTruncationParameterForCuriosity=wTapsNumber / Math.floor(a/T);
	        con("wTruncationParameterForCuriosity"+wTruncationParameterForCuriosity);
	        //Was:
	        //tapsNumber= ((int)Math.floor(a/T*1.25)) + 1;
	        //-------------------------------------------------------

			
			
			
			
		    //Convert SN from dB to numbers:
		    SN=Math.exp( SNR/10.0*Math.log(10.0) );	    
		    setSequenceMemoryInMPLevel(); 
		    assignBetaKToArray();
	    }
	    
		//To display this strings in graph background: 
	    final protected void displySpawnedParsInMPLevel(){
		    //con("subModelTile[0]="+subModelTitle[0]+" " + subModel);
			modelDescription[1]="Sub Model = "+subModelTitle[subModel];
			modelDescription[4]="A="+Amax;
			modelDescription[5]="L="+L;
			modelDescription[6]="SNR="+SNR+",db     SNR="+SN;
			
			modelDescription[7]="Symbols in Alphabet="+aphabetCount;
			modelDescription[8]="M="+M;
			modelDescription[9]="BitRate, Rb="+Rb+", bits/second.";
			
			modelDescription[10]="Average Length="+avLength;
			modelDescription[11]="Bits per chip="+bitsPerChip;
			modelDescription[12]="Chip Length="+T+", seconds.";
			modelDescription[13]="Room Height="+ceilingHeight+", meters.";
			modelDescription[14]="Dispersion Length Scale, a="+a+", seconds.";
			modelDescription[15]="Lambda="+lambda;
			
			modelDescription[16]="tapsNumber="+tapsNumber;
	    }
		//----------------------------------------------------------------
 	    //Can be used in child classes directly.
		//================================================================

	    
		//================================================================
 	    //Can not be used in child classes directly.
 	    //----------------------------------------------------------------
	    private void setSequenceMemoryInMPLevel(){
	    	    //con(" Memory");
	    	    int len=tapsNumber+L+1;
				beta=new double[tapsNumber];
	    	    b=new int[tapsNumber+L];
			    result=new double[len];

			    //----------------------------------------------------------
			    //Set memory for signals which demonstrate ISI loss or gain:
			    inputChipLost=new int[len];
			    inputChipGained=new int[len];
			    resultChipLost=new double[len];
			    resultChipGained=new double[len];
		    	for(int i=0; i<result.length; i++){
		    		inputChipLost[i]=0;
		    		inputChipGained[i]=0;
		    	}
			    //----------------------------------------------------------
	    }
	    private void setSubModel(){
	    	int i=L;
	    	switch(subModel){
	    	case 0:    //PAM
	    		       avLength=L;	         
                       aphabetCount=Amax+1;
                       while(--i>0){
           	               aphabetCount*=Amax+1;
                       }
                       M=Math.log(aphabetCount)/Math.log(2.0);
                       lambda=2.0/Amax;
                       break;
	    	case 1:    avLength=(1+L)/2; //In Chips.
	                   aphabetCount=Amax*(1+L)*L/2;
	                   M=Math.log(aphabetCount)/Math.log(2.0);
	                   lambda=1.0/((Amax+1.0)/2/avLength);
	                   break;
   	        case 2:    //PPM
	    	           avLength=L; //In Chips.
                       aphabetCount=L;
                       M=Math.log(aphabetCount)/Math.log(2.0);
                       lambda=L;
                       break;
	        case 3:    //PAPM
	    	           avLength=L; //In Chips.
                       aphabetCount=L*Amax;
                       M=Math.log(aphabetCount)/Math.log(2.0);
                       lambda=2.0*L/(Amax+1);
                       con("Amax="+Amax+" L="+L+" lambda="+lambda);
                       break;
    	    }
	    }
 	    //----------------------------------------------------------------
 	    //Can not be used in child classes directly.
		//================================================================
		//Assign Parameters
		//////////////////////////////////////////////////////////////////
	    

	    
	    
	    
	    
		
		//==========================================
		//  Algorithms
		//------------------------------------------
		protected double hh(double q){
			double power=(q+1);
			double p=power*power;
			double p4=p*p;
			p*=p4; //6
			p*=power; //7
			return 1/p*6.0;
		}
		//Returns area under h-function over interval [k*step,k*step+step]:
        //Input: step = scaled chip length.
		private double betaPortion(double k, double step){
			double t=k*step;
			double power=t+1;
			double p=power*power;
			double p4=p*p;
			double value1=1/(p*p4); //6
			t+=step;
			power=t+1;
			p=power*power;
			p4=p*p;
			return value1-1/(p*p4);
		}
		private void assignBetaKToArray(){
			for(int k=0; k<tapsNumber; k++){
				beta[k]=betaPortion((double)k, scaled_chip_length);
				//con("beta="+beta[k]);
			}
		}
		//Calculate all result chips:
	    protected void convolve(int[] inputPul,double[] outputPul, boolean drawNumbers){
	    	convolve(0, inputPul,outputPul, drawNumbers);
	    }
        //Calculate only from start result chip:
	    protected void convolve(int start, int[] inputPul,double[] outputPul, boolean drawNumbers){
	    	int nOut=outputPul.length;
	    	int nIn=inputPul.length;
	    	for(int i=start; i<nOut; i++){
	    		double s=0;
	    		for(int j=0; j<tapsNumber; j++){
	    			int tail=i-j;
	    			if(tail<0 || tail>=nIn)break;
	    			s+=inputPul[tail]*beta[j];
	    		}
	    		outputPul[i]=s;
	    		if(drawNumbers) {
	    			int drawIndex=25+i;
	    			if(drawIndex<modelDescription.length){
                       String ss=" conv "+i+" "+s+" beta="+beta[Math.min(beta.length-1,i)];
	    			   if(i<nIn)ss+=" input="+inputPul[i];
	    			   modelDescription[drawIndex]=ss;	    			
	    			}
	    		}
	    	}
	    }
	    
		//------------------------------------------
		//  Algorithms
		//==========================================
	    

	    
	    
	    
	//==========================================
	//  Simulation
	//------------------------------------------
	protected int[] b;
	//do1 replace with bh:
	protected double[] result;
	protected int tapsSimulationLimit;
	
	protected void simulateThresholding(){
    	
  	   //Demo:
  	   //Clear display:
       int i;
  	   for(i=20; i<modelDescription.length; i++){
  	    	modelDescription[i]=null;
  	   }
       	
       if(simulationBlocked)return;
	   if(subModel>0){
			  modelDescription[20]="BER for sub Model " + subModelTitle[subModel]+" is not implemented.";
			  simulationBlocked=true;
              return;
	   }
	   if(tapsSimulationLimit<=tapsNumber){
				  modelDescription[20]="BER for sub Model " + subModelTitle[subModel]+". Taps Number exceeded limit.";
				  //Set B to 1.0 to attract attention:
				  BAmb=1.0;
				  simulationBlocked=true;
				  return;
	   }

	   //PAM and OOK:
	   int unitEventsCount=2;
       int eventsCount=Amax+1;
       int tL=Math.min(tapsSimulationLimit,tapsNumber);
       i=tL;
       while(--i>0){
            	unitEventsCount*=2;
            	eventsCount*=Amax+1;
       };

	   modelDescription[22]="Chip sequences count " + eventsCount +
		         " Unit sequences count=" + unitEventsCount;
       int EVENTS_MEASURE_LIMIT=1000000;
 	   if(eventsCount>EVENTS_MEASURE_LIMIT){
	  			modelDescription[20]="Chip sequences count " + eventsCount +
	  			     " > " + EVENTS_MEASURE_LIMIT + " which seems too long to process.";
                return;
 	   }
 		    
 	   BAmb=0.0;
 	   double amInterferenceStep=amInterferencePeriodTi/amInteferenceSummationPoints;
       for(int iXAm=0; iXAm<amInteferenceSummationPoints; iXAm++){

    	    double t=amInterferenceStep*iXAm;
            double BB=0.0;
            double Bup=0;
            double Bdown=0;
            double BupISI=0;
            double BdownISI=0;
            firstLossThreshold=-1;
            double firstGainThreshold=-1;
            
            //Shortcuts:
            double teta=lambda*OOK_threshold;
            double amK=1.0/amSAR;
            //con("teta="+teta);

            boolean HeavisideMode=0==shotNoisePresented;
            
            int k=tL-1; //Probe slot to count statistics.
            for(int e=0; e<unitEventsCount; e++){
            	
            	//-------------------------------
            	//Generate signal of units:
            	int mask=e;
            	int weight=1;
            	int count=0;
          	    for(int slot=0; slot<tL; slot++){
          	    	int ampl=mask&1;
          	    	if(ampl>0){ weight*=Amax; count++; } 
          	    	b[k-slot]=ampl;
          	        mask>>>=1;
          	    }
            	//-------------------------------

          	    for(int iW=0; iW<weight; iW++){

          	    	//----------------------------------------------
          	    	//Decompose iW and assign amplitudes to a signal 
          	    	//of units:
          	    	i=tL;
          	    	int weightS=iW;
          	    	while(--i>-1){
          	    		if(b[i]>0){
          	    			int reminder=weightS%Amax;
          	    			b[i]=reminder+1;
          	    			weightS=(weightS-reminder)/Amax;
          	    		}
          	    	}
          	    	//----------------------------------------------
          	    	
          	    
                    convolve(k,b,result,false);
                    //con("e="+e);
                	//con("iW="+iW);
                    //con("bh="+result[k]);
                    double S=lambda*result[k];
            	    //con("S="+S);
                    double Z=S; //no ambient light yet.
                    if(amInteferenceSummationPoints>1){
                	    Z=Z+amK*am_vi(t, scaled_chip_length*a);
                    }
                
                    int bk=b[k];
           	        double nUp=lambda*bk+teta-Z;     //noise Up
           	        double nDown=Z-(lambda*bk-teta); //noise Down

           	        double up=0;
           	        double down=0;

                    //We have three principal cases, bk=0, bk=A, bk in the middle.
            	    if(0==bk){
            		    up=UFun.erfh(SN*nUp,HeavisideMode);
                    }else if(Amax==bk){
            		    down=UFun.erfh(SN*nDown,HeavisideMode);
                    }else{
            		   up=UFun.erfh(SN*nUp,HeavisideMode);
            		   down=UFun.erfh(SN*nDown,HeavisideMode);
           	        }

               if(0==subModel){
                	if(S>lambda*bk+teta){
                	   if(0.0==BupISI){
                		   convolve(b,result,false);
   			               copyArray(b,inputChipGained);
 			               copyArray(result,resultChipGained);
 		                   firstGainThreshold=bk+teta/lambda;
                	   }
               	       BupISI++;
        		    }
         		    if(S<lambda*bk-teta ){
                 	   if(0.0==BdownISI){
                 		   convolve(b,result,true);
       			           copyArray(b,inputChipLost);
        			       copyArray(result,resultChipLost);
 		                   firstLossThreshold=bk-teta/lambda;
        		        }
         		    BdownISI++;
         		    }
                }    
        		Bup+=up;
        		Bdown+=down;
        		BB+=up;
            	BB+=down;
            	//con(" B="+BB+" nUp="+nUp+" up="+up+" Bup="+Bup+" nDown="+nDown+" down="+down+" Bdown="+Bdown);
          	  } //for(int iW=0; iW<weight; iW++){
            }// e 
            //con("events Count="+eventsCount); 	
            BB/=eventsCount;
            Bup/=eventsCount;
            Bdown/=eventsCount;
            BdownISI/=eventsCount;
            BupISI/=eventsCount;
 
        	//con(" B="+BB+" Bup="+Bup+" Bdown="+Bdown);

            //-------------------------------------
            //Output:
        	modelDescription[21]="";
            modelDescription[22]="";
            if(1==amInteferenceSummationPoints){
       	        modelDescription[20]="BER, BERup, BERdown = "+BB+",   "+Bup+",    "+Bdown;
                if(0==subModel){
                	modelDescription[21]="BERupISI, BERdownISI = "+BupISI+",    "+BdownISI;
                    modelDescription[22]="firstLossThreshold, firstGainThresho = "+
                                         firstLossThreshold+", "+firstGainThreshold;
                }
            }
            //-------------------------------------

            
            BAmb=BAmb+BB/amInteferenceSummationPoints;
       } // for(int iXAm=0; iXAm<amInteferenceSummationPoints; iXAm++){
       if(1<amInteferenceSummationPoints)modelDescription[20]="BER "+BAmb;
	}
	//------------------------------------------
	//  Simulation
	//==========================================

    
    

		

		//==========================================
		//  Draw Graphs
		//------------------------------------------
        protected double drawSignal(double t, double[] stepSignal){
			int k=(int)(t/scaled_chip_length);
			if(k>=stepSignal.length)return 0;
			return (double) stepSignal[k];
        }
        //Integer version of the method:
        protected double drawSignal(double t, int[] stepSignal){
			int k=(int)(t/scaled_chip_length);
			if(k>=stepSignal.length)return 0;
			return (double) stepSignal[k];
        }
		
	    abstract protected double functionSwitch(int fIx, double x);
		//------------------------------------------
		//  Draw Graphs
		//==========================================
	
	
	
    //==========================================
	//  Ambient light
	//------------------------------------------

//%
//%Credit: Data has been used from the following source:
//% Adriano J.C. Moreira, Rui T. Valadas and A.M. de Oliveira Duarte.
//% Optical interference produced by artificial light.
//% Wireless Networks 3 (1997) 131–140
//%


    protected double[] am_b;
    protected double[] am_c;
    protected double[] am_zeta;
    protected double[] am_fi;
    protected double[] am_d;
    protected double[] am_teta;
    protected double am_A1_reciprocal;
    protected double am_A2_reciprocal;
    //Fundamental frequency of high frequency component in (7), Hz:
    protected double am_fh; 
    static final double PI2=2*Math.PI;
	protected void amPrepare(){
		
	    am_b=new double[21];
	    am_c=new double[21];
	    am_zeta=new double[21];
	    am_fi=new double[21];
	    am_d=new double[12];
	    am_teta=new double[12];
	    am_A1_reciprocal=1.0/5.9;
	    am_A2_reciprocal=1.0/2.1;
	    //Fundamental frequency of high frequency component in (7), Hz:
	    am_fh=37.5E3; 

		
		
		double log10 = Math.log(10);
		for(int i=1; i<21; i++){
			am_b[i]=Math.exp(  log10*(  -13.1*Math.log(100*i-50) +27.1  )/20       );
			am_c[i]=Math.exp(  log10*(  -20.8*Math.log(100*i) + 92.4    )/20       );
		}
		//[7] Table I:
		double[] amAux=new double[]{
				1, 4.65, 0.00, 11, 1.26, 6.00,
				2, 2.86, 0.08, 12, 1.29, 6.17,
				3, 5.43, 6.00, 13, 1.28, 5.69,
				4, 3.90, 5.31, 14, 0.63, 5.37,
				5, 2.00, 2.27, 15, 6.06, 4.00,
				6, 5.98, 5.70, 16, 5.49, 3.69,
				7, 2.38, 2.07, 17, 4.45, 1.86,
				8, 4.35, 3.44, 18, 3.24, 1.38,
				9, 5.87, 5.01, 19, 2.07, 5.91,
			   10, 0.70, 6.01, 20, 0.87, 4.88
		};
		for(int ii=0; ii<10; ii++){
			for(int j=0; j<2; j++){
				int pos=ii*6+j*3;
				int i=ii+10*j+1;
				am_zeta[i]=amAux[pos+1];
			    am_fi[i]=amAux[pos+2];
			}
		}
	    //Check our selves:
		for(int i=1; i<=10; i++){
			int j=i+10;
			//con(""+i+"  "+am_zeta[i]+"  "+am_fi[i]+"  "+j+"  "+am_zeta[j]+"  "+am_fi[j]);
		}
		
		//[7] Table II. Amplitude and phase parameters for high-frequency components.
        double[] amAux2=new double[]{
	       0, -22.22, 5.09, 6, -39.30, 3.55,
	       1,   0.00, 0.00, 7, -42.70, 4.15,
	       2, -11.50, 2.37, 8, -46.40, 1.64,
	       3, -30.00, 5.86, 9, -48.10, 4.51,
	       4, -33.90, 2.04, 10, -53.10, 3.55,
	       5, -35.30, 2.75, 11, -54.90, 1.78
	    };
	    for(int ii=0; ii<6; ii++){
		    for(int j=0; j<2; j++){
			    int pos=ii*6+j*3;
			    int i=ii+6*j;
			    am_d[i]=amAux2[pos+1];
		        am_teta[i]=amAux2[pos+2];
			}
		}  
	    //Check our selves:
		for(int i=0; i<=5; i++){
			int j=i+6;
			//con(""+i+"  "+am_d[i]+"  "+am_teta[i]+"  "+j+"  "+am_d[j]+"  "+am_teta[j]);
		}
	    
	}
    
	//Here we are programming formula [7 (7)]:
	//Factors R and Pm are not included into m(t) as in formula (7):
	protected double am_mi(double t){
		//Calculate member 1 in (7), RPm.
		double sum1=1.0;

		//Calculate member 2:
		double sum2=0.0;
		for(int i=1; i<=20; i++){
			sum2 += am_b[i]*Math.cos(PI2*(100*i-50)*t+am_zeta[i])+
			        am_c[i]*Math.cos(PI2*(100*i*t+am_fi[i]));
		}
		sum2*=am_A1_reciprocal;
		//Calculate member 3:
		
		double sum3=am_d[0]*Math.cos(PI2*am_fh*t+am_teta[0]);
		for(int i=1; i<=11; i++){
			sum3 += am_d[i]*Math.cos(PI2*2*i*am_fh*t+am_teta[i]);
		}
		sum3*=am_A2_reciprocal;
		
		double sum=sum1+sum2+sum3;
		//UTil.con("amV="+sum);
		
		return sum;
	}
    //Calculate ambient light average energy contribution over time t2-t1 starting
	//from time t1. 
	//Returns: 1/dt*Integral_over amV by dt.
	//Purpose: calculate noise applied to one chip:
	protected double am_vi(double t1, double dt){
		
		double a1=0;
		double delta=0;
		double alpha=0;
		
		//Calculate member 1 in (7), RPm.
		double sum1=1.0*dt;

		//Calculate member 2:
		double sum2=0.0;
		for(int i=1; i<=20; i++){
		   //Calculate member 2:
	       double a0=PI2*100*i;
	       a1=a0-PI2*50;
	       double delta0=a0*dt*0.5;
	       double delta1=a1*dt*0.5;
	       double alpha1=a1*t1+am_zeta[i]+delta1;
	       double alpha0=a0*t1+am_fi[i]+delta0;
		   sum2 = sum2+am_b[i]*2.0*(  
	                   Math.sin(delta1)*Math.cos(alpha1)/a1+
	             	   Math.sin(delta0)*Math.cos(alpha0)/a0
				            );
        }
	  
  	    sum2=sum2*am_A1_reciprocal;
	    //Calculate member 3:
		
		a1=PI2*am_fh;
		delta=a1*dt*0.5;
		alpha=a1*t1+am_teta[0]+delta;
		double sum3=am_d[0]*2.0*Math.sin(delta)*Math.cos(alpha)/a1;
		for(int i=1; i<=11; i++){
			a1=PI2*2*i*am_fh;
			alpha=a1*t1+am_teta[i]+delta;
			sum3 += am_d[i]*Math.sin(delta)*Math.cos(alpha)/a1;
		}
		sum3*=am_A2_reciprocal;
		
		double sum=(sum1+sum2+sum3)/dt;
		//UTil.con("amV="+sum);
		
		return sum;
	}
	//------------------------------------------
	//  Ambient light
    //==========================================

	
	//==========================================
	//  Auxiliary procedures
	//------------------------------------------
	protected void copyArray(int[] f, int[] t){
		int lim=Math.min(f.length,t.length);
		for(int i=0; i<lim; i++){
			t[i]=f[i];
		}
	}
	protected void copyArray(double[] f, double[] t){
		int lim=Math.min(f.length,t.length);
		for(int i=0; i<lim; i++){
			t[i]=f[i];
		}
	}
	//------------------------------------------
	//  Auxiliary procedures
	//==========================================

	//==========================================
	// Demonstration parameters.
	// To store signals which show failure
	// detection of a chip due ISI:
	// Used: in simulateThreshold() method.
	//------------------------------------------
	protected double firstLossThreshold;
	protected int[] inputChipLost;
	protected int[] inputChipGained;
	protected double[] resultChipLost;
	protected double[] resultChipGained;
	//==========================================

	
}

Copyright (C) 2009 Konstantin Kirillov