262 lines
8.6 KiB
C
262 lines
8.6 KiB
C
/* 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");
|
|
}
|