/* 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 #include #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; i0.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; kmax) /* 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"); }