1 /*---------------------------------------------------------------------------*\
2
3 FILE........: ofdm.c
4 AUTHORS.....: David Rowe & Steve Sampson
5 DATE CREATED: June 2017
6
7 A Library of functions that implement a PSK OFDM modem, C port of
8 the Octave functions in ofdm_lib.m
9
10 \*---------------------------------------------------------------------------*/
11
12 /*
13 Copyright (C) 2017-2020 David Rowe
14
15 All rights reserved.
16
17 This program is free software; you can redistribute it and/or modify
18 it under the terms of the GNU Lesser General Public License version 2.1, as
19 published by the Free Software Foundation. This program is
20 distributed in the hope that it will be useful, but WITHOUT ANY
21 WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
23 License for more details.
24
25 You should have received a copy of the GNU Lesser General Public License
26 along with this program; if not, see <http://www.gnu.org/licenses/>.
27 */
28
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <stdbool.h>
32 #include <stdint.h>
33 #include <string.h>
34 #include <math.h>
35 #include <assert.h>
36 #include <complex.h>
37
38 #include "comp.h"
39 #include "ofdm_internal.h"
40 #include "codec2_ofdm.h"
41 #include "filter.h"
42 #include "wval.h"
43 #include "debug_alloc.h"
44 #include "machdep.h"
45
46 #ifdef __EMBEDDED__
47 #include "arm_math.h"
48 #endif /* __EMBEDDED__ */
49
50 /* Static Prototypes */
51
52 static float cnormf(complex float);
53 static void allocate_tx_bpf(struct OFDM *);
54 static void deallocate_tx_bpf(struct OFDM *);
55 static void dft(struct OFDM *, complex float *, complex float *);
56 static void idft(struct OFDM *, complex float *, complex float *);
57 static complex float vector_sum(complex float *, int);
58 static int est_timing(struct OFDM *, complex float *, int, int, float *, int *, int);
59 static float est_freq_offset_pilot_corr(struct OFDM *, complex float *, int, int);
60 static int ofdm_sync_search_core(struct OFDM *);
61 static void ofdm_demod_core(struct OFDM *, int *);
62
63 /* Defines */
64
65 #define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
66 #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
67
68 /*
69 * QPSK Quadrant bit-pair values - Gray Coded
70 */
71 static const complex float qpsk[] = {
72 1.0f + 0.0f * I,
73 0.0f + 1.0f * I,
74 0.0f - 1.0f * I,
75 -1.0f + 0.0f * I
76 };
77
78 static const complex float qam16[] = {
79 1.0f + 1.0f * I,
80 1.0f + 3.0f * I,
81 3.0f + 1.0f * I,
82 3.0f + 3.0f * I,
83 1.0f - 1.0f * I,
84 1.0f - 3.0f * I,
85 3.0f - 1.0f * I,
86 3.0f - 3.0f * I,
87 -1.0f + 1.0f * I,
88 -1.0f + 3.0f * I,
89 -3.0f + 1.0f * I,
90 -3.0f + 3.0f * I,
91 -1.0f - 1.0f * I,
92 -1.0f - 3.0f * I,
93 -3.0f - 1.0f * I,
94 -3.0f - 3.0f * I
95 };
96
97 /*
98 * These pilots are compatible with Octave version
99 */
100 static const int8_t pilotvalues[] = {
101 -1,-1, 1, 1,-1,-1,-1, 1,
102 -1, 1,-1, 1, 1, 1, 1, 1,
103 1, 1, 1,-1,-1, 1,-1, 1,
104 -1, 1, 1, 1, 1, 1, 1, 1,
105 1, 1, 1,-1, 1, 1, 1, 1,
106 1,-1,-1,-1,-1,-1,-1, 1,
107 -1, 1,-1, 1,-1,-1, 1,-1,
108 1, 1, 1, 1,-1, 1,-1, 1
109 };
110
111 /* Local Functions ----------------------------------------------------------*/
112
cnormf(complex float val)113 static float cnormf(complex float val) {
114 float realf = crealf(val);
115 float imagf = cimagf(val);
116
117 return realf * realf + imagf * imagf;
118 }
119
120 /*
121 * Gray coded QPSK modulation function
122 */
qpsk_mod(int * bits)123 complex float qpsk_mod(int *bits) {
124 return qpsk[(bits[1] << 1) | bits[0]];
125 }
126
127 /*
128 * Gray coded QPSK demodulation function
129 *
130 * 01 | 00
131 * ---+---
132 * 11 | 10
133 */
qpsk_demod(complex float symbol,int * bits)134 void qpsk_demod(complex float symbol, int *bits) {
135 complex float rotate = symbol * cmplx(ROT45);
136
137 bits[0] = crealf(rotate) <= 0.0f;
138 bits[1] = cimagf(rotate) <= 0.0f;
139 }
140
qam16_mod(int * bits)141 complex float qam16_mod(int *bits) {
142 return qam16[
143 (bits[3] << 3) | (bits[2] << 2) |
144 (bits[1] << 1) | bits[0]
145 ];
146 }
147
qam16_demod(complex float symbol,int * bits)148 void qam16_demod(complex float symbol, int *bits) {
149 float dist[16];
150 int i;
151
152 for (i = 0; i < 16; i++) {
153 dist[i] = cnormf(symbol - qam16[i]);
154 }
155
156 int row = 0;
157 float mdist = 10000.0f;
158
159 for (i = 0; i < 16; i++) {
160 if (dist[i] < mdist) {
161 mdist = dist[i];
162 row = i;
163 }
164 }
165
166 bits[0] = row & 1;
167 bits[1] = (row >> 1) & 1;
168 bits[2] = (row >> 2) & 1;
169 bits[3] = (row >> 3) & 1;
170 }
171
172 /*
173 * ------------
174 * ofdm_create
175 * ------------
176 *
177 * Returns OFDM data structure on success
178 * Return NULL on fail
179 *
180 * If you want the defaults, call this with config structure
181 * and the NC setting to 0. This will fill the structure with
182 * default values of the original OFDM modem.
183 */
ofdm_create(const struct OFDM_CONFIG * config)184 struct OFDM *ofdm_create(const struct OFDM_CONFIG *config) {
185 struct OFDM *ofdm;
186 float tval;
187 int i, j;
188
189 ofdm = (struct OFDM *) CALLOC(1, sizeof (struct OFDM));
190 assert(ofdm != NULL);
191
192 if (config == NULL) {
193 /* Fill in default values */
194
195 strcpy(ofdm->mode, "700D");
196 ofdm->nc = 17; /* Number of carriers */
197 ofdm->np = 1;
198 ofdm->ns = 8; /* Number of Symbols per modem frame */
199 ofdm->ts = 0.018f;
200 ofdm->rs = (1.0f / ofdm->ts); /* Modulation Symbol Rate */
201 ofdm->tcp = .002f; /* Cyclic Prefix duration */
202 ofdm->tx_centre = 1500.0f; /* TX Carrier Frequency */
203 ofdm->rx_centre = 1500.0f; /* RX Carrier Frequency */
204 ofdm->fs = 8000.0f; /* Sample rate */
205 ofdm->ntxtbits = 4;
206 ofdm->bps = 2; /* Bits per Symbol */
207 ofdm->nuwbits = 5 * ofdm->bps; /* default is 5 symbols of Unique Word bits */
208 ofdm->bad_uw_errors = 3;
209 ofdm->ftwindowwidth = 32;
210 ofdm->timing_mx_thresh = 0.30f;
211 ofdm->state_machine = "voice1";
212 ofdm->edge_pilots = 1;
213 ofdm->codename = "HRA_112_112";
214 ofdm->amp_est_mode = 0;
215 ofdm->tx_bpf_en = true;
216 ofdm->amp_scale = 245E3;
217 ofdm->clip_gain1 = 2.0;
218 ofdm->clip_gain2 = 0.9;
219 ofdm->clip_en = false;
220 ofdm->foff_limiter = false;
221 ofdm->data_mode = "";
222 memset(ofdm->tx_uw, 0, ofdm->nuwbits);
223 } else {
224 /* Use the users values */
225
226
227 strcpy(ofdm->mode, config->mode);
228 ofdm->nc = config->nc; /* Number of carriers */
229 ofdm->np = config->np; /* Number of modem Frames per Packet */
230 ofdm->ns = config->ns; /* Number of Symbol frames */
231 ofdm->bps = config->bps; /* Bits per Symbol */
232 ofdm->ts = config->ts;
233 ofdm->tcp = config->tcp; /* Cyclic Prefix duration */
234 ofdm->tx_centre = config->tx_centre; /* TX Centre Audio Frequency */
235 ofdm->rx_centre = config->rx_centre; /* RX Centre Audio Frequency */
236 ofdm->fs = config->fs; /* Sample Frequency */
237 ofdm->rs = config->rs; /* Symbol Rate */
238 ofdm->ntxtbits = config->txtbits;
239 ofdm->nuwbits = config->nuwbits;
240 ofdm->bad_uw_errors = config->bad_uw_errors;
241 ofdm->ftwindowwidth = config->ftwindowwidth;
242 ofdm->timing_mx_thresh = config->timing_mx_thresh;
243 ofdm->state_machine = config->state_machine;
244 ofdm->edge_pilots = config->edge_pilots;
245 ofdm->codename = config->codename;
246 ofdm->amp_est_mode = config->amp_est_mode;
247 ofdm->tx_bpf_en = config->tx_bpf_en;
248 ofdm->foff_limiter = config->foff_limiter;
249 ofdm->amp_scale = config->amp_scale;
250 ofdm->clip_gain1 = config->clip_gain1;
251 ofdm->clip_gain2 = config->clip_gain2;
252 ofdm->clip_en = config->clip_en;
253 memcpy(ofdm->tx_uw, config->tx_uw, ofdm->nuwbits);
254 ofdm->data_mode = config->data_mode;
255
256 }
257
258 ofdm->rs = (1.0f / ofdm->ts); /* Modulation Symbol Rate */
259 ofdm->m = (int) (ofdm->fs / ofdm->rs); /* 700D: 144 */
260 ofdm->ncp = (int) (ofdm->tcp * ofdm->fs); /* 700D: 16 */
261 ofdm->inv_m = (1.0f / (float) ofdm->m);
262
263 /* basic sanity checks */
264 assert((int)floorf(ofdm->fs / ofdm->rs) == ofdm->m);
265 assert(!strcmp(ofdm->state_machine, "voice1") ||
266 !strcmp(ofdm->state_machine, "data") ||
267 !strcmp(ofdm->state_machine, "voice2"));
268 assert(ofdm->nuwbits <= MAX_UW_BITS);
269
270 /* Copy constants into states */
271
272 strcpy(ofdm->config.mode, ofdm->mode);
273 ofdm->config.tx_centre = ofdm->tx_centre;
274 ofdm->config.rx_centre = ofdm->rx_centre;
275 ofdm->config.fs = ofdm->fs;
276 ofdm->config.rs = ofdm->rs;
277 ofdm->config.ts = ofdm->ts;
278 ofdm->config.tcp = ofdm->tcp;
279 ofdm->config.timing_mx_thresh = ofdm->timing_mx_thresh;
280 ofdm->config.nc = ofdm->nc;
281 ofdm->config.ns = ofdm->ns;
282 ofdm->config.np = ofdm->np;
283 ofdm->config.bps = ofdm->bps;
284 ofdm->config.nuwbits = ofdm->nuwbits;
285 ofdm->config.txtbits = ofdm->ntxtbits;
286 ofdm->config.bad_uw_errors = ofdm->bad_uw_errors;
287 ofdm->config.ftwindowwidth = ofdm->ftwindowwidth;
288 ofdm->config.state_machine = ofdm->state_machine;
289 ofdm->config.edge_pilots = ofdm->edge_pilots;
290 ofdm->config.codename = ofdm->codename;
291 ofdm->config.amp_est_mode = ofdm->amp_est_mode;
292 ofdm->config.tx_bpf_en = ofdm->tx_bpf_en;
293 ofdm->config.foff_limiter = ofdm->foff_limiter;
294 ofdm->config.amp_scale = ofdm->amp_scale;
295 ofdm->config.clip_gain1 = ofdm->clip_gain1;
296 ofdm->config.clip_gain2 = ofdm->clip_gain2;
297 ofdm->config.clip_en = ofdm->clip_en;
298 memcpy(ofdm->config.tx_uw, ofdm->tx_uw, ofdm->nuwbits);
299 ofdm->config.data_mode = ofdm->data_mode;
300
301 /* Calculate sizes from config param */
302
303 ofdm->bitsperframe = (ofdm->ns - 1) * (ofdm->nc * ofdm->bps); // 238 for nc = 17
304 ofdm->bitsperpacket = ofdm->np * ofdm->bitsperframe;
305 ofdm->tpacket = (float)(ofdm->np * ofdm->ns) * (ofdm->tcp + ofdm->ts); /* time for one packet */
306 ofdm->rowsperframe = ofdm->bitsperframe / (ofdm->nc * ofdm->bps);
307 ofdm->samplespersymbol = (ofdm->m + ofdm->ncp);
308 ofdm->samplesperframe = ofdm->ns * ofdm->samplespersymbol;
309 if (*ofdm->data_mode != 0)
310 // in burst data modes we skip ahead one frame to jump over preamble
311 ofdm->max_samplesperframe = 2*ofdm->samplesperframe;
312 else
313 ofdm->max_samplesperframe = ofdm->samplesperframe + (ofdm->samplespersymbol / 4);
314 /* extra storage at start of rxbuf to allow us to step back in time */
315 if (strlen(ofdm->data_mode))
316 ofdm->nrxbufhistory = (ofdm->np+2)*ofdm->samplesperframe;
317 else
318 ofdm->nrxbufhistory = 0;
319 ofdm->rxbufst = ofdm->nrxbufhistory;
320 ofdm->nrxbufmin = 3*ofdm->samplesperframe + 3*ofdm->samplespersymbol;
321 ofdm->nrxbuf = ofdm->nrxbufhistory + ofdm->nrxbufmin;
322
323 ofdm->pilot_samples = (complex float *) MALLOC(sizeof (complex float) * ofdm->samplespersymbol);
324 assert(ofdm->pilot_samples != NULL);
325
326 ofdm->rxbuf = (complex float *) MALLOC(sizeof (complex float) * ofdm->nrxbuf);
327 assert(ofdm->rxbuf != NULL);
328 for(int i=0; i<ofdm->nrxbuf; i++) ofdm->rxbuf[i] = 0;
329
330 ofdm->pilots = (complex float *) MALLOC(sizeof (complex float) * (ofdm->nc + 2));
331 assert(ofdm->pilots != NULL);
332
333 /*
334 * rx_sym is a 2D array of variable size
335 *
336 * allocate rx_sym row storage. It is a pointer to a pointer
337 */
338 ofdm->rx_sym = MALLOC(sizeof (complex float) * (ofdm->ns + 3));
339 assert(ofdm->rx_sym != NULL);
340
341 /* allocate rx_sym column storage */
342
343 for (i = 0; i < (ofdm->ns + 3); i++) {
344 ofdm->rx_sym[i] = (complex float *) MALLOC(sizeof(complex float) * (ofdm->nc + 2));
345 assert(ofdm->rx_sym[i] != NULL);
346 }
347
348 /* The rest of these are 1D arrays of variable size */
349
350 ofdm->rx_np = MALLOC(sizeof (complex float) * (ofdm->rowsperframe * ofdm->nc));
351 assert(ofdm->rx_np != NULL);
352
353 ofdm->rx_amp = MALLOC(sizeof (float) * (ofdm->rowsperframe * ofdm->nc));
354 assert(ofdm->rx_amp != NULL);
355
356 ofdm->aphase_est_pilot_log = MALLOC(sizeof (float) * (ofdm->rowsperframe * ofdm->nc));
357 assert(ofdm->aphase_est_pilot_log != NULL);
358
359 /* Null pointers to unallocated buffers */
360 ofdm->tx_bpf = NULL;
361 if (ofdm->tx_bpf_en)
362 allocate_tx_bpf(ofdm);
363
364 /* store complex BPSK pilot symbols */
365
366 assert(sizeof (pilotvalues) >= (ofdm->nc + 2) * sizeof (int8_t));
367
368 /* There are only 64 pilot values available */
369
370 for (i = 0; i < (ofdm->nc + 2); i++) {
371 ofdm->pilots[i] = ((float) pilotvalues[i]) + 0.0f * I;
372 }
373 if (ofdm->edge_pilots == 0) {
374 ofdm->pilots[0] = ofdm->pilots[ofdm->nc + 1] = 0.0f;
375 }
376 /* carrier tables for up and down conversion */
377
378 ofdm->doc = (TAU / (ofdm->fs / ofdm->rs));
379 tval = ((float) ofdm->nc / 2.0f);
380 ofdm->tx_nlower = roundf((ofdm->tx_centre / ofdm->rs) - tval) - 1.0f;
381 ofdm->rx_nlower = roundf((ofdm->rx_centre / ofdm->rs) - tval) - 1.0f;
382
383 for (i = 0; i < ofdm->nrxbuf; i++) {
384 ofdm->rxbuf[i] = 0.0f;
385 }
386
387 for (i = 0; i < (ofdm->ns + 3); i++) {
388 for (j = 0; j < (ofdm->nc + 2); j++) {
389 ofdm->rx_sym[i][j] = 0.0f;
390 }
391 }
392
393 for (i = 0; i < ofdm->rowsperframe * ofdm->nc; i++) {
394 ofdm->rx_np[i] = 0.0f;
395 }
396
397 for (i = 0; i < ofdm->rowsperframe; i++) {
398 for (j = 0; j < ofdm->nc; j++) {
399 ofdm->aphase_est_pilot_log[ofdm->nc * i + j] = 0.0f;
400 ofdm->rx_amp[ofdm->nc * i + j] = 0.0f;
401 }
402 }
403
404 /* default settings of options and states */
405
406 ofdm->verbose = 0;
407 ofdm->timing_en = true;
408 ofdm->foff_est_en = true;
409 ofdm->phase_est_en = true;
410 ofdm->phase_est_bandwidth = high_bw;
411 ofdm->phase_est_bandwidth_mode = AUTO_PHASE_EST;
412 ofdm->packetsperburst = 0; // default: never lose syn in raw data mode
413
414 ofdm->coarse_foff_est_hz = 0.0f;
415 ofdm->foff_est_gain = 0.1f;
416 ofdm->foff_est_hz = 0.0f;
417 ofdm->sample_point = 0;
418 ofdm->timing_est = 0;
419 ofdm->timing_valid = 0;
420 ofdm->timing_mx = 0.0f;
421 ofdm->nin = ofdm->samplesperframe;
422 ofdm->mean_amp = 0.0f;
423 ofdm->foff_metric = 0.0f;
424 /*
425 * Unique Word symbol placement. Note we need to group the UW
426 * bits so they fit into symbols. The LDPC decoder works on
427 * symbols so we can't break up any symbols into UW/payload bits.
428 */
429 ofdm->uw_ind = MALLOC(sizeof (int) * ofdm->nuwbits);
430 assert(ofdm->uw_ind != NULL);
431
432 ofdm->uw_ind_sym = MALLOC(sizeof (int) * (ofdm->nuwbits / ofdm->bps));
433 assert(ofdm->uw_ind_sym != NULL);
434
435 /*
436 * The Unique Word is placed in different indexes based on
437 * the number of carriers requested.
438 */
439 int nuwsyms = ofdm->nuwbits / ofdm->bps;
440 int Ndatasymsperframe = (ofdm->ns-1)*ofdm->nc;
441 int uw_step = ofdm->nc + 1; // default step size
442 int last_sym = floorf(nuwsyms*uw_step/ofdm->bps);
443 if (last_sym >= ofdm->np*Ndatasymsperframe)
444 uw_step = ofdm->nc - 1; // try a different step
445 last_sym = floorf(nuwsyms*uw_step/ofdm->bps);
446 assert(last_sym < ofdm->np*Ndatasymsperframe);// bail if we still can't fit them all
447
448 for (i = 0, j = 0; i < nuwsyms; i++, j += ofdm->bps) {
449 int val = floorf((i + 1) * uw_step / ofdm->bps);
450
451 ofdm->uw_ind_sym[i] = val; // symbol index
452
453 for (int b = 0; b < ofdm->bps ; b++) {
454 ofdm->uw_ind[j + b] = (val * ofdm->bps) + b;
455 }
456 }
457
458 // work out how many frames UW is spread over
459 int symsperframe = ofdm->bitsperframe / ofdm->bps;
460 ofdm->nuwframes = (int) ceilf((float)ofdm->uw_ind_sym[nuwsyms-1]/symsperframe);
461
462 ofdm->tx_uw_syms = MALLOC(sizeof (complex float) * (ofdm->nuwbits / ofdm->bps));
463 assert(ofdm->tx_uw_syms != NULL);
464
465 assert(ofdm->bps == 2); // TODO generalise
466 for (int s = 0; s < (ofdm->nuwbits / ofdm->bps); s++) {
467 int dibit[2];
468 dibit[1] = ofdm->tx_uw[2*s];
469 dibit[0] = ofdm->tx_uw[2*s+1];
470 ofdm->tx_uw_syms[s] = qpsk_mod(dibit);
471 }
472
473 /* sync state machine */
474
475 ofdm->sync_state = search;
476 ofdm->last_sync_state = search;
477
478 ofdm->uw_errors = 0;
479 ofdm->sync_counter = 0;
480 ofdm->frame_count = 0;
481 ofdm->sync_start = false;
482 ofdm->sync_end = false;
483 ofdm->sync_mode = autosync;
484 ofdm->modem_frame = 0;
485
486 /* create the OFDM pilot time-domain waveform */
487
488 complex float *temp = MALLOC(sizeof (complex float) * ofdm->m);
489 assert(temp != NULL);
490
491 idft(ofdm, temp, ofdm->pilots);
492
493 /*
494 * pilot_samples is 160 samples, but timing and freq offset est
495 * were found by experiment to work better without a cyclic
496 * prefix, so we uses zeroes instead.
497 */
498
499 /* zero out Cyclic Prefix (CP) time-domain values */
500
501 for (i = 0; i < ofdm->ncp; i++) {
502 ofdm->pilot_samples[i] = 0.0f;
503 }
504
505 /* Now copy the whole thing after the above */
506
507 for (i = ofdm->ncp, j = 0; j < ofdm->m; i++, j++) {
508 ofdm->pilot_samples[i] = temp[j];
509 }
510
511 FREE(temp);
512
513 /* calculate constant used to normalise timing correlation maximum */
514 float acc = 0.0f;
515 for (i = 0; i < ofdm->samplespersymbol; i++) {
516 acc += cnormf(ofdm->pilot_samples[i]);
517 }
518
519 ofdm->timing_norm = ofdm->samplespersymbol * acc;
520 ofdm->clock_offset_counter = 0;
521 ofdm->dpsk_en = false;
522
523 if (strlen(ofdm->data_mode)) {
524 ofdm->tx_preamble = (COMP*)malloc(sizeof(COMP)*ofdm->samplesperframe);
525 assert(ofdm->tx_preamble != NULL);
526 ofdm_generate_preamble(ofdm, ofdm->tx_preamble, 2);
527 ofdm->tx_postamble = (COMP*)malloc(sizeof(COMP)*ofdm->samplesperframe);
528 assert(ofdm->tx_postamble != NULL);
529 ofdm_generate_preamble(ofdm, ofdm->tx_postamble, 3);
530 }
531 ofdm->postambledetectoren = !strcmp(ofdm->data_mode,"burst");
532
533 return ofdm; /* Success */
534 }
535
allocate_tx_bpf(struct OFDM * ofdm)536 static void allocate_tx_bpf(struct OFDM *ofdm) {
537 ofdm->tx_bpf = MALLOC(sizeof(struct quisk_cfFilter));
538 assert(ofdm->tx_bpf != NULL);
539
540 /* Transmit bandpass filter; complex coefficients, center frequency */
541
542 if (!strcmp(ofdm->mode, "700D")) {
543 quisk_filt_cfInit(ofdm->tx_bpf, filtP650S900, sizeof (filtP650S900) / sizeof (float));
544 quisk_cfTune(ofdm->tx_bpf, ofdm->tx_centre / ofdm->fs);
545 }
546 else if (!strcmp(ofdm->mode, "700E")) {
547 quisk_filt_cfInit(ofdm->tx_bpf, filtP900S1100, sizeof (filtP900S1100) / sizeof (float));
548 quisk_cfTune(ofdm->tx_bpf, ofdm->tx_centre / ofdm->fs);
549 }
550 else if (!strcmp(ofdm->mode, "datac0") || !strcmp(ofdm->mode, "datac3")) {
551 quisk_filt_cfInit(ofdm->tx_bpf, filtP400S600, sizeof (filtP400S600) / sizeof (float));
552 quisk_cfTune(ofdm->tx_bpf, ofdm->tx_centre / ofdm->fs);
553 }
554 else assert(0);
555 }
556
deallocate_tx_bpf(struct OFDM * ofdm)557 static void deallocate_tx_bpf(struct OFDM *ofdm) {
558 assert(ofdm->tx_bpf != NULL);
559 quisk_filt_destroy(ofdm->tx_bpf);
560 FREE(ofdm->tx_bpf);
561 ofdm->tx_bpf = NULL;
562 }
563
ofdm_destroy(struct OFDM * ofdm)564 void ofdm_destroy(struct OFDM *ofdm) {
565 int i;
566
567 if (strlen(ofdm->data_mode)) {
568 free(ofdm->tx_preamble);
569 free(ofdm->tx_postamble);
570 }
571 if (ofdm->tx_bpf) {
572 deallocate_tx_bpf(ofdm);
573 }
574
575 FREE(ofdm->pilot_samples);
576 FREE(ofdm->rxbuf);
577 FREE(ofdm->pilots);
578
579 for (i = 0; i < (ofdm->ns + 3); i++) { /* 2D array */
580 FREE(ofdm->rx_sym[i]);
581 }
582
583 FREE(ofdm->rx_sym);
584 FREE(ofdm->rx_np);
585 FREE(ofdm->rx_amp);
586 FREE(ofdm->aphase_est_pilot_log);
587 FREE(ofdm->tx_uw_syms);
588 FREE(ofdm->uw_ind);
589 FREE(ofdm->uw_ind_sym);
590 FREE(ofdm);
591 }
592
593 /*
594 * Convert frequency domain into time domain
595 *
596 * This algorithm was optimized for speed
597 */
idft(struct OFDM * ofdm,complex float * result,complex float * vector)598 static void idft(struct OFDM *ofdm, complex float *result, complex float *vector) {
599 int row, col;
600
601 result[0] = 0.0f;
602
603 for (col = 0; col < (ofdm->nc + 2); col++) {
604 result[0] += vector[col]; // cexp(j0) == 1
605 }
606
607 result[0] *= ofdm->inv_m;
608
609 for (row = 1; row < ofdm->m; row++) {
610 complex float c = cmplx(ofdm->tx_nlower * ofdm->doc *row);
611 complex float delta = cmplx(ofdm->doc * row);
612
613 result[row] = 0.0f;
614
615 for (col = 0; col < (ofdm->nc + 2); col++) {
616 result[row] += (vector[col] * c);
617 c *= delta;
618 }
619
620 result[row] *= ofdm->inv_m;
621 }
622 }
623
624 /*
625 * Convert time domain into frequency domain
626 *
627 * This algorithm was optimized for speed
628 */
dft(struct OFDM * ofdm,complex float * result,complex float * vector)629 static void dft(struct OFDM *ofdm, complex float *result, complex float *vector) {
630 int row, col;
631
632 for (col = 0; col < (ofdm->nc + 2); col++) {
633 result[col] = vector[0]; // conj(cexp(j0)) == 1
634 }
635
636 for (col = 0; col < (ofdm->nc + 2); col++) {
637 float tval = (ofdm->rx_nlower + col) * ofdm->doc;
638 complex float c = cmplxconj(tval);
639 complex float delta = c;
640
641 for (row = 1; row < ofdm->m; row++) {
642 result[col] += (vector[row] * c);
643 c *= delta;
644 }
645 }
646 }
647
vector_sum(complex float * a,int num_elements)648 static complex float vector_sum(complex float *a, int num_elements) {
649 complex float sum = 0.0f;
650 int i;
651
652 for (i = 0; i < num_elements; i++) {
653 sum += a[i];
654 }
655
656 return sum;
657 }
658
659
660 /*
661 * Correlates the OFDM pilot symbol samples with a window of received
662 * samples to determine the most likely timing offset. Combines two
663 * frames pilots so we need at least Nsamperframe+M+Ncp samples in rx.
664 *
665 * Can be used for acquisition (coarse timing), and fine timing.
666 *
667 * Breaks when freq offset approaches +/- symbol rate (e.g
668 * +/- 25 Hz for 700D).
669 */
est_timing(struct OFDM * ofdm,complex float * rx,int length,int fcoarse,float * timing_mx,int * timing_valid,int step)670 static int est_timing(struct OFDM *ofdm, complex float *rx, int length,
671 int fcoarse, float *timing_mx, int *timing_valid, int step) {
672 complex float corr_st, corr_en;
673 int Ncorr = length - (ofdm->samplesperframe + ofdm->samplespersymbol);
674 float corr[Ncorr];
675 int i, j;
676 float acc = 0.0f;
677
678 for (i = 0; i < length; i++) {
679 acc += cnormf(rx[i]);
680 }
681
682 float av_level = 1.0f/(2.0f * sqrtf(ofdm->timing_norm * acc / length) + 1E-12f);
683
684 /* precompute the freq shift multiplied by pilot samples outside of main loop */
685
686 PROFILE_VAR(wvecpilot);
687 PROFILE_SAMPLE(wvecpilot);
688
689 complex float wvec_pilot[ofdm->samplespersymbol];
690
691 switch(fcoarse) {
692 case -40:
693 for (j = 0; j < ofdm->samplespersymbol; j++)
694 wvec_pilot[j] = conjf(ofdm_wval[j]*ofdm->pilot_samples[j]);
695 break;
696 case 0:
697 for (j = 0; j < ofdm->samplespersymbol; j++)
698 wvec_pilot[j] = conjf(ofdm->pilot_samples[j]);
699 break;
700 case 40:
701 for (j = 0; j < ofdm->samplespersymbol; j++)
702 wvec_pilot[j] = ofdm_wval[j]*conjf(ofdm->pilot_samples[j]);
703 break;
704 default:
705 assert(0);
706 }
707
708 /* use of __REAL__ provides a speed in increase of 10ms/frame during acquisition, however complex
709 is fast enough for real time operation */
710
711 #if defined(__EMBEDDED__) && defined(__REAL__)
712 float rx_real[length];
713 float wvec_pilot_real[ofdm->samplespersymbol];
714 float wvec_pilot_imag[ofdm->samplespersymbol];
715
716 for (i = 0; i < length; i++) {
717 rx_real[i] = crealf(rx[i]);
718 }
719
720 for (i = 0; i < ofdm->samplespersymbol; i++) {
721 wvec_pilot_real[i] = crealf(wvec_pilot[i]);
722 wvec_pilot_imag[i] = cimagf(wvec_pilot[i]);
723 }
724
725 #endif
726 PROFILE_SAMPLE_AND_LOG2(wvecpilot, " wvecpilot");
727 PROFILE_VAR(corr_start);
728 PROFILE_SAMPLE(corr_start);
729
730 for (i = 0; i < Ncorr; i += step) {
731 corr_st = 0.0f;
732 corr_en = 0.0f;
733
734 #ifdef __EMBEDDED__
735 #ifdef __REAL__
736 float re,im;
737
738 arm_dot_prod_f32(&rx_real[i], wvec_pilot_real, ofdm->samplespersymbol, &re);
739 arm_dot_prod_f32(&rx_real[i], wvec_pilot_imag, ofdm->samplespersymbol, &im);
740 corr_st = re + im * I;
741
742 arm_dot_prod_f32(&rx_real[i+ ofdm->samplesperframe], wvec_pilot_real, ofdm->samplespersymbol, &re);
743 arm_dot_prod_f32(&rx_real[i+ ofdm->samplesperframe], wvec_pilot_imag, ofdm->samplespersymbol, &im);
744 corr_en = re + im * I;
745 #else
746 float re,im;
747
748 arm_cmplx_dot_prod_f32((float*)&rx[i], (float*)wvec_pilot, ofdm->samplespersymbol, &re, &im);
749 corr_st = re + im * I;
750
751 arm_cmplx_dot_prod_f32((float*)&rx[i+ ofdm->samplesperframe], (float*)wvec_pilot, ofdm->samplespersymbol, &re, &im);
752 corr_en = re + im * I;
753 #endif
754 #else
755 for (j = 0; j < ofdm->samplespersymbol; j++) {
756 int ind = i + j;
757
758 corr_st = corr_st + (rx[ind ] * wvec_pilot[j]);
759 corr_en = corr_en + (rx[ind + ofdm->samplesperframe] * wvec_pilot[j]);
760 }
761 #endif // __EMBEDDED__
762 corr[i] = (cabsf(corr_st) + cabsf(corr_en)) * av_level;
763 }
764
765 PROFILE_SAMPLE_AND_LOG2(corr_start, " corr");
766
767 /* find the max magnitude and its index */
768
769 int timing_est = 0;
770 *timing_mx = 0.0f;
771
772 for (i = 0; i < Ncorr; i+=step) {
773 if (corr[i] > *timing_mx) {
774 *timing_mx = corr[i];
775 timing_est = i;
776 }
777 }
778
779 // only declare timing valid if there are enough samples in rxbuf to demodulate a frame
780 *timing_valid = (cabsf(rx[timing_est]) > 0.0) && (*timing_mx > ofdm->timing_mx_thresh);
781
782 if (ofdm->verbose > 2) {
783 fprintf(stderr, " av_level: %f max: %f timing_est: %d timing_valid: %d\n", (double) av_level,
784 (double) *timing_mx, timing_est, *timing_valid);
785 }
786
787 return timing_est;
788 }
789
790 /*
791 * Determines frequency offset at current timing estimate, used for
792 * coarse freq offset estimation during acquisition. Works up to +/-
793 * the symbol rate, e.g. +/- 25Hz for the FreeDV 700D configuration.
794 */
est_freq_offset_pilot_corr(struct OFDM * ofdm,complex float * rx,int timing_est,int fcoarse)795 static float est_freq_offset_pilot_corr(struct OFDM *ofdm, complex float *rx, int timing_est, int fcoarse) {
796 int st = -20; int en = 20; float foff_est = 0.0f; float Cabs_max = 0.0f;
797
798 /* precompute the freq shift multiplied by pilot samples outside of main loop */
799
800 complex float wvec_pilot[ofdm->samplespersymbol];
801 int j;
802
803 switch(fcoarse) {
804 case -40:
805 for (j = 0; j < ofdm->samplespersymbol; j++)
806 wvec_pilot[j] = conjf(ofdm_wval[j]*ofdm->pilot_samples[j]);
807 break;
808 case 0:
809 for (j = 0; j < ofdm->samplespersymbol; j++)
810 wvec_pilot[j] = conjf(ofdm->pilot_samples[j]);
811 break;
812 case 40:
813 for (j = 0; j < ofdm->samplespersymbol; j++)
814 wvec_pilot[j] = ofdm_wval[j]*conjf(ofdm->pilot_samples[j]);
815 break;
816 default:
817 assert(0);
818 }
819
820 // sample sum of DFT magnitude of correlated signals at each freq offset and look for peak
821 for (int f = st; f < en; f++) {
822 complex float corr_st = 0.0f;
823 complex float corr_en = 0.0f;
824 float tmp = TAU * f / ofdm->fs;
825 complex float delta = cmplxconj(tmp);
826 complex float w = cmplxconj(0.0f);
827 int i;
828
829 for (i = 0; i < ofdm->samplespersymbol; i++) {
830 // "mix" down (correlate) the pilot sequences from frame with 0 Hz offset pilot samples
831 complex float csam = wvec_pilot[i] * w;
832 int est = timing_est + i;
833
834 corr_st += rx[est ] * csam;
835 corr_en += rx[est + ofdm->samplesperframe] * csam;
836 w = w * delta;
837 }
838
839 float Cabs = cabsf(corr_st) + cabsf(corr_en);
840
841 if (Cabs > Cabs_max) {
842 Cabs_max = Cabs;
843 foff_est = f;
844 }
845 }
846
847 ofdm->foff_metric = 0.0f; // not used in this version of freq est algorithm
848
849 if (ofdm->verbose > 2) {
850 fprintf(stderr, "cabs_max: %f foff_est: %f\n", (double) Cabs_max, (double) foff_est);
851 }
852
853 return foff_est;
854 }
855
856
857 /*
858 * ----------------------------------------------
859 * ofdm_txframe - modulates one frame of symbols
860 * ----------------------------------------------
861 */
ofdm_txframe(struct OFDM * ofdm,complex float * tx,complex float * tx_sym_lin)862 void ofdm_txframe(struct OFDM *ofdm, complex float *tx, complex float *tx_sym_lin) {
863 complex float aframe[ofdm->np * ofdm->ns][ofdm->nc + 2];
864 complex float asymbol[ofdm->m];
865 complex float asymbol_cp[ofdm->samplespersymbol];
866 int i, j, k, m;
867
868 /* initialize aframe to complex zero */
869
870 for (i = 0; i < (ofdm->np * ofdm->ns); i++) {
871 for (j = 0; j < (ofdm->nc + 2); j++) {
872 aframe[i][j] = 0.0f;
873 }
874 }
875
876 /*
877 * Place symbols in multi-carrier frame with pilots
878 * This will place boundary values of complex zero around data
879 */
880 int s = 0;
881 for (int r = 0; r < ofdm->np*ofdm->ns; r++) {
882
883 if ((r % ofdm->ns) == 0) {
884 /* copy in a row of complex pilots to first row of each frame */
885 for (i = 0; i < (ofdm->nc + 2); i++) {
886 aframe[r][i] = ofdm->pilots[i];
887 }
888 }
889 else {
890 /* copy in the Nc complex data symbols with [0 Nc 0] or (Nc + 2) total */
891 for (j = 1; j < (ofdm->nc + 1); j++) {
892 aframe[r][j] = tx_sym_lin[s++];
893 if (ofdm->dpsk_en == true) {
894 aframe[r][j] *= aframe[r-1][j];
895 }
896 }
897 }
898 }
899
900 /* OFDM up-convert symbol by symbol so we can add CP */
901
902 for (i = 0, m = 0; i < (ofdm->np * ofdm->ns); i++, m += ofdm->samplespersymbol) {
903 idft(ofdm, asymbol, aframe[i]);
904
905 /* Copy the last Ncp samples to the front */
906
907 for (j = (ofdm->m - ofdm->ncp), k = 0; j < ofdm->m; j++, k++) {
908 asymbol_cp[k] = asymbol[j];
909 }
910
911 /* Now copy the all samples for this row after it */
912
913 for (j = ofdm->ncp, k = 0; k < ofdm->m; j++, k++) {
914 asymbol_cp[j] = asymbol[k];
915 }
916
917 /* Now move row to the tx output */
918
919 for (j = 0; j < ofdm->samplespersymbol; j++) {
920 tx[m + j] = asymbol_cp[j];
921 }
922 }
923
924 size_t samplesperpacket = ofdm->np*ofdm->samplesperframe;
925 ofdm_hilbert_clipper(ofdm, tx, samplesperpacket);
926 }
927
928
929 /* Scale Tx signal and optionally apply two stage Hilbert clipper to improve PAPR */
ofdm_hilbert_clipper(struct OFDM * ofdm,complex float * tx,size_t n)930 void ofdm_hilbert_clipper(struct OFDM *ofdm, complex float *tx, size_t n) {
931
932 /* vanilla Tx output waveform should be about OFDM_PEAK */
933 for(int i=0; i<n; i++) tx[i] *= ofdm->amp_scale;
934
935 if (ofdm->clip_en) {
936 // this gain set the drive into the Hilbert Clipper and sets PAPR
937 for(int i=0; i<n; i++) tx[i] *= ofdm->clip_gain1;
938 ofdm_clip(tx, OFDM_PEAK, n);
939 }
940
941 /* BPF to remove out of band energy clipper introduces */
942 if (ofdm->tx_bpf_en) {
943 assert(!strcmp(ofdm->mode, "700D") || !strcmp(ofdm->mode, "700E")
944 || !strcmp(ofdm->mode, "datac0") || !strcmp(ofdm->mode, "datac3"));
945 assert(ofdm->tx_bpf != NULL);
946 complex float tx_filt[n];
947
948 quisk_ccfFilter(tx, tx_filt, n, ofdm->tx_bpf);
949 memmove(tx, tx_filt, n * sizeof (complex float));
950 }
951
952 /* BPF messes up peak levels, this gain gets back to approx OFDM_PEAK */
953 if (ofdm->tx_bpf_en && ofdm->clip_en)
954 for(int i=0; i<n; i++) tx[i] *= ofdm->clip_gain2;
955
956 /* a very small percentage of samples may still exceed OFDM_PEAK, in
957 clipped or unclipped mode. Lets remove them so we present consistent
958 levels to the transmitter */
959
960 ofdm_clip(tx, OFDM_PEAK, n);
961 }
962
963
ofdm_get_config_param(struct OFDM * ofdm)964 struct OFDM_CONFIG *ofdm_get_config_param(struct OFDM *ofdm) { return &ofdm->config; }
ofdm_get_nin(struct OFDM * ofdm)965 int ofdm_get_nin(struct OFDM *ofdm) {return ofdm->nin;}
ofdm_get_samples_per_frame(struct OFDM * ofdm)966 int ofdm_get_samples_per_frame(struct OFDM *ofdm) { return ofdm->samplesperframe;}
ofdm_get_samples_per_packet(struct OFDM * ofdm)967 int ofdm_get_samples_per_packet(struct OFDM *ofdm) { return ofdm->samplesperframe*ofdm->np;}
ofdm_get_max_samples_per_frame(struct OFDM * ofdm)968 int ofdm_get_max_samples_per_frame(struct OFDM *ofdm) {return ofdm->max_samplesperframe; }
ofdm_get_bits_per_frame(struct OFDM * ofdm)969 int ofdm_get_bits_per_frame(struct OFDM *ofdm) {return ofdm->bitsperframe; }
ofdm_get_bits_per_packet(struct OFDM * ofdm)970 int ofdm_get_bits_per_packet(struct OFDM *ofdm) {return ofdm->bitsperpacket; }
ofdm_set_verbose(struct OFDM * ofdm,int level)971 void ofdm_set_verbose(struct OFDM *ofdm, int level) { ofdm->verbose = level; }
972
ofdm_set_timing_enable(struct OFDM * ofdm,bool val)973 void ofdm_set_timing_enable(struct OFDM *ofdm, bool val) {
974 ofdm->timing_en = val;
975
976 if (ofdm->timing_en == false) {
977 /* manually set ideal timing instant */
978
979 ofdm->sample_point = (ofdm->ncp - 1);
980 }
981 }
982
ofdm_get_phase_est_bandwidth_mode(struct OFDM * ofdm)983 int ofdm_get_phase_est_bandwidth_mode(struct OFDM *ofdm) {
984 return ofdm->phase_est_bandwidth_mode; /* int version of enum */
985 }
986
ofdm_set_phase_est_bandwidth_mode(struct OFDM * ofdm,int val)987 void ofdm_set_phase_est_bandwidth_mode(struct OFDM *ofdm, int val) {
988 assert((val == AUTO_PHASE_EST) || (val == LOCKED_PHASE_EST));
989 ofdm->phase_est_bandwidth_mode = val;
990 }
991
ofdm_set_foff_est_enable(struct OFDM * ofdm,bool val)992 void ofdm_set_foff_est_enable(struct OFDM *ofdm, bool val) {
993 ofdm->foff_est_en = val;
994 }
995
ofdm_set_phase_est_enable(struct OFDM * ofdm,bool val)996 void ofdm_set_phase_est_enable(struct OFDM *ofdm, bool val) {
997 ofdm->phase_est_en = val;
998 }
999
ofdm_set_off_est_hz(struct OFDM * ofdm,float val)1000 void ofdm_set_off_est_hz(struct OFDM *ofdm, float val) {
1001 ofdm->foff_est_hz = val;
1002 }
1003
ofdm_set_tx_bpf(struct OFDM * ofdm,bool val)1004 void ofdm_set_tx_bpf(struct OFDM *ofdm, bool val) {
1005 if (val == true) {
1006 if (ofdm->tx_bpf == NULL)
1007 allocate_tx_bpf(ofdm);
1008
1009 ofdm->tx_bpf_en = true;
1010 }
1011 else {
1012 if (ofdm->tx_bpf != NULL)
1013 deallocate_tx_bpf(ofdm);
1014
1015 ofdm->tx_bpf_en = false;
1016 }
1017 }
1018
ofdm_set_dpsk(struct OFDM * ofdm,bool val)1019 void ofdm_set_dpsk(struct OFDM *ofdm, bool val) {
1020 ofdm->dpsk_en = val;
1021 }
1022
1023 // select burst mode, and set packets per burst
ofdm_set_packets_per_burst(struct OFDM * ofdm,int packetsperburst)1024 void ofdm_set_packets_per_burst(struct OFDM *ofdm, int packetsperburst) {
1025 ofdm->data_mode = "burst";
1026 ofdm->packetsperburst = packetsperburst;
1027 ofdm->postambledetectoren = true;
1028 }
1029
1030 /*
1031 * --------------------------------------
1032 * ofdm_mod - modulates one frame of bits
1033 * --------------------------------------
1034 */
ofdm_mod(struct OFDM * ofdm,COMP * result,const int * tx_bits)1035 void ofdm_mod(struct OFDM *ofdm, COMP *result, const int *tx_bits) {
1036 int length = ofdm->bitsperpacket / ofdm->bps;
1037 complex float *tx = (complex float *) result; // complex has same memory layout
1038 complex float tx_sym_lin[length];
1039 int dibit[2];
1040 int s, i;
1041
1042 if (ofdm->bps == 1) {
1043 /* Here we will have Nbitsperpacket / 1 */
1044
1045 for (s = 0; s < length; s++) {
1046 tx_sym_lin[s] = (float) (2 * tx_bits[s] - 1);
1047 }
1048 } else if (ofdm->bps == 2) {
1049 /* Here we will have Nbitsperpacket / 2 */
1050
1051 for (s = 0, i = 0; i < length; s += 2, i++) {
1052 dibit[0] = tx_bits[s + 1] & 0x1;
1053 dibit[1] = tx_bits[s ] & 0x1;
1054
1055 tx_sym_lin[i] = qpsk_mod(dibit);
1056 }
1057 } /* else if (ofdm->bps == 3) { } TODO */
1058
1059 ofdm_txframe(ofdm, tx, tx_sym_lin);
1060 }
1061
1062 /*
1063 * ----------------------------------------------------------------------------------
1064 * ofdm_sync_search - attempts to find coarse sync parameters for modem initial sync
1065 * ----------------------------------------------------------------------------------
1066 */
1067
1068 /*
1069 * This is a wrapper to maintain the older functionality
1070 * with an array of COMPs as input
1071 */
ofdm_sync_search(struct OFDM * ofdm,COMP * rxbuf_in)1072 int ofdm_sync_search(struct OFDM *ofdm, COMP *rxbuf_in) {
1073 /*
1074 * insert latest input samples into rxbuf
1075 * so it is primed for when we have to call ofdm_demod()
1076 */
1077
1078 /* note can't use memcpy when src and dest overlap */
1079 memmove(&ofdm->rxbuf[0], &ofdm->rxbuf[ofdm->nin],
1080 (ofdm->nrxbuf - ofdm->nin) * sizeof (complex float));
1081 memmove(&ofdm->rxbuf[(ofdm->nrxbuf - ofdm->nin)],
1082 rxbuf_in, ofdm->nin * sizeof (complex float));
1083
1084 return(ofdm_sync_search_core(ofdm));
1085 }
1086
1087 /*
1088 * This is a wrapper to reduce memory allocated.
1089 * This works with ofdm_demod and freedv_api. Gain is not used here.
1090 */
ofdm_sync_search_shorts(struct OFDM * ofdm,short * rxbuf_in,float gain)1091 int ofdm_sync_search_shorts(struct OFDM *ofdm, short *rxbuf_in, float gain) {
1092 int i, j;
1093
1094 /* shift the buffer left based on nin */
1095
1096 memmove(&ofdm->rxbuf[0], &ofdm->rxbuf[ofdm->nin],
1097 (ofdm->nrxbuf - ofdm->nin) * sizeof (complex float));
1098
1099 /* insert latest input samples onto tail of rxbuf */
1100
1101 for (j = 0, i = (ofdm->nrxbuf - ofdm->nin); i < ofdm->nrxbuf; j++, i++) {
1102 ofdm->rxbuf[i] = ((float)rxbuf_in[j] / 32767.0f);
1103 }
1104
1105 return ofdm_sync_search_core(ofdm);
1106 }
1107
1108 /* Joint estimation of timing and freq used for burst data acquisition */
1109
est_timing_and_freq(struct OFDM * ofdm,int * t_est,float * foff_est,complex float * rx,int Nrx,complex float * known_samples,int Npsam,int tstep,float fmin,float fmax,float fstep)1110 static float est_timing_and_freq(struct OFDM *ofdm,
1111 int *t_est, float *foff_est,
1112 complex float *rx, int Nrx,
1113 complex float *known_samples, int Npsam,
1114 int tstep, float fmin, float fmax, float fstep) {
1115 int Ncorr = Nrx - Npsam + 1;
1116 float max_corr = 0;
1117 *t_est = 0; *foff_est = 0.0;
1118 for (float afcoarse=fmin; afcoarse<=fmax; afcoarse += fstep) {
1119 float w = TAU * afcoarse / ofdm->fs;
1120 complex float mvec[Npsam];
1121 for(int i=0; i<Npsam; i++) {
1122 complex float ph = cmplx(w*i);
1123 mvec[i] = known_samples[i]*ph;
1124 }
1125 for(int t=0; t<Ncorr; t+=tstep) {
1126 complex float corr = 0;
1127 for(int i=0; i<Npsam; i++)
1128 corr += rx[i+t]*conjf(mvec[i]);
1129 if (cabsf(corr) > max_corr) {
1130 max_corr = cabsf(corr);
1131 *t_est = t;
1132 *foff_est = afcoarse;
1133 }
1134 }
1135 }
1136
1137 /* obtain normalised real number for timing_mx */
1138 float mag1=0, mag2=0;
1139 for(int i=0; i<Npsam; i++) {
1140 mag1 += cabsf(known_samples[i]*conjf(known_samples[i]));
1141 mag2 += cabsf(rx[i+*t_est]*conjf(rx[i+*t_est]));
1142 }
1143 float timing_mx = max_corr*max_corr/(mag1*mag2+1E-12);
1144 if (ofdm->verbose > 2) {
1145 fprintf(stderr, " t_est: %4d timing:mx: %f foff_est: %f\n", *t_est, (double)timing_mx, (double)*foff_est);
1146 }
1147
1148 return timing_mx;
1149 }
1150
1151 /* Two stage burst mode acquisition */
1152
burst_acquisition_detector(struct OFDM * ofdm,complex float * rx,int n,complex float * known_sequence,int * ct_est,float * foff_est,float * timing_mx)1153 static void burst_acquisition_detector(struct OFDM *ofdm,
1154 complex float *rx, int n,
1155 complex float *known_sequence,
1156 int *ct_est, float *foff_est, float *timing_mx)
1157 {
1158
1159 float fmin, fmax, fstep;
1160 int tstep;
1161
1162 // initial search over coarse grid
1163 tstep = 4; fstep = 5; fmin = -50.0; fmax = 50.0;
1164 *timing_mx = est_timing_and_freq(ofdm, ct_est, foff_est,
1165 &rx[n], 2*ofdm->samplesperframe,
1166 known_sequence, ofdm->samplesperframe,
1167 tstep, fmin, fmax, fstep);
1168
1169 // refine estimate over finer grid
1170 fmin = *foff_est - ceilf(fstep/2.0); fmax = *foff_est + ceilf(fstep/2.0);
1171 int fine_st = n + *ct_est - tstep/2.0;
1172 *timing_mx = est_timing_and_freq(ofdm, ct_est, foff_est,
1173 &rx[fine_st], ofdm->samplesperframe + tstep,
1174 known_sequence, ofdm->samplesperframe,
1175 1, fmin, fmax, 1.0);
1176
1177 // refer ct_est to nominal start of frame rx[n]
1178 *ct_est += fine_st - n;
1179 }
1180
ofdm_sync_search_burst(struct OFDM * ofdm)1181 static int ofdm_sync_search_burst(struct OFDM *ofdm) {
1182
1183 int st = ofdm->rxbufst + ofdm->m + ofdm->ncp + ofdm->samplesperframe;
1184 char *pre_post = "";
1185
1186 int pre_ct_est; float pre_foff_est, pre_timing_mx;
1187 burst_acquisition_detector(ofdm, ofdm->rxbuf, st, (complex float*)ofdm->tx_preamble,
1188 &pre_ct_est, &pre_foff_est, &pre_timing_mx);
1189
1190 int post_ct_est; float post_foff_est, post_timing_mx;
1191 if (ofdm->postambledetectoren)
1192 burst_acquisition_detector(ofdm, ofdm->rxbuf, st, (complex float*)ofdm->tx_postamble,
1193 &post_ct_est, &post_foff_est, &post_timing_mx);
1194
1195 int ct_est; float foff_est, timing_mx;
1196 if (!ofdm->postambledetectoren || (pre_timing_mx > post_timing_mx)) {
1197 timing_mx = pre_timing_mx; ct_est = pre_ct_est; foff_est = pre_foff_est;
1198 pre_post = "pre";
1199 } else {
1200 timing_mx = post_timing_mx; ct_est = post_ct_est; foff_est = post_foff_est;
1201 pre_post = "post";
1202 }
1203
1204 int timing_valid = timing_mx > ofdm->timing_mx_thresh;
1205 if (timing_valid) {
1206 if (!strcmp(pre_post, "post")) {
1207 ofdm->post++;
1208 // we won't be need any new samples for a while ....
1209 ofdm->nin = 0;
1210 // backup to first modem frame in packet
1211 ofdm->rxbufst -= ofdm->np*ofdm->samplesperframe;
1212 ofdm->rxbufst += ct_est;
1213 } else {
1214 ofdm->pre++;
1215 // ct_est is start of preamble, so advance past that to start of first modem frame
1216 ofdm->nin = ofdm->samplesperframe + ct_est - 1;
1217 }
1218 } else
1219 ofdm->nin = ofdm->samplesperframe;
1220
1221 ofdm->ct_est = ct_est;
1222 ofdm->foff_est_hz = foff_est;
1223 ofdm->timing_mx = timing_mx;
1224 ofdm->timing_valid = timing_valid;
1225
1226 if (ofdm->verbose > 1) {
1227 fprintf(stderr, " ct_est: %4d nin: %4d mx: %3.2f foff_est: % 5.1f timing_valid: %d %4s\n",
1228 ct_est, ofdm->nin, (double)timing_mx, (double)foff_est, timing_valid, pre_post);
1229 }
1230
1231 return ofdm->timing_valid;
1232 }
1233
1234 /*
1235 * Attempts to find coarse sync parameters for modem initial sync (streaming mode)
1236 */
ofdm_sync_search_stream(struct OFDM * ofdm)1237 static int ofdm_sync_search_stream(struct OFDM *ofdm) {
1238 int act_est, afcoarse;
1239
1240 /* Attempt coarse timing estimate (i.e. detect start of frame) at a range of frequency offsets */
1241
1242 int st = ofdm->rxbufst + ofdm->samplesperframe + ofdm->samplespersymbol;
1243 int en = st + 2 * ofdm->samplesperframe + ofdm->samplespersymbol;
1244
1245 int fcoarse = 0;
1246 float atiming_mx, timing_mx = 0.0f;
1247 int ct_est = 0;
1248 int atiming_valid, timing_valid = 0;
1249
1250 PROFILE_VAR(timing_start);
1251 PROFILE_SAMPLE(timing_start);
1252
1253 for (afcoarse = -40; afcoarse <= 40; afcoarse += 40) {
1254 act_est = est_timing(ofdm, &ofdm->rxbuf[st], (en - st), afcoarse, &atiming_mx, &atiming_valid, 2);
1255
1256 if (atiming_mx > timing_mx) {
1257 ct_est = act_est;
1258 timing_mx = atiming_mx;
1259 fcoarse = afcoarse;
1260 timing_valid = atiming_valid;
1261 }
1262 }
1263
1264 PROFILE_SAMPLE_AND_LOG2(timing_start, " timing");
1265
1266 /* refine freq est within -/+ 20 Hz window */
1267
1268 PROFILE_VAR(freq_start);
1269 PROFILE_SAMPLE(freq_start);
1270
1271 ofdm->coarse_foff_est_hz = est_freq_offset_pilot_corr(ofdm, &ofdm->rxbuf[st], ct_est, fcoarse);
1272 ofdm->coarse_foff_est_hz += fcoarse;
1273
1274 PROFILE_SAMPLE_AND_LOG2(freq_start, " freq");
1275
1276 if (ofdm->verbose > 1) {
1277 fprintf(stderr, " ct_est: %4d foff_est: %4.1f timing_valid: %d timing_mx: %5.4f\n",
1278 ct_est, (double) ofdm->coarse_foff_est_hz, timing_valid,
1279 (double)timing_mx);
1280 }
1281
1282 ofdm->timing_valid = timing_valid;
1283 if (ofdm->timing_valid != 0) {
1284 /* potential candidate found .... */
1285
1286 /* calculate number of samples we need on next buffer to get into sync */
1287
1288 ofdm->nin = ct_est;
1289
1290 /* reset modem states */
1291
1292 ofdm->sample_point = ofdm->timing_est = 0;
1293 ofdm->foff_est_hz = ofdm->coarse_foff_est_hz;
1294 ofdm->timing_mx = timing_mx;
1295 } else {
1296 ofdm->nin = ofdm->samplesperframe;
1297 }
1298
1299 ofdm->timing_mx = timing_mx;
1300
1301 return ofdm->timing_valid;
1302 }
1303
ofdm_sync_search_core(struct OFDM * ofdm)1304 static int ofdm_sync_search_core(struct OFDM *ofdm) {
1305 if (!strcmp(ofdm->data_mode, "burst"))
1306 return ofdm_sync_search_burst(ofdm);
1307 else
1308 return ofdm_sync_search_stream(ofdm);
1309 }
1310
1311 /*
1312 * ------------------------------------------
1313 * ofdm_demod - Demodulates one frame of bits
1314 * ------------------------------------------
1315 */
1316
1317 /*
1318 * This is a wrapper to maintain the older functionality with an
1319 * array of COMPs as input
1320 */
ofdm_demod(struct OFDM * ofdm,int * rx_bits,COMP * rxbuf_in)1321 void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
1322 complex float *rx = (complex float *) &rxbuf_in[0]; // complex has same memory layout
1323 int i, j;
1324
1325 /* shift the buffer left based on nin */
1326 for (i = 0, j = ofdm->nin; i < (ofdm->nrxbuf - ofdm->nin); i++, j++) {
1327 ofdm->rxbuf[i] = ofdm->rxbuf[j];
1328 }
1329
1330 /* insert latest input samples onto tail of rxbuf */
1331 for (j = 0, i = (ofdm->nrxbuf - ofdm->nin); i < ofdm->nrxbuf; j++, i++) {
1332 ofdm->rxbuf[i] = rx[j];
1333 }
1334
1335 ofdm_demod_core(ofdm, rx_bits);
1336 }
1337
1338 /*
1339 * This is a wrapper with a new interface to reduce memory allocated.
1340 * This works with ofdm_demod and freedv_api. Gain is not used here.
1341 */
ofdm_demod_shorts(struct OFDM * ofdm,int * rx_bits,short * rxbuf_in,float gain)1342 void ofdm_demod_shorts(struct OFDM *ofdm, int *rx_bits, short *rxbuf_in, float gain) {
1343 int i, j;
1344
1345 /* shift the buffer left based on nin */
1346
1347 for (i = 0, j = ofdm->nin; i < (ofdm->nrxbuf - ofdm->nin); i++, j++) {
1348 ofdm->rxbuf[i] = ofdm->rxbuf[j];
1349 }
1350
1351 /* insert latest input samples onto tail of rxbuf */
1352
1353 for (j = 0, i = (ofdm->nrxbuf - ofdm->nin); i < ofdm->nrxbuf; j++, i++) {
1354 ofdm->rxbuf[i] = ((float)rxbuf_in[j] / 32767.0f);
1355 }
1356
1357 ofdm_demod_core(ofdm, rx_bits);
1358 }
1359
1360 /*
1361 * This is the rest of the function which expects that the data is
1362 * already in ofdm->rxbuf
1363 */
ofdm_demod_core(struct OFDM * ofdm,int * rx_bits)1364 static void ofdm_demod_core(struct OFDM *ofdm, int *rx_bits) {
1365 int prev_timing_est = ofdm->timing_est;
1366 int i, j, k, rr, st, en;
1367
1368 /*
1369 * get user and calculated freq offset
1370 */
1371 float woff_est = TAU * ofdm->foff_est_hz / ofdm->fs;
1372
1373 /* update timing estimate ---------------------------------------------- */
1374
1375 if (ofdm->timing_en == true) {
1376 /* update timing at start of every frame */
1377
1378 st = ofdm->rxbufst + ofdm->samplespersymbol + ofdm->samplesperframe - floorf(ofdm->ftwindowwidth / 2) + ofdm->timing_est;
1379 en = st + ofdm->samplesperframe - 1 + ofdm->samplespersymbol + ofdm->ftwindowwidth;
1380
1381 complex float work[(en - st)];
1382
1383 /*
1384 * Adjust for the frequency error by shifting the phase
1385 * using a conjugate multiply
1386 */
1387 for (j = 0, i = st; i < en; j++, i++) {
1388 work[j] = ofdm->rxbuf[i] * cmplxconj(woff_est * i);
1389 }
1390
1391 int ft_est = est_timing(ofdm, work, (en - st), 0.0f, &ofdm->timing_mx, &ofdm->timing_valid, 1);
1392
1393 ofdm->timing_est += ft_est - ceilf((float)ofdm->ftwindowwidth / 2) + 1;
1394
1395 if (ofdm->verbose > 2) {
1396 fprintf(stderr, " ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, ofdm->timing_est,
1397 ofdm->sample_point);
1398 }
1399
1400 /* Black magic to keep sample_point inside cyclic prefix. Or something like that. */
1401
1402 ofdm->sample_point = max(ofdm->timing_est + 4, ofdm->sample_point);
1403 ofdm->sample_point = min(ofdm->timing_est + ofdm->ncp-4, ofdm->sample_point);
1404 }
1405
1406 /*
1407 * Convert the time-domain samples to the frequency-domain using the rx_sym
1408 * data matrix. This will be Nc+2 carriers of 11 symbols.
1409 *
1410 * You will notice there are Nc+2 BPSK symbols for each pilot symbol, and
1411 * that there are Nc QPSK symbols for each data symbol.
1412 *
1413 * XXXXXXXXXXXXXXXXX <-- Timing Slip
1414 * PPPPPPPPPPPPPPPPPPP <-- Previous Frames Pilot
1415 * DDDDDDDDDDDDDDDDD
1416 * DDDDDDDDDDDDDDDDD
1417 * DDDDDDDDDDDDDDDDD
1418 * DDDDDDDDDDDDDDDDD Ignore these past data symbols
1419 * DDDDDDDDDDDDDDDDD
1420 * DDDDDDDDDDDDDDDDD
1421 * DDDDDDDDDDDDDDDDD
1422 * PPPPPPPPPPPPPPPPPPP <-- This Frames Pilot
1423 * DDDDDDDDDDDDDDDDD
1424 * DDDDDDDDDDDDDDDDD
1425 * DDDDDDDDDDDDDDDDD
1426 * DDDDDDDDDDDDDDDDD These are the current data symbols to be decoded
1427 * DDDDDDDDDDDDDDDDD
1428 * DDDDDDDDDDDDDDDDD
1429 * DDDDDDDDDDDDDDDDD
1430 * PPPPPPPPPPPPPPPPPPP <-- Next Frames Pilot
1431 * DDDDDDDDDDDDDDDDD
1432 * DDDDDDDDDDDDDDDDD
1433 * DDDDDDDDDDDDDDDDD
1434 * DDDDDDDDDDDDDDDDD Ignore these next data symbols
1435 * DDDDDDDDDDDDDDDDD
1436 * DDDDDDDDDDDDDDDDD
1437 * DDDDDDDDDDDDDDDDD
1438 * PPPPPPPPPPPPPPPPPPP <-- Future Frames Pilot
1439 * XXXXXXXXXXXXXXXXX <-- Timing Slip
1440 *
1441 * So this algorithm will have seven data symbols and four pilot symbols to process.
1442 * The average of the four pilot symbols is our phase estimation.
1443 */
1444 for (i = 0; i < (ofdm->ns + 3); i++) {
1445 for (j = 0; j < (ofdm->nc + 2); j++) {
1446 ofdm->rx_sym[i][j] = 0.0f;
1447 }
1448 }
1449
1450 /*
1451 * "Previous" pilot symbol is one modem frame above.
1452 */
1453 st = ofdm->rxbufst + ofdm->samplespersymbol + 1 + ofdm->sample_point;
1454 en = st + ofdm->m;
1455
1456 complex float work[ofdm->m];
1457
1458 /* down-convert at current timing instant------------------------------- */
1459
1460 for (k = 0, j = st; j < en; k++, j++) {
1461 work[k] = ofdm->rxbuf[j] * cmplxconj(woff_est * j);
1462 }
1463
1464 /*
1465 * Each symbol is of course ofdm->samplespersymbol samples long and
1466 * becomes Nc+2 carriers after DFT.
1467 *
1468 * We put this carrier pilot symbol at the top of our matrix:
1469 *
1470 * 1 .................. Nc+2
1471 *
1472 * +----------------------+
1473 * | Previous Pilot | rx_sym[0]
1474 * +----------------------+
1475 * | |
1476 *
1477 */
1478 dft(ofdm, ofdm->rx_sym[0], work);
1479
1480 /*
1481 * "This" pilot comes after the extra symbol allotted at the top, and after
1482 * the "previous" pilot and previous data symbols (let's call it, the previous
1483 * modem frame).
1484 *
1485 * So we will now be starting at "this" pilot symbol, and continuing to the
1486 * "next" pilot symbol.
1487 *
1488 * In this routine we also process the current data symbols.
1489 */
1490 for (rr = 0; rr < (ofdm->ns + 1); rr++) {
1491 st = ofdm->rxbufst + ofdm->samplespersymbol + ofdm->samplesperframe + (rr * ofdm->samplespersymbol) + 1 + ofdm->sample_point;
1492 en = st + ofdm->m;
1493
1494 /* down-convert at current timing instant---------------------------------- */
1495
1496 for (k = 0, j = st; j < en; k++, j++) {
1497 work[k] = ofdm->rxbuf[j] * cmplxconj(woff_est * j);
1498 }
1499
1500 /*
1501 * We put these Nc+2 carrier symbols into our matrix after the previous pilot:
1502 *
1503 * 1 .................. Nc+2
1504 * | Previous Pilot | rx_sym[0]
1505 * +----------------------+
1506 * | This Pilot | rx_sym[1]
1507 * +----------------------+
1508 * | Data | rx_sym[2]
1509 * +----------------------+
1510 * | Data | rx_sym[3]
1511 * +----------------------+
1512 * | Data | rx_sym[4]
1513 * +----------------------+
1514 * | Data | rx_sym[5]
1515 * +----------------------+
1516 * | Data | rx_sym[6]
1517 * +----------------------+
1518 * | Data | rx_sym[7]
1519 * +----------------------+
1520 * | Data | rx_sym[8]
1521 * +----------------------+
1522 * | Next Pilot | rx_sym[9]
1523 * +----------------------+
1524 * | | rx_sym[10]
1525 */
1526 dft(ofdm, ofdm->rx_sym[rr + 1], work);
1527 }
1528
1529 /*
1530 * OK, now we want to process to the "future" pilot symbol. This is after
1531 * the "next" modem frame.
1532 *
1533 * We are ignoring the data symbols between the "next" pilot and "future" pilot.
1534 * We only want the "future" pilot symbol, to perform the averaging of all pilots.
1535 */
1536 st = ofdm->rxbufst + ofdm->samplespersymbol + (3 * ofdm->samplesperframe) + 1 + ofdm->sample_point;
1537 en = st + ofdm->m;
1538
1539 /* down-convert at current timing instant------------------------------- */
1540
1541 for (k = 0, j = st; j < en; k++, j++) {
1542 work[k] = ofdm->rxbuf[j] * cmplxconj(woff_est * j);
1543 }
1544
1545 /*
1546 * We put the future pilot after all the previous symbols in the matrix:
1547 *
1548 * 1 .................. Nc+2
1549 *
1550 * | | rx_sym[9]
1551 * +----------------------+
1552 * | Future Pilot | rx_sym[10]
1553 * +----------------------+
1554 */
1555 dft(ofdm, ofdm->rx_sym[ofdm->ns + 2], work);
1556
1557 /*
1558 * We are finished now with the DFT and down conversion
1559 * From here on down we are in the frequency domain
1560 */
1561
1562 /* est freq err based on all carriers ---------------------------------- */
1563
1564 if (ofdm->foff_est_en == true) {
1565 /*
1566 * sym[1] is 'this' pilot symbol, sym[9] is 'next' pilot symbol.
1567 *
1568 * By subtracting the two averages of these pilots, we find the frequency
1569 * by the change in phase over time.
1570 */
1571 complex float freq_err_rect =
1572 conjf(vector_sum(ofdm->rx_sym[1], ofdm->nc + 2)) *
1573 vector_sum(ofdm->rx_sym[ofdm->ns + 1], ofdm->nc + 2);
1574
1575 /* prevent instability in atan(im/re) when real part near 0 */
1576
1577 freq_err_rect += 1E-6f;
1578
1579 float freq_err_hz = cargf(freq_err_rect) * ofdm->rs / (TAU * ofdm->ns);
1580 if (ofdm->foff_limiter) {
1581 /* optionally tame updates in low SNR channels */
1582 if (freq_err_hz > 1.0) freq_err_hz = 1.0;
1583 if (freq_err_hz < -1.0) freq_err_hz = -1.0;
1584 }
1585 ofdm->foff_est_hz += (ofdm->foff_est_gain * freq_err_hz);
1586 }
1587
1588 /* OK - now estimate and correct pilot phase -------------------------- */
1589
1590 complex float aphase_est_pilot_rect;
1591 float aphase_est_pilot[ofdm->nc + 2];
1592 float aamp_est_pilot[ofdm->nc + 2];
1593
1594 for (i = 0; i < (ofdm->nc + 2); i++) {
1595 aphase_est_pilot[i] = 10.0f;
1596 aamp_est_pilot[i] = 0.0f;
1597 }
1598
1599 for (i = 1; i < (ofdm->nc + 1); i++) { /* ignore first and last carrier for count */
1600 if (ofdm->phase_est_bandwidth == low_bw) {
1601 complex float symbol[3];
1602
1603 /*
1604 * Use all pilots normally, results in low SNR performance,
1605 * but will fall over in high Doppler propagation
1606 *
1607 * Basically we divide the Nc+2 pilots into groups of 3
1608 * Then average the phase surrounding each of the data symbols.
1609 */
1610 for (k = 0, j = (i - 1); k < 3; k++, j++) {
1611 symbol[k] = ofdm->rx_sym[1][j] * conjf(ofdm->pilots[j]); /* this pilot conjugate */
1612 }
1613
1614 aphase_est_pilot_rect = vector_sum(symbol, 3);
1615
1616 for (k = 0, j = (i - 1); k < 3; k++, j++) {
1617 symbol[k] = ofdm->rx_sym[ofdm->ns + 1][j] * conjf(ofdm->pilots[j]); /* next pilot conjugate */
1618 }
1619
1620 aphase_est_pilot_rect += vector_sum(symbol, 3);
1621
1622 /* use pilots in past and future */
1623
1624 for (k = 0, j = (i - 1); k < 3; k++, j++) {
1625 symbol[k] = ofdm->rx_sym[0][j] * conjf(ofdm->pilots[j]); /* previous pilot */
1626 }
1627
1628 aphase_est_pilot_rect += vector_sum(symbol, 3);
1629
1630 for (k = 0, j = (i - 1); k < 3; k++, j++) {
1631 symbol[k] = ofdm->rx_sym[ofdm->ns + 2][j] * conjf(ofdm->pilots[j]); /* future pilot */
1632 }
1633
1634 aphase_est_pilot_rect += vector_sum(symbol, 3);
1635
1636 /* amplitude is estimated over 12 pilots */
1637 aphase_est_pilot_rect /= 12.0f;
1638
1639 aphase_est_pilot[i] = cargf(aphase_est_pilot_rect);
1640 aamp_est_pilot[i] = cabsf(aphase_est_pilot_rect);
1641 } else {
1642 assert(ofdm->phase_est_bandwidth == high_bw);
1643
1644 /*
1645 * Use only symbols at 'this' and 'next' to quickly track changes
1646 * in phase due to high Doppler spread in propagation (no neighbor averaging).
1647 *
1648 * As less pilots are averaged, low SNR performance will be poorer
1649 */
1650 aphase_est_pilot_rect = ofdm->rx_sym[1][i] * conjf(ofdm->pilots[i]); /* "this" pilot conjugate */
1651 aphase_est_pilot_rect += ofdm->rx_sym[ofdm->ns + 1][i] * conjf(ofdm->pilots[i]); /* "next" pilot conjugate */
1652
1653 /* we estimate over 2 pilots */
1654 aphase_est_pilot_rect /= 2.0f;
1655 aphase_est_pilot[i] = cargf(aphase_est_pilot_rect);
1656
1657 if (ofdm->amp_est_mode == 0) {
1658 // legacy 700D ampl est method
1659 aamp_est_pilot[i] = cabsf(aphase_est_pilot_rect);
1660 } else {
1661 aamp_est_pilot[i] = cabsf(ofdm->rx_sym[1][i]) + cabsf(ofdm->rx_sym[ofdm->ns + 1][i])/2.0;
1662 }
1663 }
1664
1665 aphase_est_pilot[i] = cargf(aphase_est_pilot_rect);
1666 aamp_est_pilot[i] = cabsf(aphase_est_pilot_rect);
1667 }
1668
1669 /*
1670 * correct the phase offset using phase estimate, and demodulate
1671 * bits, separate loop as it runs across cols (carriers) to get
1672 * frame bit ordering correct
1673 */
1674 complex float rx_corr;
1675 int abit[2];
1676 int bit_index = 0;
1677 float sum_amp = 0.0f;
1678
1679 for (rr = 0; rr < ofdm->rowsperframe; rr++) {
1680 /*
1681 * Note the i starts with the second carrier, ends with Nc+1.
1682 * so we ignore the first and last carriers.
1683 *
1684 * Also note we are using sym[2..8] or the seven data symbols.
1685 */
1686 for (i = 1; i < (ofdm->nc + 1); i++) {
1687 if (ofdm->phase_est_en == true) {
1688 if (ofdm->dpsk_en == true) {
1689 /* differential detection, using pilot as reference at start of frame */
1690 rx_corr = ofdm->rx_sym[rr + 2][i] * cmplxconj(cargf(ofdm->rx_sym[rr + 1][i]));
1691 } else {
1692 /* regular coherent detection */
1693 rx_corr = ofdm->rx_sym[rr + 2][i] * cmplxconj(aphase_est_pilot[i]);
1694 }
1695 } else {
1696 rx_corr = ofdm->rx_sym[rr + 2][i];
1697 }
1698
1699 /*
1700 * Output complex data symbols after phase correction;
1701 * (_np = no pilots) the pilot symbols have been removed
1702 */
1703 ofdm->rx_np[(rr * ofdm->nc) + (i - 1)] = rx_corr;
1704
1705 /*
1706 * Note even though amp ests are the same for each col,
1707 * the FEC decoder likes to have one amplitude per symbol
1708 * so convenient to log them all
1709 */
1710 ofdm->rx_amp[(rr * ofdm->nc) + (i - 1)] = aamp_est_pilot[i];
1711 sum_amp += aamp_est_pilot[i];
1712
1713 /*
1714 * Note like amps in this implementation phase ests are the
1715 * same for each col, but we log them for each symbol anyway
1716 */
1717 ofdm->aphase_est_pilot_log[(rr * ofdm->nc) + (i - 1)] = aphase_est_pilot[i];
1718
1719 if (ofdm->bps == 1) {
1720 rx_bits[bit_index++] = crealf(rx_corr) > 0.0f;
1721 } else if (ofdm->bps == 2) {
1722 /*
1723 * Only one final task, decode what quadrant the phase
1724 * is in, and return the dibits
1725 */
1726 qpsk_demod(rx_corr, abit);
1727 rx_bits[bit_index++] = abit[1];
1728 rx_bits[bit_index++] = abit[0];
1729 }
1730 }
1731 }
1732
1733 /* update mean amplitude estimate for LDPC decoder scaling */
1734
1735 ofdm->mean_amp = 0.9f * ofdm->mean_amp + 0.1f * sum_amp / (ofdm->rowsperframe * ofdm->nc);
1736
1737 /* Adjust nin to take care of sample clock offset */
1738
1739 ofdm->nin = ofdm->samplesperframe;
1740
1741 if (ofdm->timing_en == true) {
1742 ofdm->clock_offset_counter += (prev_timing_est - ofdm->timing_est);
1743
1744 int thresh = ofdm->samplespersymbol / 8;
1745 int tshift = ofdm->samplespersymbol / 4;
1746
1747 if (ofdm->timing_est > thresh) {
1748 ofdm->nin = ofdm->samplesperframe + tshift;
1749 ofdm->timing_est -= tshift;
1750 ofdm->sample_point -= tshift;
1751 } else if (ofdm->timing_est < -thresh) {
1752 ofdm->nin = ofdm->samplesperframe - tshift;
1753 ofdm->timing_est += tshift;
1754 ofdm->sample_point += tshift;
1755 }
1756 }
1757
1758 // use internal rxbuf samples if they are available
1759 int rxbufst_next = ofdm->rxbufst + ofdm->nin;
1760 if (rxbufst_next + ofdm->nrxbufmin <= ofdm->nrxbuf) {
1761 ofdm->rxbufst = rxbufst_next;
1762 ofdm->nin = 0;
1763 }
1764 }
1765
1766
1767 /*
1768 * Returns an estimate of Es/No in dB - see esno_est.m for more info
1769 */
ofdm_esno_est_calc(complex float * rx_sym,int nsym)1770 float ofdm_esno_est_calc(complex float *rx_sym, int nsym) {
1771
1772 float sig_var = 0;
1773 float step = 1.0f/nsym;
1774 for (int i = 0; i < nsym; i++)
1775 sig_var += (cnormf(rx_sym[i]) * step);
1776 float sig_rms = sqrtf(sig_var);
1777
1778 float sum_x = 0.0f; float sum_xx = 0.0f; int n = 0;
1779 for (int i = 0; i < nsym; i++) {
1780 complex float s = rx_sym[i];
1781
1782 if (cabsf(s) > sig_rms) {
1783 if (fabs(crealf(s)) > fabs(cimagf(s))) {
1784 sum_x += cimagf(s);
1785 sum_xx += cimagf(s) * cimagf(s);
1786 } else {
1787 sum_x += crealf(s);
1788 sum_xx += crealf(s) * crealf(s);
1789 }
1790 n++;
1791 }
1792 }
1793
1794 float noise_var;
1795 if (n > 1)
1796 noise_var = (n * sum_xx - sum_x * sum_x) / (n * (n - 1));
1797 else
1798 noise_var = sig_var;
1799 noise_var *= 2.0f;
1800
1801 float EsNodB = 10.0f * log10f((1E-12f + sig_var) / (1E-12f + noise_var));
1802 assert(isnan(EsNodB) == 0);
1803 return EsNodB;
1804 }
1805
1806
ofdm_snr_from_esno(struct OFDM * ofdm,float EsNodB)1807 float ofdm_snr_from_esno(struct OFDM *ofdm, float EsNodB) {
1808 float cyclic_power = 10.0f * log10f((float)(ofdm->ncp + ofdm->m) / ofdm->m);
1809 return EsNodB + 10.0f * log10f((float)(ofdm->nc * ofdm->rs) / 3000.0f) + cyclic_power;
1810 }
1811
1812 /*
1813 * state machine for 700D/2020
1814 */
ofdm_sync_state_machine_voice1(struct OFDM * ofdm,uint8_t * rx_uw)1815 void ofdm_sync_state_machine_voice1(struct OFDM *ofdm, uint8_t *rx_uw) {
1816 int i;
1817
1818 State next_state = ofdm->sync_state;
1819
1820 ofdm->sync_start = false;
1821 ofdm->sync_end = false;
1822
1823 if (ofdm->sync_state == search) {
1824 if (ofdm->timing_valid) {
1825 ofdm->frame_count = 0;
1826 ofdm->sync_counter = 0;
1827 ofdm->sync_start = true;
1828 ofdm->clock_offset_counter = 0;
1829 next_state = trial;
1830 }
1831 }
1832
1833 if ((ofdm->sync_state == synced) || (ofdm->sync_state == trial)) {
1834 ofdm->frame_count++;
1835
1836 /*
1837 * freq offset est may be too far out, and has aliases every 1/Ts, so
1838 * we use a Unique Word to get a really solid indication of sync.
1839 */
1840 ofdm->uw_errors = 0;
1841
1842 for (i = 0; i < ofdm->nuwbits; i++) {
1843 ofdm->uw_errors += ofdm->tx_uw[i] ^ rx_uw[i];
1844 }
1845
1846 /*
1847 * during trial sync we don't tolerate errors so much, we look
1848 * for 3 consecutive frames with low error rate to confirm sync
1849 */
1850 if (ofdm->sync_state == trial) {
1851 if (ofdm->uw_errors > 2) {
1852 /* if we exceed thresh stay in trial sync */
1853
1854 ofdm->sync_counter++;
1855 ofdm->frame_count = 0;
1856 }
1857
1858 if (ofdm->sync_counter == 2) {
1859 /* if we get two bad frames drop sync and start again */
1860
1861 next_state = search;
1862 ofdm->phase_est_bandwidth = high_bw;
1863 }
1864
1865 if (ofdm->frame_count == 4) {
1866 /* three good frames, sync is OK! */
1867
1868 next_state = synced;
1869 /* change to low bandwidth, but more accurate phase estimation */
1870 /* but only if not locked to high */
1871
1872 if (ofdm->phase_est_bandwidth_mode != LOCKED_PHASE_EST) {
1873 ofdm->phase_est_bandwidth = low_bw;
1874 }
1875 }
1876 }
1877
1878 /* once we have synced up we tolerate a higher error rate to wait out fades */
1879
1880 if (ofdm->sync_state == synced) {
1881 if (ofdm->uw_errors > 2) {
1882 ofdm->sync_counter++;
1883 } else {
1884 ofdm->sync_counter = 0;
1885 }
1886
1887 if ((ofdm->sync_mode == autosync) && (ofdm->sync_counter > 6)) {
1888 /* run of consecutive bad frames ... drop sync */
1889
1890 next_state = search;
1891 ofdm->phase_est_bandwidth = high_bw;
1892 }
1893 }
1894 }
1895
1896 ofdm->last_sync_state = ofdm->sync_state;
1897 ofdm->sync_state = next_state;
1898 }
1899
1900 /*
1901 * data (streaming mode) state machine
1902 */
ofdm_sync_state_machine_data_streaming(struct OFDM * ofdm,uint8_t * rx_uw)1903 void ofdm_sync_state_machine_data_streaming(struct OFDM *ofdm, uint8_t *rx_uw) {
1904 State next_state = ofdm->sync_state;
1905 int i;
1906
1907 ofdm->sync_start = ofdm->sync_end = 0;
1908
1909 if (ofdm->sync_state == search) {
1910 if (ofdm->timing_valid != 0) {
1911 ofdm->sync_start = true;
1912 ofdm->sync_counter = 0;
1913 next_state = trial;
1914 }
1915 }
1916
1917 ofdm->uw_errors = 0;
1918 for (i = 0; i < ofdm->nuwbits; i++) {
1919 ofdm->uw_errors += ofdm->tx_uw[i] ^ rx_uw[i];
1920 }
1921
1922 if (ofdm->sync_state == trial) {
1923 if (ofdm->uw_errors < ofdm->bad_uw_errors) {
1924 next_state = synced;
1925 ofdm->packet_count = 0;
1926 ofdm->modem_frame = ofdm->nuwframes;
1927 } else {
1928 ofdm->sync_counter++;
1929
1930 if (ofdm->sync_counter > ofdm->np) {
1931 next_state = search;
1932 }
1933 }
1934 }
1935
1936 // Note if frameperburst==0 we don't ever lose sync, which is useful for
1937 // stream based testing or external control of state machine
1938
1939 if (ofdm->sync_state == synced) {
1940 ofdm->modem_frame++;
1941
1942 if (ofdm->modem_frame >= ofdm->np) {
1943 ofdm->modem_frame = 0;
1944 ofdm->packet_count++;
1945 if (ofdm->packetsperburst) {
1946 if (ofdm->packet_count >= ofdm->packetsperburst)
1947 next_state = search;
1948 }
1949 }
1950
1951 }
1952
1953 ofdm->last_sync_state = ofdm->sync_state;
1954 ofdm->sync_state = next_state;
1955 }
1956
1957 /*
1958 * data (burst mode) state machine
1959 */
ofdm_sync_state_machine_data_burst(struct OFDM * ofdm,uint8_t * rx_uw)1960 void ofdm_sync_state_machine_data_burst(struct OFDM *ofdm, uint8_t *rx_uw) {
1961 State next_state = ofdm->sync_state;
1962 int i;
1963
1964 ofdm->sync_start = ofdm->sync_end = 0;
1965
1966 if (ofdm->sync_state == search) {
1967 if (ofdm->timing_valid != 0) {
1968 ofdm->sync_start = true;
1969 ofdm->sync_counter = 0;
1970 next_state = trial;
1971 }
1972 }
1973
1974 ofdm->uw_errors = 0;
1975 for (i = 0; i < ofdm->nuwbits; i++) {
1976 ofdm->uw_errors += ofdm->tx_uw[i] ^ rx_uw[i];
1977 }
1978
1979 /* pre or post-amble has told us this is the start of the packet. Confirm we
1980 have a valid frame by checking the UW after the modem frames containing
1981 the UW have been received */
1982 if (ofdm->sync_state == trial) {
1983 ofdm->sync_counter++;
1984 if (ofdm->sync_counter == ofdm->nuwframes) {
1985 if (ofdm->uw_errors < ofdm->bad_uw_errors) {
1986 next_state = synced;
1987 ofdm->packet_count = 0;
1988 ofdm->modem_frame = ofdm->nuwframes;
1989 } else {
1990 next_state = search;
1991 // reset rxbuf to make sure we only ever do a postamble loop once through same samples
1992 ofdm->rxbufst = ofdm->nrxbufhistory;
1993 for(int i=0; i<ofdm->nrxbuf; i++) ofdm->rxbuf[i] = 0;
1994 ofdm->uw_fails++;
1995 }
1996 }
1997 }
1998
1999 if (ofdm->sync_state == synced) {
2000 ofdm->modem_frame++;
2001 if (ofdm->modem_frame >= ofdm->np) {
2002 ofdm->modem_frame = 0;
2003 ofdm->packet_count++;
2004 if (ofdm->packetsperburst) {
2005 if (ofdm->packet_count >= ofdm->packetsperburst) {
2006 next_state = search;
2007 // reset rxbuf to make sure we only ever do a postamble loop once through same samples
2008 ofdm->rxbufst = ofdm->nrxbufhistory;
2009 for(int i=0; i<ofdm->nrxbuf; i++) ofdm->rxbuf[i] = 0;
2010 }
2011 }
2012 }
2013 }
2014
2015 ofdm->last_sync_state = ofdm->sync_state;
2016 ofdm->sync_state = next_state;
2017 }
2018
2019
ofdm_sync_state_machine_voice2(struct OFDM * ofdm,uint8_t * rx_uw)2020 void ofdm_sync_state_machine_voice2(struct OFDM *ofdm, uint8_t *rx_uw) {
2021 int i;
2022
2023 State next_state = ofdm->sync_state;
2024
2025 ofdm->sync_start = false;
2026 ofdm->sync_end = false;
2027
2028 if (ofdm->sync_state == search) {
2029 if (ofdm->timing_valid) {
2030 ofdm->frame_count = 0;
2031 ofdm->sync_counter = 0;
2032 ofdm->sync_start = true;
2033 ofdm->clock_offset_counter = 0;
2034 next_state = trial;
2035 }
2036 }
2037
2038 if ((ofdm->sync_state == synced) || (ofdm->sync_state == trial)) {
2039 ofdm->frame_count++;
2040
2041 ofdm->uw_errors = 0;
2042 for (i = 0; i < ofdm->nuwbits; i++) {
2043 ofdm->uw_errors += ofdm->tx_uw[i] ^ rx_uw[i];
2044 }
2045
2046 if (ofdm->sync_state == trial) {
2047 if (ofdm->uw_errors <= ofdm->bad_uw_errors)
2048 next_state = synced;
2049 else
2050 next_state = search;
2051 }
2052
2053 if (ofdm->sync_state == synced) {
2054 if (ofdm->uw_errors > ofdm->bad_uw_errors) {
2055 ofdm->sync_counter++;
2056 } else {
2057 ofdm->sync_counter = 0;
2058 }
2059
2060 if (ofdm->sync_counter == 6) {
2061 /* run of consecutive bad frames ... drop sync */
2062 next_state = search;
2063 }
2064 }
2065 }
2066
2067 ofdm->last_sync_state = ofdm->sync_state;
2068 ofdm->sync_state = next_state;
2069 }
2070
2071
2072 /* mode based dispatcher for sync state machines */
ofdm_sync_state_machine(struct OFDM * ofdm,uint8_t * rx_uw)2073 void ofdm_sync_state_machine(struct OFDM *ofdm, uint8_t *rx_uw) {
2074 if (!strcmp(ofdm->state_machine, "voice1"))
2075 ofdm_sync_state_machine_voice1(ofdm, rx_uw);
2076 if (!strcmp(ofdm->state_machine, "data")) {
2077 if (strcmp(ofdm->data_mode,"streaming") == 0)
2078 ofdm_sync_state_machine_data_streaming(ofdm, rx_uw);
2079 else
2080 ofdm_sync_state_machine_data_burst(ofdm, rx_uw);
2081 }
2082 if (!strcmp(ofdm->state_machine, "voice2"))
2083 ofdm_sync_state_machine_voice2(ofdm, rx_uw);
2084 }
2085
2086
2087 /*---------------------------------------------------------------------------* \
2088
2089 FUNCTIONS...: ofdm_set_sync
2090 AUTHOR......: David Rowe
2091 DATE CREATED: May 2018
2092
2093 External control of sync state machine.
2094 Ensure this is called in the same thread as ofdm_sync_state_machine().
2095
2096 \*---------------------------------------------------------------------------*/
2097
ofdm_set_sync(struct OFDM * ofdm,int sync_cmd)2098 void ofdm_set_sync(struct OFDM *ofdm, int sync_cmd) {
2099 assert(ofdm != NULL);
2100
2101 switch (sync_cmd) {
2102 case UN_SYNC:
2103 /* force manual unsync, which will cause sync state machine to
2104 have re-attempt sync */
2105 ofdm->sync_state = search;
2106 /* clear rxbuf so we don't try to sync on any existing OFDM signals
2107 in buffer */
2108 for (int i = 0; i < ofdm->nrxbuf; i++) ofdm->rxbuf[i] = 0.0f;
2109 break;
2110 case AUTO_SYNC:
2111 /* normal operating mode - sync state machine decides when to unsync */
2112 ofdm->sync_mode = autosync;
2113 break;
2114 case MANUAL_SYNC:
2115 /*
2116 * allow sync state machine to sync, but not to unsync, the
2117 * operator will decide that manually
2118 */
2119 ofdm->sync_mode = manualsync;
2120 break;
2121 default:
2122 assert(0);
2123 }
2124 }
2125
2126 /*---------------------------------------------------------------------------*\
2127
2128 FUNCTION....: ofdm_get_demod_stats()
2129 AUTHOR......: David Rowe
2130 DATE CREATED: May 2018
2131
2132 Fills stats structure with a bunch of demod information. Call once per
2133 packet.
2134
2135 \*---------------------------------------------------------------------------*/
2136
ofdm_get_demod_stats(struct OFDM * ofdm,struct MODEM_STATS * stats,complex float * rx_syms,int Nsymsperpacket)2137 void ofdm_get_demod_stats(struct OFDM *ofdm, struct MODEM_STATS *stats, complex float *rx_syms, int Nsymsperpacket) {
2138 stats->Nc = ofdm->nc;
2139 assert(stats->Nc <= MODEM_STATS_NC_MAX);
2140
2141 float EsNodB = ofdm_esno_est_calc(rx_syms, Nsymsperpacket);
2142 float SNR3kdB = ofdm_snr_from_esno(ofdm, EsNodB);
2143
2144 if (strlen(ofdm->data_mode))
2145 /* no smoothing as we have a large number of symbols per packet */
2146 stats->snr_est = SNR3kdB;
2147 else {
2148 /* in voice modes we further smooth SNR est, fast attack, slow decay */
2149 if (SNR3kdB > stats->snr_est)
2150 stats->snr_est = SNR3kdB;
2151 else
2152 stats->snr_est = 0.9f * stats->snr_est + 0.1f * SNR3kdB;
2153 }
2154 stats->sync = ((ofdm->sync_state == synced) || (ofdm->sync_state == trial));
2155 stats->foff = ofdm->foff_est_hz;
2156 stats->rx_timing = ofdm->timing_est;
2157
2158 float total = ofdm->frame_count * ofdm->samplesperframe;
2159 stats->clock_offset = 0;
2160 if (total != 0.0f) {
2161 stats->clock_offset = ofdm->clock_offset_counter / total;
2162 }
2163
2164 stats->sync_metric = ofdm->timing_mx;
2165 stats->pre = ofdm->pre;
2166 stats->post = ofdm->post;
2167 stats->uw_fails = ofdm->uw_fails;
2168
2169 #ifndef __EMBEDDED__
2170 assert(Nsymsperpacket % ofdm->nc == 0);
2171 int Nrowsperpacket = Nsymsperpacket/ofdm->nc;
2172 assert(Nrowsperpacket <= MODEM_STATS_NR_MAX);
2173 stats->nr = Nrowsperpacket;
2174 for (int c = 0; c < ofdm->nc; c++) {
2175 for (int r = 0; r < Nrowsperpacket; r++) {
2176 complex float rot = rx_syms[r * ofdm->nc + c] * cmplx(ROT45);
2177 stats->rx_symbols[r][c].real = crealf(rot);
2178 stats->rx_symbols[r][c].imag = cimagf(rot);
2179 }
2180 }
2181 #endif
2182 }
2183
2184 /*
2185 * Assemble packet of bits from UW, payload bits, and txt bits
2186 */
ofdm_assemble_qpsk_modem_packet(struct OFDM * ofdm,uint8_t modem_frame[],uint8_t payload_bits[],uint8_t txt_bits[])2187 void ofdm_assemble_qpsk_modem_packet(struct OFDM *ofdm, uint8_t modem_frame[],
2188 uint8_t payload_bits[], uint8_t txt_bits[]) {
2189 int s, t;
2190
2191 int p = 0;
2192 int u = 0;
2193
2194 for (s = 0; s < (ofdm->bitsperpacket - ofdm->ntxtbits); s++) {
2195 if ((u < ofdm->nuwbits) && (s == ofdm->uw_ind[u])) {
2196 modem_frame[s] = ofdm->tx_uw[u++];
2197 } else {
2198 modem_frame[s] = payload_bits[p++];
2199 }
2200 }
2201
2202 assert(u == ofdm->nuwbits);
2203 assert(p == (ofdm->bitsperpacket - ofdm->nuwbits - ofdm->ntxtbits));
2204
2205 for (t = 0; s < ofdm->bitsperframe; s++, t++) {
2206 modem_frame[s] = txt_bits[t];
2207 }
2208
2209 assert(t == ofdm->ntxtbits);
2210 }
2211
2212 /*
2213 * Assemble packet of symbols from UW, payload symbols, and txt bits
2214 */
ofdm_assemble_qpsk_modem_packet_symbols(struct OFDM * ofdm,complex float modem_packet[],COMP payload_syms[],uint8_t txt_bits[])2215 void ofdm_assemble_qpsk_modem_packet_symbols(struct OFDM *ofdm, complex float modem_packet[],
2216 COMP payload_syms[], uint8_t txt_bits[]) {
2217 complex float *payload = (complex float *) &payload_syms[0]; // complex has same memory layout
2218 int Nsymsperpacket = ofdm->bitsperpacket / ofdm->bps;
2219 int Nuwsyms = ofdm->nuwbits / ofdm->bps;
2220 int Ntxtsyms = ofdm->ntxtbits / ofdm->bps;
2221 int dibit[2];
2222 int s, t;
2223
2224 int p = 0;
2225 int u = 0;
2226
2227 assert(ofdm->bps == 2); /* this only works for QPSK at this stage (e.g. modem packet mod) */
2228
2229 for (s = 0; s < (Nsymsperpacket - Ntxtsyms); s++) {
2230 if ((u < Nuwsyms) && (s == ofdm->uw_ind_sym[u])) {
2231 modem_packet[s] = ofdm->tx_uw_syms[u++];
2232 } else {
2233 modem_packet[s] = payload[p++];
2234 }
2235 }
2236
2237 assert(u == Nuwsyms);
2238 assert(p == (Nsymsperpacket - Nuwsyms - Ntxtsyms));
2239
2240 for (t = 0; s < Nsymsperpacket; s++, t += 2) {
2241 dibit[1] = txt_bits[t ] & 0x1;
2242 dibit[0] = txt_bits[t + 1] & 0x1;
2243 modem_packet[s] = qpsk_mod(dibit);
2244 }
2245
2246 assert(t == ofdm->ntxtbits);
2247 }
2248
2249 /*
2250 * Disassemble a received packet of symbols into UW bits and payload data symbols
2251 */
ofdm_disassemble_qpsk_modem_packet(struct OFDM * ofdm,complex float rx_syms[],float rx_amps[],COMP codeword_syms[],float codeword_amps[],short txt_bits[])2252 void ofdm_disassemble_qpsk_modem_packet(struct OFDM *ofdm, complex float rx_syms[], float rx_amps[],
2253 COMP codeword_syms[], float codeword_amps[], short txt_bits[])
2254 {
2255 complex float *codeword = (complex float *) &codeword_syms[0]; // complex has same memory layout
2256 int Nsymsperpacket = ofdm->bitsperpacket / ofdm->bps;
2257 int Nuwsyms = ofdm->nuwbits / ofdm->bps;
2258 int Ntxtsyms = ofdm->ntxtbits / ofdm->bps;
2259 int dibit[2];
2260 int s, t;
2261
2262 int p = 0;
2263 int u = 0;
2264
2265 assert(ofdm->bps == 2); /* this only works for QPSK at this stage */
2266
2267 for (s = 0; s < (Nsymsperpacket - Ntxtsyms); s++) {
2268 if ((u < Nuwsyms) && (s == ofdm->uw_ind_sym[u])) {
2269 u++;
2270 } else {
2271 codeword[p] = rx_syms[s];
2272 codeword_amps[p] = rx_amps[s];
2273 p++;
2274 }
2275 }
2276
2277 assert(u == Nuwsyms);
2278 assert(p == (Nsymsperpacket - Nuwsyms - Ntxtsyms));
2279
2280 for (t = 0; s < Nsymsperpacket; s++, t += 2) {
2281 qpsk_demod(rx_syms[s], dibit);
2282
2283 txt_bits[t ] = dibit[1];
2284 txt_bits[t + 1] = dibit[0];
2285 }
2286
2287 assert(t == ofdm->ntxtbits);
2288 }
2289
2290 /*
2291 * Disassemble a received packet of symbols into UW bits and payload data symbols
2292 */
ofdm_disassemble_qpsk_modem_packet_with_text_amps(struct OFDM * ofdm,complex float rx_syms[],float rx_amps[],COMP codeword_syms[],float codeword_amps[],short txt_bits[],int * textIndex)2293 void ofdm_disassemble_qpsk_modem_packet_with_text_amps(
2294 struct OFDM *ofdm, complex float rx_syms[], float rx_amps[],
2295 COMP codeword_syms[], float codeword_amps[], short txt_bits[],
2296 int* textIndex)
2297 {
2298 complex float *codeword = (complex float *) &codeword_syms[0]; // complex has same memory layout
2299 int Nsymsperpacket = ofdm->bitsperpacket / ofdm->bps;
2300 int Nuwsyms = ofdm->nuwbits / ofdm->bps;
2301 int Ntxtsyms = ofdm->ntxtbits / ofdm->bps;
2302 int dibit[2];
2303 int s, t;
2304
2305 int p = 0;
2306 int u = 0;
2307
2308 assert(ofdm->bps == 2); /* this only works for QPSK at this stage */
2309 assert(textIndex != NULL);
2310
2311 for (s = 0; s < (Nsymsperpacket - Ntxtsyms); s++) {
2312 if ((u < Nuwsyms) && (s == ofdm->uw_ind_sym[u])) {
2313 u++;
2314 } else {
2315 codeword[p] = rx_syms[s];
2316 codeword_amps[p] = rx_amps[s];
2317 p++;
2318 }
2319 }
2320
2321 assert(u == Nuwsyms);
2322 assert(p == (Nsymsperpacket - Nuwsyms - Ntxtsyms));
2323
2324 *textIndex = s;
2325 for (t = 0; s < Nsymsperpacket; s++, t += 2) {
2326 qpsk_demod(rx_syms[s], dibit);
2327
2328 txt_bits[t ] = dibit[1];
2329 txt_bits[t + 1] = dibit[0];
2330 }
2331
2332 assert(t == ofdm->ntxtbits);
2333 }
2334
2335 /*
2336 * Extract just the UW from the packet
2337 */
ofdm_extract_uw(struct OFDM * ofdm,complex float rx_syms[],float rx_amps[],uint8_t rx_uw[])2338 void ofdm_extract_uw(struct OFDM *ofdm, complex float rx_syms[], float rx_amps[], uint8_t rx_uw[]) {
2339 int Nsymsperframe = ofdm->bitsperframe / ofdm->bps;
2340 int Nuwsyms = ofdm->nuwbits / ofdm->bps;
2341 int dibit[2];
2342 int s,u;
2343
2344 assert(ofdm->bps == 2); /* this only works for QPSK at this stage (e.g. UW demod) */
2345
2346 for (s = 0, u = 0; s < Nsymsperframe*ofdm->nuwframes; s++) {
2347 if ((u < Nuwsyms) && (s == ofdm->uw_ind_sym[u])) {
2348 qpsk_demod(rx_syms[s], dibit);
2349 rx_uw[2 * u ] = dibit[1];
2350 rx_uw[2 * u + 1] = dibit[0];
2351 u++;
2352 }
2353 }
2354
2355 assert(u == Nuwsyms);
2356 }
2357
2358 /*
2359 * Pseudo-random number generator that we can implement in C with
2360 * identical results to Octave. Returns an unsigned int between 0
2361 * and 32767. Used for generating test frames of various lengths.
2362 */
ofdm_rand(uint16_t r[],int n)2363 void ofdm_rand(uint16_t r[], int n) {
2364 ofdm_rand_seed(r, n, 1);
2365 }
2366
ofdm_rand_seed(uint16_t r[],int n,uint64_t seed)2367 void ofdm_rand_seed(uint16_t r[], int n, uint64_t seed) {
2368 for (int i = 0; i < n; i++) {
2369 seed = (1103515245l * seed + 12345) % 32768;
2370 r[i] = seed;
2371 }
2372 }
2373
ofdm_generate_payload_data_bits(uint8_t payload_data_bits[],int n)2374 void ofdm_generate_payload_data_bits(uint8_t payload_data_bits[], int n) {
2375 uint16_t r[n];
2376 int i;
2377
2378 ofdm_rand(r, n);
2379
2380 for (i = 0; i < n; i++) {
2381 payload_data_bits[i] = r[i] > 16384;
2382 }
2383 }
2384
ofdm_generate_preamble(struct OFDM * ofdm,COMP * tx_preamble,int seed)2385 void ofdm_generate_preamble(struct OFDM *ofdm, COMP *tx_preamble, int seed) {
2386 // need to modify bits per packet to set up pre-amble of a few modem frames in length
2387 struct OFDM ofdm_preamble;
2388 memcpy(&ofdm_preamble, ofdm, sizeof(struct OFDM));
2389 ofdm_preamble.np = 1;
2390 ofdm_preamble.bitsperpacket = ofdm_preamble.bitsperframe;
2391 uint16_t r[ofdm_preamble.bitsperpacket];
2392 ofdm_rand_seed(r, ofdm_preamble.bitsperpacket, seed);
2393 int preamble_bits[ofdm_preamble.bitsperpacket];
2394 for(int i=0; i<ofdm_preamble.bitsperpacket; i++)
2395 preamble_bits[i] = r[i] > 16384;
2396 // ensures the signal passes through hilbert clipper unchanged
2397 ofdm_preamble.amp_scale = 1.0; ofdm_preamble.tx_bpf_en = false;
2398 ofdm_mod(&ofdm_preamble, tx_preamble, preamble_bits);
2399 }
2400
ofdm_print_info(struct OFDM * ofdm)2401 void ofdm_print_info(struct OFDM *ofdm) {
2402 char *syncmode[] = {
2403 "unsync",
2404 "autosync",
2405 "manualsync"
2406 };
2407 char *phase_est_bandwidth_mode[] = {
2408 "auto",
2409 "locked_high"
2410 };
2411
2412 fprintf(stderr, "ofdm->tx_centre = %g\n", (double)ofdm->tx_centre);
2413 fprintf(stderr, "ofdm->rx_centre = %g\n", (double)ofdm->rx_centre);
2414 fprintf(stderr, "ofdm->fs = %g\n", (double)ofdm->fs);
2415 fprintf(stderr, "ofdm->ts = %g\n", (double)ofdm->ts);
2416 fprintf(stderr, "ofdm->rs = %g\n", (double)ofdm->rs);
2417 fprintf(stderr, "ofdm->tcp = %g\n", (double)ofdm->tcp);
2418 fprintf(stderr, "ofdm->inv_m = %g\n", (double)ofdm->inv_m);
2419 fprintf(stderr, "ofdm->tx_nlower = %g\n", (double)ofdm->tx_nlower);
2420 fprintf(stderr, "ofdm->rx_nlower = %g\n", (double)ofdm->rx_nlower);
2421 fprintf(stderr, "ofdm->doc = %g\n", (double)ofdm->doc);
2422 fprintf(stderr, "ofdm->timing_mx_thresh = %g\n", (double)ofdm->timing_mx_thresh);
2423 fprintf(stderr, "ofdm->nc = %d\n", ofdm->nc);
2424 fprintf(stderr, "ofdm->np = %d\n", ofdm->np);
2425 fprintf(stderr, "ofdm->ns = %d\n", ofdm->ns);
2426 fprintf(stderr, "ofdm->bps = %d\n", ofdm->bps);
2427 fprintf(stderr, "ofdm->m = %d\n", ofdm->m);
2428 fprintf(stderr, "ofdm->ncp = %d\n", ofdm->ncp);
2429 fprintf(stderr, "ofdm->ftwindowwidth = %d\n", ofdm->ftwindowwidth);
2430 fprintf(stderr, "ofdm->bitsperframe = %d\n", ofdm->bitsperframe);
2431 fprintf(stderr, "ofdm->bitsperpacket = %d\n", ofdm->bitsperpacket);
2432 fprintf(stderr, "ofdm->rowsperframe = %d\n", ofdm->rowsperframe);
2433 fprintf(stderr, "ofdm->samplespersymbol = %d\n", ofdm->samplespersymbol);
2434 fprintf(stderr, "ofdm->samplesperframe = %d\n", ofdm->samplesperframe);
2435 fprintf(stderr, "ofdm->max_samplesperframe = %d\n", ofdm->max_samplesperframe);
2436 fprintf(stderr, "ofdm->nrxbuf = %d\n", ofdm->nrxbuf);
2437 fprintf(stderr, "ofdm->ntxtbits = %d\n", ofdm->ntxtbits);
2438 fprintf(stderr, "ofdm->nuwbits = %d\n", ofdm->nuwbits);
2439 fprintf(stderr, "ofdm->foff_est_gain = %g\n", (double)ofdm->foff_est_gain);
2440 fprintf(stderr, "ofdm->foff_est_hz = %g\n", (double)ofdm->foff_est_hz);
2441 fprintf(stderr, "ofdm->timing_mx = %g\n", (double)ofdm->timing_mx);
2442 fprintf(stderr, "ofdm->coarse_foff_est_hz = %g\n", (double)ofdm->coarse_foff_est_hz);
2443 fprintf(stderr, "ofdm->timing_norm = %g\n", (double)ofdm->timing_norm);
2444 fprintf(stderr, "ofdm->mean_amp = %g\n", (double)ofdm->mean_amp);
2445 fprintf(stderr, "ofdm->clock_offset_counter = %d\n", ofdm->clock_offset_counter);
2446 fprintf(stderr, "ofdm->verbose = %d\n", ofdm->verbose);
2447 fprintf(stderr, "ofdm->sample_point = %d\n", ofdm->sample_point);
2448 fprintf(stderr, "ofdm->timing_est = %d\n", ofdm->timing_est);
2449 fprintf(stderr, "ofdm->timing_valid = %d\n", ofdm->timing_valid);
2450 fprintf(stderr, "ofdm->nin = %d\n", ofdm->nin);
2451 fprintf(stderr, "ofdm->uw_errors = %d\n", ofdm->uw_errors);
2452 fprintf(stderr, "ofdm->sync_counter = %d\n", ofdm->sync_counter);
2453 fprintf(stderr, "ofdm->frame_count = %d\n", ofdm->frame_count);
2454 fprintf(stderr, "ofdm->sync_start = %s\n", ofdm->sync_start ? "true" : "false");
2455 fprintf(stderr, "ofdm->sync_end = %s\n", ofdm->sync_end ? "true" : "false");
2456 fprintf(stderr, "ofdm->sync_mode = %s\n", syncmode[ofdm->sync_mode]);
2457 fprintf(stderr, "ofdm->timing_en = %s\n", ofdm->timing_en ? "true" : "false");
2458 fprintf(stderr, "ofdm->foff_est_en = %s\n", ofdm->foff_est_en ? "true" : "false");
2459 fprintf(stderr, "ofdm->phase_est_en = %s\n", ofdm->phase_est_en ? "true" : "false");
2460 fprintf(stderr, "ofdm->tx_bpf_en = %s\n", ofdm->tx_bpf_en ? "true" : "false");
2461 fprintf(stderr, "ofdm->dpsk_en = %s\n", ofdm->dpsk_en ? "true" : "false");
2462 fprintf(stderr, "ofdm->phase_est_bandwidth_mode = %s\n", phase_est_bandwidth_mode[ofdm->phase_est_bandwidth_mode]);
2463 }
2464
2465 // hilbert clipper
ofdm_clip(complex float tx[],float clip_thresh,int n)2466 void ofdm_clip(complex float tx[], float clip_thresh, int n) {
2467 complex float sam;
2468 float mag;
2469 int i;
2470
2471 for(i=0; i<n; i++) {
2472 sam = tx[i];
2473 mag = cabsf(sam);
2474 if (mag > clip_thresh) {
2475 sam *= clip_thresh/mag;
2476 }
2477 tx[i] = sam;
2478 }
2479 }
2480