- initial import
git-svn-id: http://moon:8086/svn/mips@1 a8ebac50-d88d-4704-bea3-6648445a41b3
This commit is contained in:
+261
@@ -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");
|
||||
}
|
||||
Reference in New Issue
Block a user