/* compile with "cc bbsim.c -o bbsim -lm" */

/* this is a 2nd-order BB loop simulator program */
/* pipe stdout to autoplot (ap) to graph */
/* Rick Walker (walker@opus.hpl.hp.com), 1998 */

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <unistd.h>

/* for prbs.c */
#define PLEN  23             /* shift register length */
#define TAP   18             /* xor tap location */

double gauss();

#define NUM_CYCLES 40000
#define MAX_LATENCY 128

int main(argc,argv)
int argc;
char **argv;
{
    int c;
    int i;
    extern char *optarg;
    extern int optind;
    extern int opterr;
    char *progname;

    int cycle;
    time_t time_val;
    double vco_phase=0.0;
    double rms_sum_vco=0.0;
    double data_phase=0.0;
    double o_data_phase=0.0;  
    double rms_sum_data=0.0;
    int dir;

    int pd_output;	       /* intermediate variable for phase detector */
    int pd_delay[MAX_LATENCY]; /* array for introducing latency into p.d. */
    int pd_count;              /* counter for p.d latency array */


    double loop_filter=0.0;
    double amp=0.0;
    double period=100.0;

    double psi=1.0;	      /* (2*BETA*TAU)/tupdate */
    double jitter=0.0;	      /* data jitter temp variable */
    double djitter=1.0;	      /* data rms input jitter */
    double vjitter=0.0;	      /* vco rms jitter */
    double vcojitter=0.0;     /* vco jitter running integral */
    double bangbang=1.0;      /* vco phase advance/retard per frame */
    double tdensity=1.0;      /* transition probability */
    int transition=0;		

    double npoints=4000.0;    /* number of points calculated */
    double offset=0.0;        /* initial phase offset */
    double cycletime=1.0;         /* initial normalized update time */

    double mismatch=0.0;      /* vco mismatch in % of BB */
			      /* mismatches != 0.0 turn off loop integrator */

    double x,y,z;
    double pdel, fdel;

    double window;	      /* triangular windowing fxn for xfer gain */
    double real,imag;	      /* dft scratch variables */

    int random=0;	      /* flag to plot input jitter */
    int errflag=0;
    int verbose=0;
    int titleflag=0;
    int plot_error=0;
    int statonly=0;
    int tristate=0;
    int do_windowing=0;

    int interpolate=1;	      /* how many times to subcompute tupdate */
    int latency=0;            /* number of bits phase detect latency  */

    int maxzero,zerorun;
    int maxone,onerun;
    int olddir;

    time(&time_val);	    /* will be used as seed for rand48 */
    
    opterr = 1;
    progname = argv[0];

    while ((c = getopt(argc,argv,"a:b:c:ef:i:j:l:m:n:o:p:q:rstvx:wz")) != EOF)
	switch (c) {
	case 'a':
	    amp = atof(optarg);
	    break;
	case 'c':
	    cycletime= atof(optarg);
	    break;
	case 'b':
	    bangbang = atof(optarg);
	    break;
	case 'e':
	    plot_error++;
	    break;
	case 'f':
	    period = atof(optarg);
	    break;
	case 'i':
	    interpolate = (int) atof(optarg);
	    break;
	case 'j':
	    djitter = atof(optarg);
	    break;
	case 'l':
	    latency = (int) atof(optarg);
	    if (latency < 0 || latency > MAX_LATENCY-1) {
		fprintf(stderr,"ERROR: -l <latency> out of range (0-15)\n");
		exit(1);
	    }
	    break;
	case 'm':
	    mismatch = atof(optarg);
	    if (mismatch < -1.0 || mismatch > 1.0) {
		fprintf(stderr,"ERROR: -m <vco_mismatch> out of range\n");
		exit(1);
	    }
	    break;
	case 'n':
	    npoints = atof(optarg);
	    break;
	case 'o':
	    offset=atof(optarg);
	    break;
	case 'p':
	    psi = atof(optarg);
	    break;
	case 'q':
	    vjitter = atof(optarg);
	    break;
	case 'r':
	    random++;
	    break;
	case 's':
	    statonly++;
	    break;
	case 't':
	    titleflag++;
	    break;
	case 'v':
	    verbose++;
	    break;
	case 'x':
	    tdensity = atof(optarg);
	    if (tdensity < 0.0) {
		fprintf(stderr,"ERROR: -x <tdensity> out of range\n");
		exit(1);
	    }
	    break;
	case 'w':
	    do_windowing++;
	    break;
	case 'z':
	    tristate++;
	    break;
	case '?':
	    errflag++;
	}

    if (errflag) {
	fprintf(stderr, "usage: %s [options]\n",progname);
	fprintf(stderr, "   -a <input phase mod amplitude>\n");
	fprintf(stderr, "   -b <bangbang time>\n");
	fprintf(stderr, "   -c <absolute cycle time> (default is normalized to 1)\n");
	fprintf(stderr, "   -e (plot phase error rather than phase)\n");
	fprintf(stderr, "   -f <input phase mod period> (in update times)\n");
	fprintf(stderr, "   -i <n> (interpolate data points by n)\n");
	fprintf(stderr, "   -j <gaussian input jitter normalized to bb delta T> (rms)\n");
	fprintf(stderr, "   -l <phase detector latency> (bit-times)\n");
	fprintf(stderr, "   -m <vco mismatch normalized to bb>\n");
	fprintf(stderr, "   -n <number of timesteps to simulate>\n");
	fprintf(stderr, "   -o <vco starting phase> (in bb units)\n");
	fprintf(stderr, "   -p <psi> (2*Beta*Tau/tupdate)\n");
	fprintf(stderr, "   -q <sigma per cycle normalized to bb delta T> (vco phase noise)\n");
	fprintf(stderr, "   -r (plot data phase also)\n");
	fprintf(stderr, "   -s (statistics only)\n");
	fprintf(stderr, "   -t (output tplot title)\n");
	fprintf(stderr, "   -v (verbose flag)\n");
	fprintf(stderr, "   -w (do windowed dft freq & gain calc.)\n");
	fprintf(stderr, "   -x <transition probability>\n");
	fprintf(stderr, "       (if xarg < 1 interpreted as probability)\n");
	fprintf(stderr, "       (if xarg > 1 interpreted as run length)\n");
	fprintf(stderr, "   -z (tristate during non-transition times)\n");
	exit(2);
    }

    if (titleflag) {
	printf("title options: ");
    } else {
	printf("# ");
    }
    for (i=1; i< argc; i++) {
	printf(" %s", argv[i]);
    }
    printf("\n");

    if (titleflag) {
	printf("title BB loop response\n");
	printf("yscale 1 jitter [normalized to BB]\n");
	printf("xscale 1 timestep [normalized to bit time] #\n");
    }

    srand48((long) time_val);	    /* reseed generator */

    /* this next statement sets up steady state phase conditions */
    /* plus any needed delta offset from command line */

    vco_phase = offset+amp-bangbang/2.0;
    dir=0;
    pd_output=0;

    /* initialize pd latency array buffer */
    for (pd_count=0; pd_count<=latency; pd_count++) {
	pd_delay[pd_count]=0;
    }
    pd_count=0;

    maxzero = zerorun=0;
    maxone = onerun=0;
    
    /* for windowing dft analysis */
    real = imag = 0.0;

    jitter=0.0;
    vcojitter=0.0;
    for (cycle=1; cycle<= (int) npoints; cycle++) {

	o_data_phase = data_phase;

	vcojitter += gauss()*vjitter;		/* running integral for VCO phase noise */
	jitter = gauss()*djitter + vcojitter;	/* add in gaussian data phase noise */

	data_phase = jitter +
		     amp*cos(2.0*M_PI*((double) cycle)/period);

	if (!statonly) {
	    for (x=0.0/interpolate; x<1.0; x+=1.0/interpolate) {
		pdel=(dir+mismatch)*bangbang*x;

		if (mismatch == 0.0) {
		    fdel=(((dir+mismatch)*x*x+2.0*loop_filter*x)*bangbang/psi);
		} else {
		    fdel=0.0;
		}

		if (plot_error) { /* plot phase error */
		    printf("%g %g\n", (((double) cycle)-1.0+x)*cycletime, 
		o_data_phase+(data_phase-o_data_phase)*x-(vco_phase+pdel+fdel)); 
			/* data_phase-(vco_phase+pdel+fdel)); */
		} else {       /* plot actual phase trajectory */
		    printf("%g %g\n", (((double) cycle)-1.0+x)*cycletime, 
			vco_phase+pdel+fdel);
		}
	    }
	}

	vco_phase += (dir+mismatch)*bangbang;	 
	if (mismatch == 0.0) { /* turn off integrator when mismatched */
	    vco_phase += 2.0*loop_filter*bangbang/psi;	
	    vco_phase += (dir+mismatch)*bangbang/psi;    
	}

	fflush(stdout);

	olddir=dir;

	if (tdensity >= 1.0) {		    /* run length spec */
	    if (cycle%((int) tdensity) == 0) {
		transition = 1;
	    } else {
		transition = 0;
	    }
	} else if (tdensity < 1.0) {	   /* probabilistic defininition */
	    if (drand48() <= tdensity) {      /* got a transition? */
		transition = 1;
	    } else {
		transition = 0;
	    }
	}

	if (transition) {		/* got a transition? */
	    if (vco_phase >= data_phase) {
		pd_output = -1;
	    } else {
		pd_output = 1;
	    }
	} else {				/* nope */
	    if (tristate) {
		pd_output = 0;			/* tristate */
	    } else {
		; 				/* leave direction alone */
	    }
	}

	pd_delay[pd_count] = pd_output;	
	   /* fprintf(stderr,"1: %d, %d\t", pd_count, pd_output); */
	/* complicated incrementing because (pd_count%0) dumps core */
	pd_count++;
	if (pd_count >latency) {
	    pd_count=0;
	}
	    /* fprintf(stderr,"2: %d, %d\n", pd_count, pd_delay[pd_count]); */
	dir=pd_delay[pd_count];

	loop_filter += (double) dir;

	if (olddir==dir) {
	    if (dir == 1) {
                onerun++;
	    } else {
		zerorun++;
	    }
	} else {
	    if (zerorun > maxzero) {
		maxzero=zerorun;
	    }
	    if (onerun > maxone) {
		maxone=onerun;
	    }
	    zerorun=0;
	    onerun=0;
	}

	rms_sum_vco += vco_phase*vco_phase;
	rms_sum_data += data_phase*data_phase;

	if (do_windowing) {
	    if (cycle > ((int) npoints)/2 ) {
	    window = 1.0-fabs(((double)cycle-3.0*npoints/4.0)/(npoints/4.0));
		real+=cos(2.0*M_PI*((double) cycle)/period)*window*vco_phase;
		imag+=sin(2.0*M_PI*((double) cycle)/period)*window*vco_phase;
	    }
	}

    }

    if (random && !statonly ) {
	srand48((long) time_val);   /* seed generator */

	printf("pen\n");
	jitter=0.0;
	vcojitter=0.0;
	for (cycle=1; cycle<= (int) npoints; cycle++) {

	    vcojitter += gauss()*vjitter;		/* running integral for VCO phase noise */
	    jitter = gauss()*djitter + vcojitter;	/* add in gaussian data phase noise */

	    data_phase = jitter + 
		amp*cos(2.0*M_PI*((double) cycle)/period);
	    printf("%g %g\n", ((double) cycle)*cycletime, data_phase);

	    if (tdensity < 1.0) {	 
	    	drand48();	/* simply called to keep in sync */
	    }
	}

	/* windowing function */
	for (cycle=1; cycle<= (int) npoints; cycle++) {
	    if (cycle > ((int) npoints)/2 ) {
		window = (((double)cycle-3.0*npoints/4.0)/(npoints/4.0));
		/* printf("%d %g\n", cycle, 1.0-fabs(window)); */
	    }
	}
    }

    if (do_windowing) {
	
	/* triangular window normalization */
	x = npoints/8.0;

    printf("# DFT: %g\t%g\t%g\t%g\n",
	    period,
	    amp,
	    sqrt(real*real+imag*imag)/x,
	    360.0*atan2(imag,real))/(2.0*M_PI);
    }

    printf("# STAT: %g\t%g\t%d\t%d\t%g\t%g\n", 
    	djitter, 
	psi, 
	maxone, maxzero,
	sqrt(rms_sum_data/(npoints-1)),
    	sqrt(rms_sum_vco/(npoints-1)) );
    fflush(stdout);

    exit(0);
}

