Initial version

Committed on the Free edition of March Hare Software CVSNT Server.
Upgrade to CVS Suite for more features and support:
http://march-hare.com/cvsnt/


git-svn-id: http://moon:8086/svn/vhdl/trunk@309 cc03376c-175c-47c8-b038-4cd826a8556b
This commit is contained in:
2009-02-01 19:37:54 +00:00
parent 3e6b68932f
commit f2ea795156
2 changed files with 669 additions and 0 deletions
+261
View File
@@ -0,0 +1,261 @@
/* dttl.c -- simulation of Cassini/Huygens symbol synchronizer loop
*
* Author: Lorenzo Simone's (matlab script)
* Modifications by: Jon Hamkins (conversion to C, support routines,
* memory requirement reductions)
* Last Revised: Wed Jan 31 10:47:52 PST 2001
*/
#define M_PI 3.1415926535897932384626433832795
#include <math.h>
#include <stdio.h>
#include "random.h"
/* global variables */
int delay,
L; /* samples per symbol */
double Pt, /* transition density */
sigma; /* noise standard deviation */
long seed=-123;
/* This function simulates a large 1-dimensional array without requiring
* storage for the whole array. This is possible because the array is
* accessed in roughly increasing indices. Thus, we only need to
* store a few symbols worth of samples in the array, and we can reuse
* the array indices as time goes on. */
int drint(double v)
{
return (int)(v+0.5);
}
double
rec(int idx)
{
static int firsttime = 1,
max, /* maximum index containing data */
mod, /* length of ring buffer */
b; /* value of last transmitted bit */
static double r[2000]; /* need several symbols worth of samples */
int i,j,k;
/* First time, store samples for random delay and first few symbols */
if (firsttime==1) {
firsttime=0;
for (i=0; i<delay; i++) r[i] = 0; /* random delay */
b = (ran1(&seed)>0.5) ? 1 : -1; /* first random symbol */
for (j=0; j<3; j++) { /* 3 binary symbols with transition Pt */
b = (ran1(&seed)>Pt) ? b : -b; /* next symbol value */
for (k=0; k<L; k++) r[i++] = b + sigma*gaussian(&seed);
}
max = i-1;
mod = i;
}
if (idx>max) /* check if we need to create another sample */
{
b = (ran1(&seed)>Pt) ? b : -b; /* next symbol value */
/* generate samples for one symbol and store in ring buffer */
for (i=1; i<=L; i++) r[(int)(fmod(max+i,mod))] = b+sigma*gaussian(&seed);
max += L; /* maximum index for which samples exist */
}
/* return appropriate sample from ring buffer */
return(r[(int)(fmod(idx,mod))]);
}
main()
{
double
detout, /* detector output */
I, old_I, /* in-phase accumulations */
Q, old_Q, /* quadrature accumulations */
T, /* transition detection */
theta, /* baseband NCO output (radians) */
lambda, /* normalized timing error (symbols) */
A, alfa, Ampl, Bl, Df,
DR, Es_No, Es_No_lin, Fc1, Fc2, Fvco, Fs, Kd, Kv,
lambda_ss_sim = 0., lambda_ss_th, No, P, P_No,
sigma_lambda_sim = 0., sigma_lambda_th, SNR_sim, SNR_th, tau, Tc1, Tc2,
time, Ts, x1, Y, bit_I, bit_Q, bit;
int i, j, k, N, offset, k1, k2, EOB;
FILE *p1, *p2, *p3, *p4, *p5;
/* Data settings */
Fc1 = 48000; /* sample rate at DTTL input (F1/2 = 8.215MHz/2) */
Tc1 = 1./Fc1;
DR = 4800; /* data rate (symb/sec) */
Fc2 = DR; /* sampling rate pre-detection */
Tc2 = 1./Fc2;
L = drint(Tc2/Tc1); /* samples per symbol in the arm filters */
offset = drint(L/2.);
Es_No = 600.; /* dB */
P_No = Es_No+10.*log10(DR); /* data power-over-noise spectral density */
/* ratio (dBHz) */
N = 2500; /* number of symbols to simulate */
Pt=.5; /* Transition density */
Bl=0.01*DR; /* loop bandwidth (Hz) */
Df=0.001*DR; /* frequency offset */
Fvco=DR+Df;
tau=ran1(&seed)-.5; /* set symbol timing offset, -.5 to .5 */
delay=drint(2.*L+tau*L+1); /* channel delay */
P=1.; /* data power */
Ampl=sqrt(P);
No=10.*log10(P*Tc2)-Es_No; /* noise spectral density (dBW/Hz) */
sigma=sqrt(pow(10.,0.1*No)/(2.*Tc1));
printf("sigma=%f\n",sigma);
/* DPLL settings */
Fs=DR; /* loop sampling frequency (Hz) */
Ts=1./Fs; /* loop sampling time (s) */
/* Loop gain */
Kd=2.*Ampl*Pt; /* phase detector gain (quants/rad) */
/* @ high Es/No */
Kv=1.; /* VCO gain (rad/quants) */
/* Analog loop filter components */
alfa=4.*Bl/(Kd*Kv);
/* Digital loop filter components */
A=alfa*Ts;
/* Initial condition */
theta = 0;
old_I=1;
old_Q=1;
x1=0;
Y=0;
/* Plot data files */
// p1 = fopen("plot1.dat","w"); /* theta */
// p2 = fopen("plot2.dat","w"); /* lambda */
// p3 = fopen("plot3.dat","w"); /* detout */
// p4 = fopen("plot4.dat","w"); /* I-Data */
// p5 = fopen("plot5.dat","w"); /* Q-Data */
k1 = drint(0.5*N);
k2 = N;
printf("Data rate: %.0f symbols/sec\n",DR);
printf("Sample rate at DTTL input: %.0f samples/sec\n",Fc1);
printf("Samples per symbol: %d\n",L);
printf("Random delay set to: %f symbols (%d samples)\n",tau,delay);
printf("Transition density used in random bit generation: %f\n",Pt);
printf("Es/No: %f dB\n",Es_No);
printf("Loop bandwidth: %f Hz\n",Bl);
printf("Doppler: %f Hz\n",Df);
printf("Simulation length: %d symbols (%d samples)\n\n",N,N*L+delay);
/* Simulation Loop */
for (i=2; i<=N; i++)
{
/* loop over symbols */
/* find start of in-phase integration */
for (j=1; j < L; j++)
{
if (fmod(theta + 2*M_PI * Fvco*(j+1+(i-1)*(Fc1/DR))/Fc1, 2*M_PI)
< fmod(theta + 2*M_PI * Fvco*(j +(i-1)*(Fc1/DR))/Fc1, 2*M_PI))
break;
/* one could replace loop above with a direct calculation such as: */
/* j=ceil((1.-modulo(theta/(2*M_PI),1.))/(Fvco*Tc1))-1.; */
/* j+1 is the first index in the block */
EOB=j+(i-1)*L; /* End of Block */
}
/* In-phase integrator */
bit_I = 0.;
printf("Symbol # %d\n", i);
printf("I[%d:%d] :\n", EOB-L+1, EOB);
for (j=EOB-L+1; j<=EOB; j++)
{
bit = rec(j);
if (bit > 0)
printf("-");
else
printf("_");
bit_I += bit;
}
printf("\n");
I = bit_I;
/* Mid-Phase integrator */
bit_Q = 0.;
printf("Q[%d:%d] :\n", EOB-L+offset+1, EOB+offset);
for (j=EOB-L+offset+1; j<=EOB+offset; j++)
{
bit = rec(j);
if (bit > 0)
printf("-");
else
printf("_");
bit_Q += bit;
}
printf("\n\n");
Q = bit_Q;
Q *= 2.*M_PI/(Tc2/Tc1);
/* Transition Detector */
I = (I > 0) ? 1. : -1.; /* hard limiter */
T = 0.5*(I - old_I);
/* Detector Output */
detout = fmod(old_Q*T,2*M_PI);
/* Save accumulations for next time through loop */
old_I = I;
old_Q = Q;
if (fabs(detout)>M_PI)
detout = -detout+2*M_PI*((detout>0) ? 1. : -1.);
/* NCO */
Y += detout*A; /* Loop filter */
theta=Kv*Y; /* phase */
// lambda = theta/(2.*M_PI)+(i*Tc2*Df+tau); /* normalized error */
lambda = theta/(2.*M_PI)+(i*(Fvco/DR -1) + tau); /* normalized error */
/* error mean and variance calculation */
lambda_ss_sim += lambda; /* partial sum */
sigma_lambda_sim += lambda*lambda; /* partial sum of squares */
/* store phase and normalized error to data files */
// fprintf(p1,"%f %f\n",(double)i,theta);
// fprintf(p2,"%f %f\n",(double)i,lambda);
// fprintf(p3,"%f %f\n",(double)i,detout);
// fprintf(p4,"%f %f\n",(double)i,I);
// fprintf(p5,"%f %f\n",(double)i,Q);
}
// fclose(p1);
// fclose(p2);
// fclose(p3);
// fclose(p4);
// fclose(p5);
/* Statistics Analysis */
/* ------------------- */
/* complete the mean and variance calculation */
lambda_ss_sim /= N;
sigma_lambda_sim = sigma_lambda_sim/N - lambda_ss_sim*lambda_ss_sim;
Es_No_lin = pow(10.,0.1*Es_No);
// sigma_lambda_th = Bl*Tc2/(4.*Pt*Es_No_lin*pow(erf(sqrt(Es_No_lin)),2.));
sigma_lambda_th = Bl*Tc2/(4.*Pt*Es_No_lin*Es_No_lin);
SNR_sim = -10*log10(sigma_lambda_sim);
SNR_th = -10*log10(sigma_lambda_th);
/* Steady-State error */
lambda_ss_th=Df/(4.*Bl);
printf("Simulation: var: %f rad^2 loop SNR: %f dB ss error: %f%%\n",
sigma_lambda_sim, SNR_sim, 100.*lambda_ss_sim);
printf("Theory: var: %f rad^2 loop SNR: %f dB ss error: %f%%\n",
sigma_lambda_th, SNR_th, 100.*lambda_ss_th);
// system("gnuplot < plot1.gp");
// system("gnuplot < plot2.gp");
printf("Done.\n");
}
@@ -0,0 +1,408 @@
/* C version of the systems programming language benchmark
** Author: M. J. Jordan Cambridge Computer Laboratory.
**
** Modified by: M. Richards, Nov 1996
** to be ANSI C and runnable on 64 bit machines + other minor changes
** Modified by: M. Richards, 20 Oct 1998
** made minor corrections to improve ANSI compliance (suggested
** by David Levine)
**
** Compile with, say
**
** gcc -o bench bench.c
**
** or
**
** gcc -o bench100 -Dbench100 bench.c (for a version that obeys
** the main loop 100x more often)
*/
#include <stdio.h>
#include <stdlib.h>
#ifdef SMALL_PROBLEM_SIZE
#define Count 1000*1000
#define Qpktcountval 2326389
#define Holdcountval 930555
#else
#define Count 10000*1000
#define Qpktcountval 23263894
#define Holdcountval 9305557
#endif
#define TRUE 1
#define FALSE 0
#define MAXINT 32767
#define BUFSIZE 3
#define I_IDLE 1
#define I_WORK 2
#define I_HANDLERA 3
#define I_HANDLERB 4
#define I_DEVA 5
#define I_DEVB 6
#define PKTBIT 1
#define WAITBIT 2
#define HOLDBIT 4
#define NOTPKTBIT !1
#define NOTWAITBIT !2
#define NOTHOLDBIT 0XFFFB
#define S_RUN 0
#define S_RUNPKT 1
#define S_WAIT 2
#define S_WAITPKT 3
#define S_HOLD 4
#define S_HOLDPKT 5
#define S_HOLDWAIT 6
#define S_HOLDWAITPKT 7
#define K_DEV 1000
#define K_WORK 1001
struct packet
{
struct packet *p_link;
int p_id;
int p_kind;
int p_a1;
char p_a2[4];
};
struct task
{
struct task *t_link;
int t_id;
int t_pri;
struct packet *t_wkq;
int t_state;
struct task *(*t_fn)(struct packet *);
long t_v1;
long t_v2;
};
char alphabet[28] = "0ABCDEFGHIJKLMNOPQRSTUVWXYZ";
struct task *tasktab[11] = {0,0,0,0,0,0,0,0,0,0,0};
struct task *tasklist = 0;
struct task *tcb;
long taskid;
long v1;
long v2;
int qpktcount = 0;
int holdcount = 0;
int tracing = 1;
int layout = 0;
void append(struct packet *pkt, struct packet *ptr);
void createtask(int id,
int pri,
struct packet *wkq,
int state,
struct task *(*fn)(struct packet *),
long v1,
long v2)
{
struct task *t = (struct task *)malloc(sizeof(struct task));
tasktab[id] = t;
t->t_link = tasklist;
t->t_id = id;
t->t_pri = pri;
t->t_wkq = wkq;
t->t_state = state;
t->t_fn = fn;
t->t_v1 = v1;
t->t_v2 = v2;
tasklist = t;
}
struct packet *pkt(struct packet *link, int id, int kind)
{
int i;
struct packet *p = (struct packet *)malloc(sizeof(struct packet));
for (i=0; i<=BUFSIZE; i++)
p->p_a2[i] = 0;
p->p_link = link;
p->p_id = id;
p->p_kind = kind;
p->p_a1 = 0;
return (p);
}
void trace(char a)
{
if ( --layout <= 0 )
{
printf("\n");
layout = 50;
}
printf("%c", a);
}
void schedule()
{
while ( tcb != 0 )
{
struct packet *pkt;
struct task *newtcb;
pkt=0;
switch ( tcb->t_state )
{
case S_WAITPKT:
pkt = tcb->t_wkq;
tcb->t_wkq = pkt->p_link;
tcb->t_state = tcb->t_wkq == 0 ? S_RUN : S_RUNPKT;
case S_RUN:
case S_RUNPKT:
taskid = tcb->t_id;
v1 = tcb->t_v1;
v2 = tcb->t_v2;
if (tracing==TRUE) trace(taskid+'0');
newtcb = (*(tcb->t_fn))(pkt);
tcb->t_v1 = v1;
tcb->t_v2 = v2;
tcb = newtcb;
break;
case S_WAIT:
case S_HOLD:
case S_HOLDPKT:
case S_HOLDWAIT:
case S_HOLDWAITPKT:
tcb = tcb->t_link;
break;
default:
return;
}
}
}
struct task *Wait(void)
{
tcb->t_state |= WAITBIT;
return (tcb);
}
struct task *holdself(void)
{
++holdcount;
tcb->t_state |= HOLDBIT;
return (tcb->t_link) ;
}
struct task *findtcb(int id)
{
struct task *t = 0;
if (1<=id && id<=(long)10)
t = tasktab[id];
if (t==0) printf("\nBad task id %d\n", id);
return(t);
}
struct task *release(int id)
{
struct task *t;
t = findtcb(id);
if ( t==0 ) return (0);
t->t_state &= NOTHOLDBIT;
if ( t->t_pri > tcb->t_pri ) return (t);
return (tcb) ;
}
struct task *qpkt(struct packet *pkt)
{
struct task *t;
t = findtcb(pkt->p_id);
if (t==0) return (t);
qpktcount++;
pkt->p_link = 0;
pkt->p_id = taskid;
if (t->t_wkq==0)
{
t->t_wkq = pkt;
t->t_state |= PKTBIT;
if (t->t_pri > tcb->t_pri) return (t);
}
else
{
append(pkt, (struct packet *)&(t->t_wkq));
}
return (tcb);
}
struct task *idlefn(struct packet *pkt)
{
--v2;
if ( v2==0 ) return ( holdself() );
if ( (v1&1) == 0 )
{
v1 = ( v1>>1) & MAXINT;
return ( release(I_DEVA) );
}
else
{
v1 = ( (v1>>1) & MAXINT) ^ 0XD008;
return ( release(I_DEVB) );
}
}
struct task *workfn(struct packet *pkt)
{
if ( pkt==0 ) return ( Wait() );
else
{
int i;
v1 = I_HANDLERA + I_HANDLERB - v1;
pkt->p_id = v1;
pkt->p_a1 = 0;
for (i=0; i<=BUFSIZE; i++)
{
v2++;
if ( v2 > 26 ) v2 = 1;
(pkt->p_a2)[i] = alphabet[v2];
}
return ( qpkt(pkt) );
}
}
struct task *handlerfn(struct packet *pkt)
{
if ( pkt!=0) append(pkt,
(struct packet *)(pkt->p_kind==K_WORK ? &v1 : &v2));
if ( v1!=0 )
{
int count;
struct packet *workpkt = (struct packet *)v1;
count = workpkt->p_a1;
if ( count > BUFSIZE )
{
v1 = (long)(((struct packet *)v1)->p_link);
return ( qpkt(workpkt) );
}
if ( v2!=0 )
{
struct packet *devpkt;
devpkt = (struct packet *)v2;
v2 = (long)(((struct packet *)v2)->p_link);
devpkt->p_a1 = workpkt->p_a2[count];
workpkt->p_a1 = count+1;
return( qpkt(devpkt) );
}
}
return ( Wait() );
}
struct task *devfn(struct packet *pkt)
{
if ( pkt==0 )
{
if ( v1==0 ) return ( Wait() );
pkt = (struct packet *)v1;
v1 = 0;
return ( qpkt(pkt) );
}
else
{
v1 = (long)pkt;
if (tracing==TRUE) trace(pkt->p_a1);
return ( holdself() );
}
}
void append(struct packet *pkt, struct packet *ptr)
{
pkt->p_link = 0;
while ( ptr->p_link ) ptr = ptr->p_link;
ptr->p_link = pkt;
}
int main()
{
struct packet *wkq = 0;
int retval;
printf("Bench mark starting\n");
createtask(I_IDLE, 0, wkq, S_RUN, idlefn, 1, Count);
wkq = pkt(0, 0, K_WORK);
wkq = pkt(wkq, 0, K_WORK);
createtask(I_WORK, 1000, wkq, S_WAITPKT, workfn, I_HANDLERA, 0);
wkq = pkt(0, I_DEVA, K_DEV);
wkq = pkt(wkq, I_DEVA, K_DEV);
wkq = pkt(wkq, I_DEVA, K_DEV);
createtask(I_HANDLERA, 2000, wkq, S_WAITPKT, handlerfn, 0, 0);
wkq = pkt(0, I_DEVB, K_DEV);
wkq = pkt(wkq, I_DEVB, K_DEV);
wkq = pkt(wkq, I_DEVB, K_DEV);
createtask(I_HANDLERB, 3000, wkq, S_WAITPKT, handlerfn, 0, 0);
wkq = 0;
createtask(I_DEVA, 4000, wkq, S_WAIT, devfn, 0, 0);
createtask(I_DEVB, 5000, wkq, S_WAIT, devfn, 0, 0);
tcb = tasklist;
qpktcount = holdcount = 0;
printf("Starting\n");
tracing = FALSE;
layout = 0;
schedule();
printf("finished\n");
printf("qpkt count = %d holdcount = %d\n",
qpktcount, holdcount);
printf("These results are ");
if (qpktcount == Qpktcountval && holdcount == Holdcountval) {
printf("correct");
retval = 0;
} else {
printf("incorrect");
retval = 1;
}
printf("\nend of run\n");
return retval;
}