1 //
2 // bpresync_example.c
3 //
4 // This example demonstrates the binary pre-demodulator synchronizer. A random
5 // binary sequence is generated, modulated with BPSK, and then interpolated.
6 // The resulting sequence is used to generate a bpresync object which in turn
7 // is used to detect a signal in the presence of carrier frequency and timing
8 // offsets and additive white Gauss noise.
9 //
10 
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <getopt.h>
14 #include <math.h>
15 #include <time.h>
16 #include "liquid.h"
17 
18 #define OUTPUT_FILENAME "bpresync_example.m"
19 
20 // print usage/help message
usage()21 void usage()
22 {
23     printf("bpresync_example -- test binary pre-demodulation synchronization\n");
24     printf("options (default values in <>):\n");
25     printf("  h     : print usage/help\n");
26     printf("  k     : samples/symbol, default: 2\n");
27     printf("  m     : filter delay [symbols], default: 5\n");
28     printf("  n     : number of data symbols, default: 64\n");
29     printf("  b     : bandwidth-time product beta in (0,1), default: 0.3\n");
30     printf("  F     : carrier frequency offset, default: 0.02\n");
31     printf("  t     : fractional sample offset dt in [-0.5, 0.5], default: 0\n");
32     printf("  S     : SNR [dB], default: 20\n");
33 }
34 
main(int argc,char * argv[])35 int main(int argc, char*argv[])
36 {
37     srand(time(NULL));
38 
39     // options
40     unsigned int k=2;                   // filter samples/symbol
41     unsigned int m=5;                   // filter delay (symbols)
42     float beta=0.3f;                    // bandwidth-time product
43     float dt = 0.0f;                    // fractional sample timing offset
44     unsigned int num_sync_symbols = 64; // number of synchronization symbols
45     float SNRdB = 20.0f;                // signal-to-noise ratio [dB]
46     float dphi = 0.02f;                 // carrier frequency offset
47     float phi  = 2*M_PI*randf();        // carrier phase offset
48 
49     int dopt;
50     while ((dopt = getopt(argc,argv,"uhk:m:n:b:t:F:t:S:")) != EOF) {
51         switch (dopt) {
52         case 'h': usage();                          return 0;
53         case 'k': k = atoi(optarg);                 break;
54         case 'm': m = atoi(optarg);                 break;
55         case 'n': num_sync_symbols = atoi(optarg);  break;
56         case 'b': beta = atof(optarg);              break;
57         case 'F': dphi = atof(optarg);              break;
58         case 't': dt = atof(optarg);                break;
59         case 'S': SNRdB = atof(optarg);             break;
60         default:
61             exit(1);
62         }
63     }
64 
65     unsigned int i;
66 
67     // validate input
68     if (beta <= 0.0f || beta >= 1.0f) {
69         fprintf(stderr,"error: %s, bandwidth-time product must be in (0,1)\n", argv[0]);
70         exit(1);
71     } else if (dt < -0.5f || dt > 0.5f) {
72         fprintf(stderr,"error: %s, fractional sample offset must be in [-0.5,0.5]\n", argv[0]);
73         exit(1);
74     }
75 
76     // derived values
77     unsigned int num_symbols = num_sync_symbols + 2*m + 10;
78     unsigned int num_samples = k*num_symbols;
79     float nstd = powf(10.0f, -SNRdB/20.0f);
80 
81     // arrays
82     float complex seq[num_sync_symbols];    // synchronization pattern (symbols)
83     float complex s0[k*num_sync_symbols];   // synchronization pattern (samples)
84     float complex x[num_samples];           // transmitted signal
85     float complex y[num_samples];           // received signal
86     float complex rxy[num_samples];         // pre-demod correlation output
87     float dphi_hat[num_samples];            // carrier offset estimate
88 
89     // create transmit/receive interpolator/decimator
90     firinterp_crcf interp = firinterp_crcf_create_prototype(LIQUID_FIRFILT_RRC,k,m,beta,dt);
91 
92     // generate synchronization pattern (BPSK) and interpolate
93     for (i=0; i<num_sync_symbols + 2*m; i++) {
94         float complex sym = 0.0f;
95 
96         if (i < num_sync_symbols) {
97             sym = rand() % 2 ? -1.0f : 1.0f;
98             seq[i] = sym;
99         }
100 
101         if (i < 2*m) firinterp_crcf_execute(interp, sym, s0);
102         else         firinterp_crcf_execute(interp, sym, &s0[k*(i-2*m)]);
103     }
104 
105     // reset interpolator
106     firinterp_crcf_reset(interp);
107 
108     // interpolate input
109     for (i=0; i<num_symbols; i++) {
110         float complex sym = i < num_sync_symbols ? seq[i] : 0.0f;
111 
112         firinterp_crcf_execute(interp, sym, &x[k*i]);
113     }
114 
115     // push through channel
116     for (i=0; i<num_samples; i++)
117         y[i] = x[i]*cexpf(_Complex_I*(dphi*i + phi)) + nstd*(randnf() + _Complex_I*randnf())*M_SQRT1_2;
118 
119     // create cross-correlator
120     bpresync_cccf sync = bpresync_cccf_create(s0, k*num_sync_symbols, 0.05f, 11);
121     bpresync_cccf_print(sync);
122 
123     // push signal through cross-correlator
124     float rxy_max  = 0.0f;  // maximum cross-correlation
125     float dphi_est = 0.0f;  // carrier frequency offset estimate
126     int delay_est  = 0;     // delay estimate
127     for (i=0; i<num_samples; i++) {
128 
129         // correlate
130         bpresync_cccf_push(sync, y[i]);
131         bpresync_cccf_execute(sync, &rxy[i], &dphi_hat[i]);
132 
133         // detect...
134         if (cabsf(rxy[i]) > 0.6f) {
135             printf("****** preamble found, rxy = %12.8f (dphi-hat: %12.8f), i=%3u ******\n",
136                     cabsf(rxy[i]), dphi_hat[i], i);
137         }
138 
139         // retain maximum
140         if (cabsf(rxy[i]) > rxy_max) {
141             rxy_max   = cabsf(rxy[i]);
142             dphi_est  = dphi_hat[i];
143             delay_est = (int)i - (int)2*k*m + 1;
144         }
145     }
146 
147     // destroy objects
148     firinterp_crcf_destroy(interp);
149     bpresync_cccf_destroy(sync);
150 
151     // print results
152     printf("\n");
153     printf("rxy (max) : %12.8f\n", rxy_max);
154     printf("dphi est. : %12.8f ,error=%12.8f\n",      dphi_est, dphi-dphi_est);
155     printf("delay est.: %12d ,error=%3d sample(s)\n", delay_est, k*num_sync_symbols - delay_est);
156     printf("\n");
157 
158     //
159     // export results
160     //
161     FILE * fid = fopen(OUTPUT_FILENAME,"w");
162     fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
163     fprintf(fid,"clear all\n");
164     fprintf(fid,"close all\n");
165     fprintf(fid,"num_samples = %u;\n", num_samples);
166     fprintf(fid,"num_symbols = %u;\n", num_symbols);
167     fprintf(fid,"k           = %u;\n", k);
168 
169     fprintf(fid,"x   = zeros(1,num_samples);\n");
170     fprintf(fid,"y   = zeros(1,num_samples);\n");
171     fprintf(fid,"rxy = zeros(1,num_samples);\n");
172     for (i=0; i<num_samples; i++) {
173         fprintf(fid,"x(%4u)     = %12.8f + j*%12.8f;\n", i+1, crealf(x[i]),   cimagf(x[i]));
174         fprintf(fid,"y(%4u)     = %12.8f + j*%12.8f;\n", i+1, crealf(y[i]),   cimagf(y[i]));
175         fprintf(fid,"rxy(%4u)   = %12.8f + j*%12.8f;\n", i+1, crealf(rxy[i]), cimagf(rxy[i]));
176     }
177 
178     fprintf(fid,"t=[0:(num_samples-1)]/k;\n");
179     fprintf(fid,"figure;\n");
180     fprintf(fid,"subplot(2,1,1);\n");
181     fprintf(fid,"  plot(t,real(y), t,imag(y));\n");
182     fprintf(fid,"  axis([0 num_symbols -2 2]);\n");
183     fprintf(fid,"  grid on;\n");
184     fprintf(fid,"  xlabel('time');\n");
185     fprintf(fid,"  ylabel('received signal');\n");
186     fprintf(fid,"subplot(2,1,2);\n");
187     fprintf(fid,"  plot(t,abs(rxy));\n");
188     fprintf(fid,"  axis([0 num_symbols 0 1.5]);\n");
189     fprintf(fid,"  xlabel('time');\n");
190     fprintf(fid,"  ylabel('correlator output');\n");
191     fprintf(fid,"  grid on;\n");
192 
193     fclose(fid);
194     printf("results written to '%s'\n", OUTPUT_FILENAME);
195 
196     return 0;
197 }
198