double gauss()
{
    static int iset=0;
    static float gset;
    double fac,r,v1,v2;
    double drand48();

    if(iset == 0) {	/* we don't have an extra deviate handy, so */
	do {
	    v1=2.0*drand48()-1.0;   /* pick two uniform numbers in the */
	    v2=2.0*drand48()-1.0;   /* square from -1 to +1 in each axis */
	    r=v1*v1+v2*v2;	    /* see if they are in unit circle */
	} while (r >= 1.0);	/* and if not, try again... */
	fac=sqrt(-2.0*log(r)/r);

	/* now make the box-muller transformation to get two normal */
	/* deviates. return one and save the other for next time */

	gset=v1*fac;
	iset=1;			/* set flag */
	return v2*fac;
    } else {		/* we have extra deviate handy */
	iset=0;		/* so unset flag */
	return gset;	/* and return it */
    }
}

int prbs23() /* 2^23-1 sequence generator */
{

    static int bitval;
    static int reg[PLEN]    ={ 1,0,0,0,0, 0,0,0,0,0,
			       0,0,0,0,0, 0,0,0,0,0, 0,0,0 }; 
    static int loop,temp;

    bitval = reg[0];
    if (reg[PLEN-1] == reg[TAP-1])  /* do xor */
        temp = 0;
    else
        temp = 1;
    for(loop = 0; loop <= PLEN-1; loop++)
        reg[PLEN-loop-1]=reg[PLEN-loop-2];
    reg[0] = temp;

    return(bitval);
}


/* cyclical pattern generator */

int patgen(index, str) 
int *index	/* index into argument string */;
char *str;	/* pattern string 
		** 0=low, (any other char)=high,
		** [try to use '1' for high... 
		** when EOS is hit, the output cycles over
		** from the beginning.
		*/
{
    *index = (++*index%(strlen(str)));

    if(str[*index] == '0') {
	return(0);
    } else {
	return(1);
    }
}