1 //
2 // predemod_sync_test.c
3 //
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <getopt.h>
8 #include <math.h>
9 #include <assert.h>
10 #include <time.h>
11 #include "liquid.h"
12
13 #define OUTPUT_FILENAME "predemod_sync_test.m"
14
15 // print usage/help message
usage()16 void usage()
17 {
18 printf("predemod_sync_test -- test pre-demodulation synchronization\n");
19 printf("options:\n");
20 printf(" h : print usage/help\n");
21 printf(" k : samples/symbol, default: 2\n");
22 printf(" m : filter delay [symbols], default: 4\n");
23 printf(" n : number of data symbols, default: 64\n");
24 printf(" b : bandwidth-time product, (0,1), default: 0.3\n");
25 printf(" t : fractional sample offset, (-0.5,0.5), default: 0\n");
26 printf(" F : frequency offset, default: 0\n");
27 printf(" P : phase offset, default: 0\n");
28 printf(" s : SNR [dB], default: 30\n");
29 }
30
main(int argc,char * argv[])31 int main(int argc, char*argv[])
32 {
33 srand(time(NULL));
34 // options
35 unsigned int k=2; // filter samples/symbol
36 unsigned int m=4; // filter delay (symbols)
37 float beta=0.3f; // bandwidth-time product
38 float dt = 0.0f; // fractional sample timing offset
39 unsigned int num_sync_symbols = 64; // number of data symbols
40 float SNRdB = 30.0f; // signal-to-noise ratio [dB]
41 float dphi = 0.0f; // carrier frequency offset
42 float phi = 0.0f; // carrier phase offset
43
44 unsigned int num_delay_symbols = 12;
45 unsigned int num_dphi_hat = 21; // number of frequency offset estimates
46 float dphi_hat_step = 0.01f; // frequency offset step size
47
48 int dopt;
49 while ((dopt = getopt(argc,argv,"uhk:m:n:b:t:F:P:s:")) != EOF) {
50 switch (dopt) {
51 case 'h': usage(); return 0;
52 case 'k': k = atoi(optarg); break;
53 case 'm': m = atoi(optarg); break;
54 case 'n': num_sync_symbols = atoi(optarg); break;
55 case 'b': beta = atof(optarg); break;
56 case 't': dt = atof(optarg); break;
57 case 'F': dphi = atof(optarg); break;
58 case 'P': phi = 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.5 || dt > 0.5) {
72 fprintf(stderr,"error: %s, fractional sample offset must be in (0,1)\n", argv[0]);
73 exit(1);
74 }
75
76 // derived values
77 unsigned int num_symbols = num_delay_symbols + num_sync_symbols + 2*m;
78 unsigned int num_samples = k*num_symbols;
79 unsigned int num_sync_samples = k*num_sync_symbols;
80 float nstd = powf(10.0f, -SNRdB/20.0f);
81
82 // arrays
83 float complex seq[num_sync_symbols]; // data sequence (symbols)
84 float complex s0[num_sync_samples]; // data sequence (interpolated samples)
85 float complex x[num_samples]; // transmitted signal
86 float complex y[num_samples]; // received signal
87 float rxy[num_dphi_hat][num_samples]; // pre-demod output matrix
88
89 // generate sequence
90 for (i=0; i<num_sync_symbols; i++) {
91 float sym_i = rand() % 2 ? M_SQRT1_2 : -M_SQRT1_2;
92 float sym_q = rand() % 2 ? M_SQRT1_2 : -M_SQRT1_2;
93 seq[i] = sym_i + _Complex_I*sym_q;
94 }
95
96 // create interpolated sequence, compensating for filter delay
97 firinterp_crcf interp_seq = firinterp_crcf_create_prototype(LIQUID_FIRFILT_RRC,k,m,beta,0.0f);
98 for (i=0; i<num_sync_symbols+m; i++) {
99 if (i < m) firinterp_crcf_execute(interp_seq, seq[i], &s0[0]);
100 else if (i < num_sync_symbols) firinterp_crcf_execute(interp_seq, seq[i], &s0[k*(i-m)]);
101 else firinterp_crcf_execute(interp_seq, 0, &s0[k*(i-m)]);
102 }
103 firinterp_crcf_destroy(interp_seq);
104
105 // compute g = E{ |s0|^2 }
106 float g = 0.0f;
107 for (i=0; i<num_sync_samples; i++)
108 g += crealf( s0[i]*conjf(s0[i]) );
109
110 // create transmit interpolator and generate sequence
111 firinterp_crcf interp_tx = firinterp_crcf_create_prototype(LIQUID_FIRFILT_RRC,k,m,beta,dt);
112 unsigned int n=0;
113 for (i=0; i<num_delay_symbols; i++) {
114 firinterp_crcf_execute(interp_tx, 0, &x[k*n]);
115 n++;
116 }
117 for (i=0; i<num_sync_symbols; i++) {
118 firinterp_crcf_execute(interp_tx, seq[i], &x[k*n]);
119 n++;
120 }
121 for (i=0; i<2*m; i++) {
122 firinterp_crcf_execute(interp_tx, 0, &x[k*n]);
123 n++;
124 }
125 assert(n==num_symbols);
126 firinterp_crcf_destroy(interp_tx);
127
128 // add channel impairments
129 for (i=0; i<num_samples; i++) {
130 y[i] = x[i]*cexp(_Complex_I*(dphi*i + phi)) + nstd*( randnf() + _Complex_I*randnf() );
131 }
132
133 float complex z; // filter output sample
134 for (n=0; n<num_dphi_hat; n++) {
135 float dphi_hat = ((float)n - 0.5*(float)(num_dphi_hat-1)) * dphi_hat_step;
136 printf(" dphi_hat : %12.8f\n", dphi_hat);
137
138 // create flipped, conjugated coefficients
139 float complex s1[num_sync_samples];
140 for (i=0; i<num_sync_samples; i++)
141 s1[i] = conjf( s0[num_sync_samples-i-1]*cexpf(_Complex_I*(dphi_hat*i)) );
142
143 // create matched filter and detect signal
144 firfilt_cccf fsync = firfilt_cccf_create(s1, num_sync_samples);
145 for (i=0; i<num_samples; i++) {
146 firfilt_cccf_push(fsync, y[i]);
147 firfilt_cccf_execute(fsync, &z);
148
149 rxy[n][i] = cabsf(z) / g;
150 }
151 // destroy filter
152 firfilt_cccf_destroy(fsync);
153 }
154
155 // print results
156 //printf("rxy (max) : %12.8f\n", rxy_max);
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,"k = %u;\n", k);
166 fprintf(fid,"m = %u;\n", m);
167 fprintf(fid,"beta = %f;\n", beta);
168 fprintf(fid,"num_sync_symbols = %u;\n", num_sync_symbols);
169 fprintf(fid,"num_sync_samples = k*num_sync_symbols;\n");
170 fprintf(fid,"num_symbols = %u;\n", num_symbols);
171 fprintf(fid,"num_samples = %u;\n", num_samples);
172 fprintf(fid,"num_dphi_hat = %u;\n", num_dphi_hat);
173 fprintf(fid,"dphi_hat_step = %f;\n", dphi_hat_step);
174
175 // save sequence symbols
176 fprintf(fid,"seq = zeros(1,num_sync_symbols);\n");
177 for (i=0; i<num_sync_symbols; i++)
178 fprintf(fid,"seq(%4u) = %12.8f + j*%12.8f;\n", i+1, crealf(seq[i]), cimagf(seq[i]));
179
180 // save interpolated sequence
181 fprintf(fid,"s = zeros(1,num_sync_samples);\n");
182 for (i=0; i<num_sync_samples; i++)
183 fprintf(fid,"s(%4u) = %12.8f + j*%12.8f;\n", i+1, crealf(s0[i]), cimagf(s0[i]));
184
185 fprintf(fid,"x = zeros(1,num_samples);\n");
186 fprintf(fid,"y = zeros(1,num_samples);\n");
187 for (i=0; i<num_samples; i++) {
188 fprintf(fid,"x(%6u) = %12.8f + j*%12.8f;\n", i+1, crealf(x[i]), cimagf(x[i]));
189 fprintf(fid,"y(%6u) = %12.8f + j*%12.8f;\n", i+1, crealf(y[i]), cimagf(y[i]));
190 }
191
192 // save cross-correlation output
193 fprintf(fid,"rxy = zeros(num_dphi_hat,num_samples);\n");
194 for (n=0; n<num_dphi_hat; n++) {
195 for (i=0; i<num_samples; i++) {
196 fprintf(fid,"rxy(%6u,%6u) = %12.8f;\n", n+1, i+1, rxy[n][i]);
197 }
198 }
199 fprintf(fid,"t=[0:(num_samples-1)]/k;\n");
200 fprintf(fid,"figure;\n");
201 fprintf(fid,"plot(1:length(s),real(s), 1:length(s),imag(s));\n");
202
203 fprintf(fid,"dphi_hat = ( [0:(num_dphi_hat-1)] - (num_dphi_hat-1)/2 ) * dphi_hat_step;\n");
204 fprintf(fid,"mesh(dphi_hat, t, rxy');\n");
205
206 #if 0
207 fprintf(fid,"z = abs( z );\n");
208 fprintf(fid,"[zmax i] = max(z);\n");
209 fprintf(fid,"plot(1:length(z),z,'-x');\n");
210 fprintf(fid,"axis([(i-8*k) (i+8*k) 0 zmax*1.2]);\n");
211 fprintf(fid,"grid on\n");
212 #endif
213
214 fclose(fid);
215 printf("results written to '%s'\n", OUTPUT_FILENAME);
216
217 return 0;
218 }
219