1 /*
2  * Copyright (c) 2007 - 2017 Joseph Gaeddert
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to deal
6  * in the Software without restriction, including without limitation the rights
7  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8  * copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20  * THE SOFTWARE.
21  */
22 
23 //
24 // ofdmframesync.c
25 //
26 // OFDM frame synchronizer
27 //
28 
29 #include <stdlib.h>
30 #include <stdio.h>
31 #include <string.h>
32 #include <math.h>
33 #include <assert.h>
34 
35 #include "liquid.internal.h"
36 
37 #define DEBUG_OFDMFRAMESYNC             1
38 #define DEBUG_OFDMFRAMESYNC_PRINT       0
39 #define DEBUG_OFDMFRAMESYNC_FILENAME    "ofdmframesync_internal_debug.m"
40 #define DEBUG_OFDMFRAMESYNC_BUFFER_LEN  (2048)
41 
42 #define OFDMFRAMESYNC_ENABLE_SQUELCH    0
43 
44 struct ofdmframesync_s {
45     unsigned int M;         // number of subcarriers
46     unsigned int M2;        // number of subcarriers (divided by 2)
47     unsigned int cp_len;    // cyclic prefix length
48     unsigned char * p;      // subcarrier allocation (null, pilot, data)
49 
50     // constants
51     unsigned int M_null;    // number of null subcarriers
52     unsigned int M_pilot;   // number of pilot subcarriers
53     unsigned int M_data;    // number of data subcarriers
54     unsigned int M_S0;      // number of enabled subcarriers in S0
55     unsigned int M_S1;      // number of enabled subcarriers in S1
56 
57     // scaling factors
58     float g_data;           // data symbols gain
59     float g_S0;             // S0 training symbols gain
60     float g_S1;             // S1 training symbols gain
61 
62     // transform object
63     FFT_PLAN fft;           // ifft object
64     float complex * X;      // frequency-domain buffer
65     float complex * x;      // time-domain buffer
66     windowcf input_buffer;  // input sequence buffer
67 
68     // PLCP sequences
69     float complex * S0;     // short sequence (freq)
70     float complex * s0;     // short sequence (time)
71     float complex * S1;     // long sequence (freq)
72     float complex * s1;     // long sequence (time)
73 
74     // gain
75     float g0;               // nominal gain (coarse initial estimate)
76     float complex * G0a;    // complex subcarrier gain estimate, S0[a]
77     float complex * G0b;    // complex subcarrier gain estimate, S0[b]
78     float complex * G1;     // complex subcarrier gain estimate, S1
79     float complex * G;      // complex subcarrier gain estimate
80     float complex * B;      // subcarrier phase rotation due to backoff
81     float complex * R;      //
82 
83     // receiver state
84     enum {
85         OFDMFRAMESYNC_STATE_SEEKPLCP=0,   // seek initial PLCP
86         OFDMFRAMESYNC_STATE_PLCPSHORT0,   // seek first PLCP short sequence
87         OFDMFRAMESYNC_STATE_PLCPSHORT1,   // seek second PLCP short sequence
88         OFDMFRAMESYNC_STATE_PLCPLONG,     // seek PLCP long sequence
89         OFDMFRAMESYNC_STATE_RXSYMBOLS     // receive payload symbols
90     } state;
91 
92     // synchronizer objects
93     nco_crcf nco_rx;        // numerically-controlled oscillator
94     msequence ms_pilot;     // pilot sequence generator
95     float phi_prime;        // ...
96     float p1_prime;         // filtered pilot phase slope
97 
98 #if OFDMFRAMESYNC_ENABLE_SQUELCH
99     // coarse signal detection
100     float squelch_threshold;
101     int squelch_enabled;
102 #endif
103 
104     // timing
105     unsigned int timer;         // input sample timer
106     unsigned int num_symbols;   // symbol counter
107     unsigned int backoff;       // sample timing backoff
108     float complex s_hat_0;      // first S0 symbol metrics estimate
109     float complex s_hat_1;      // second S0 symbol metrics estimate
110 
111     // detection thresholds
112     float plcp_detect_thresh;   // plcp detection threshold, nominally 0.35
113     float plcp_sync_thresh;     // long symbol threshold, nominally 0.30
114 
115     // callback
116     ofdmframesync_callback callback;
117     void * userdata;
118 
119 #if DEBUG_OFDMFRAMESYNC
120     int debug_enabled;
121     int debug_objects_created;
122     windowcf debug_x;
123     windowf  debug_rssi;
124     windowcf debug_framesyms;
125     float complex * G_hat;  // complex subcarrier gain estimate, S1
126     float * px;             // pilot x-value
127     float * py;             // pilot y-value
128     float p_phase[2];       // pilot polyfit
129     windowf debug_pilot_0;  // pilot polyfit history, p[0]
130     windowf debug_pilot_1;  // pilot polyfit history, p[1]
131 #endif
132 };
133 
134 // create OFDM framing synchronizer object
135 //  _M          :   number of subcarriers, >10 typical
136 //  _cp_len     :   cyclic prefix length
137 //  _taper_len  :   taper length (OFDM symbol overlap)
138 //  _p          :   subcarrier allocation (null, pilot, data), [size: _M x 1]
139 //  _callback   :   user-defined callback function
140 //  _userdata   :   user-defined data pointer
ofdmframesync_create(unsigned int _M,unsigned int _cp_len,unsigned int _taper_len,unsigned char * _p,ofdmframesync_callback _callback,void * _userdata)141 ofdmframesync ofdmframesync_create(unsigned int           _M,
142                                    unsigned int           _cp_len,
143                                    unsigned int           _taper_len,
144                                    unsigned char *        _p,
145                                    ofdmframesync_callback _callback,
146                                    void *                 _userdata)
147 {
148     ofdmframesync q = (ofdmframesync) malloc(sizeof(struct ofdmframesync_s));
149 
150     // validate input
151     if (_M < 8) {
152         fprintf(stderr,"warning: ofdmframesync_create(), less than 8 subcarriers\n");
153     } else if (_M % 2) {
154         fprintf(stderr,"error: ofdmframesync_create(), number of subcarriers must be even\n");
155         exit(1);
156     } else if (_cp_len > _M) {
157         fprintf(stderr,"error: ofdmframesync_create(), cyclic prefix length cannot exceed number of subcarriers\n");
158         exit(1);
159     }
160     q->M = _M;
161     q->cp_len = _cp_len;
162 
163     // derived values
164     q->M2 = _M/2;
165 
166     // subcarrier allocation
167     q->p = (unsigned char*) malloc((q->M)*sizeof(unsigned char));
168     if (_p == NULL) {
169         ofdmframe_init_default_sctype(q->M, q->p);
170     } else {
171         memmove(q->p, _p, q->M*sizeof(unsigned char));
172     }
173 
174     // validate and count subcarrier allocation
175     ofdmframe_validate_sctype(q->p, q->M, &q->M_null, &q->M_pilot, &q->M_data);
176     if ( (q->M_pilot + q->M_data) == 0) {
177         fprintf(stderr,"error: ofdmframesync_create(), must have at least one enabled subcarrier\n");
178         exit(1);
179     } else if (q->M_data == 0) {
180         fprintf(stderr,"error: ofdmframesync_create(), must have at least one data subcarriers\n");
181         exit(1);
182     } else if (q->M_pilot < 2) {
183         fprintf(stderr,"error: ofdmframesync_create(), must have at least two pilot subcarriers\n");
184         exit(1);
185     }
186 
187     // create transform object
188     q->X = (float complex*) malloc((q->M)*sizeof(float complex));
189     q->x = (float complex*) malloc((q->M)*sizeof(float complex));
190     q->fft = FFT_CREATE_PLAN(q->M, q->x, q->X, FFT_DIR_FORWARD, FFT_METHOD);
191 
192     // create input buffer the length of the transform
193     q->input_buffer = windowcf_create(q->M + q->cp_len);
194 
195     // allocate memory for PLCP arrays
196     q->S0 = (float complex*) malloc((q->M)*sizeof(float complex));
197     q->s0 = (float complex*) malloc((q->M)*sizeof(float complex));
198     q->S1 = (float complex*) malloc((q->M)*sizeof(float complex));
199     q->s1 = (float complex*) malloc((q->M)*sizeof(float complex));
200     ofdmframe_init_S0(q->p, q->M, q->S0, q->s0, &q->M_S0);
201     ofdmframe_init_S1(q->p, q->M, q->S1, q->s1, &q->M_S1);
202 
203     // compute scaling factor
204     q->g_data = sqrtf(q->M) / sqrtf(q->M_pilot + q->M_data);
205     q->g_S0   = sqrtf(q->M) / sqrtf(q->M_S0);
206     q->g_S1   = sqrtf(q->M) / sqrtf(q->M_S1);
207 
208     // gain
209     q->g0 = 1.0f;
210     q->G0a = (float complex*) malloc((q->M)*sizeof(float complex));
211     q->G0b = (float complex*) malloc((q->M)*sizeof(float complex));
212     q->G   = (float complex*) malloc((q->M)*sizeof(float complex));
213     q->B   = (float complex*) malloc((q->M)*sizeof(float complex));
214     q->R   = (float complex*) malloc((q->M)*sizeof(float complex));
215 
216 #if 1
217     memset(q->G0a, 0x00, q->M*sizeof(float complex));
218     memset(q->G0b, 0x00, q->M*sizeof(float complex));
219     memset(q->G ,  0x00, q->M*sizeof(float complex));
220     memset(q->B,   0x00, q->M*sizeof(float complex));
221 #endif
222 
223     // timing backoff
224     q->backoff = q->cp_len < 2 ? q->cp_len : 2;
225     float phi = (float)(q->backoff)*2.0f*M_PI/(float)(q->M);
226     unsigned int i;
227     for (i=0; i<q->M; i++)
228         q->B[i] = liquid_cexpjf(i*phi);
229 
230     // set callback data
231     q->callback = _callback;
232     q->userdata = _userdata;
233 
234     //
235     // synchronizer objects
236     //
237 
238     // numerically-controlled oscillator
239     q->nco_rx = nco_crcf_create(LIQUID_NCO);
240 
241     // set pilot sequence
242     q->ms_pilot = msequence_create_default(8);
243 
244 #if OFDMFRAMESYNC_ENABLE_SQUELCH
245     // coarse detection
246     q->squelch_threshold = -25.0f;
247     q->squelch_enabled = 0;
248 #endif
249 
250     // reset object
251     ofdmframesync_reset(q);
252 
253 #if DEBUG_OFDMFRAMESYNC
254     q->debug_enabled = 0;
255     q->debug_objects_created = 0;
256 
257     q->debug_x =        NULL;
258     q->debug_rssi =     NULL;
259     q->debug_framesyms =NULL;
260 
261     q->G_hat = NULL;
262     q->px    = NULL;
263     q->py    = NULL;
264 
265     q->debug_pilot_0 = NULL;
266     q->debug_pilot_1 = NULL;
267 #endif
268 
269     // return object
270     return q;
271 }
272 
ofdmframesync_destroy(ofdmframesync _q)273 void ofdmframesync_destroy(ofdmframesync _q)
274 {
275 #if DEBUG_OFDMFRAMESYNC
276     // destroy debugging objects
277     if (_q->debug_x         != NULL) windowcf_destroy(_q->debug_x);
278     if (_q->debug_rssi      != NULL) windowf_destroy(_q->debug_rssi);
279     if (_q->debug_framesyms != NULL) windowcf_destroy(_q->debug_framesyms);
280     if (_q->G_hat           != NULL) free(_q->G_hat);
281     if (_q->px              != NULL) free(_q->px);
282     if (_q->py              != NULL) free(_q->py);
283     if (_q->debug_pilot_0   != NULL) windowf_destroy(_q->debug_pilot_0);
284     if (_q->debug_pilot_1   != NULL) windowf_destroy(_q->debug_pilot_1);
285 #endif
286 
287     // free subcarrier type array memory
288     free(_q->p);
289 
290     // free transform object
291     windowcf_destroy(_q->input_buffer);
292     free(_q->X);
293     free(_q->x);
294     FFT_DESTROY_PLAN(_q->fft);
295 
296     // clean up PLCP arrays
297     free(_q->S0);
298     free(_q->s0);
299     free(_q->S1);
300     free(_q->s1);
301 
302     // free gain arrays
303     free(_q->G0a);
304     free(_q->G0b);
305     free(_q->G);
306     free(_q->B);
307     free(_q->R);
308 
309     // destroy synchronizer objects
310     nco_crcf_destroy(_q->nco_rx);           // numerically-controlled oscillator
311     msequence_destroy(_q->ms_pilot);
312 
313     // free main object memory
314     free(_q);
315 }
316 
ofdmframesync_print(ofdmframesync _q)317 void ofdmframesync_print(ofdmframesync _q)
318 {
319     printf("ofdmframesync:\n");
320     printf("    num subcarriers     :   %-u\n", _q->M);
321     printf("    cyclic prefix len   :   %-u\n", _q->cp_len);
322     //printf("    taper len           :   %-u\n", _q->taper_len);
323 }
324 
ofdmframesync_reset(ofdmframesync _q)325 void ofdmframesync_reset(ofdmframesync _q)
326 {
327 #if 0
328     // reset gain parameters
329     unsigned int i;
330     for (i=0; i<_q->M; i++)
331         _q->G[i] = 1.0f;
332 #endif
333 
334     // reset synchronizer objects
335     nco_crcf_reset(_q->nco_rx);
336     msequence_reset(_q->ms_pilot);
337 
338     // reset timers
339     _q->timer = 0;
340     _q->num_symbols = 0;
341     _q->s_hat_0 = 0.0f;
342     _q->s_hat_1 = 0.0f;
343     _q->phi_prime = 0.0f;
344     _q->p1_prime = 0.0f;
345 
346     // set thresholds (increase for small number of subcarriers)
347     _q->plcp_detect_thresh = (_q->M > 44) ? 0.35f : 0.35f + 0.01f*(44 - _q->M);
348     _q->plcp_sync_thresh   = (_q->M > 44) ? 0.30f : 0.30f + 0.01f*(44 - _q->M);
349 
350     // reset state
351     _q->state = OFDMFRAMESYNC_STATE_SEEKPLCP;
352 }
353 
ofdmframesync_is_frame_open(ofdmframesync _q)354 int ofdmframesync_is_frame_open(ofdmframesync _q)
355 {
356     return (_q->state == OFDMFRAMESYNC_STATE_SEEKPLCP) ? 0 : 1;
357 }
358 
ofdmframesync_execute(ofdmframesync _q,float complex * _x,unsigned int _n)359 void ofdmframesync_execute(ofdmframesync _q,
360                            float complex * _x,
361                            unsigned int _n)
362 {
363     unsigned int i;
364     float complex x;
365     for (i=0; i<_n; i++) {
366         x = _x[i];
367 
368         // correct for carrier frequency offset
369         if (_q->state != OFDMFRAMESYNC_STATE_SEEKPLCP) {
370             nco_crcf_mix_down(_q->nco_rx, x, &x);
371             nco_crcf_step(_q->nco_rx);
372         }
373 
374         // save input sample to buffer
375         windowcf_push(_q->input_buffer,x);
376 
377 #if DEBUG_OFDMFRAMESYNC
378         if (_q->debug_enabled) {
379             windowcf_push(_q->debug_x, x);
380             windowf_push(_q->debug_rssi, crealf(x)*crealf(x) + cimagf(x)*cimagf(x));
381         }
382 #endif
383 
384         switch (_q->state) {
385         case OFDMFRAMESYNC_STATE_SEEKPLCP:
386             ofdmframesync_execute_seekplcp(_q);
387             break;
388         case OFDMFRAMESYNC_STATE_PLCPSHORT0:
389             ofdmframesync_execute_S0a(_q);
390             break;
391         case OFDMFRAMESYNC_STATE_PLCPSHORT1:
392             ofdmframesync_execute_S0b(_q);
393             break;
394         case OFDMFRAMESYNC_STATE_PLCPLONG:
395             ofdmframesync_execute_S1(_q);
396             break;
397         case OFDMFRAMESYNC_STATE_RXSYMBOLS:
398             ofdmframesync_execute_rxsymbols(_q);
399             break;
400         default:;
401         }
402 
403     } // for (i=0; i<_n; i++)
404 } // ofdmframesync_execute()
405 
406 // get receiver RSSI
ofdmframesync_get_rssi(ofdmframesync _q)407 float ofdmframesync_get_rssi(ofdmframesync _q)
408 {
409     return -10.0f*log10(_q->g0);
410 }
411 
412 // get receiver carrier frequency offset estimate
ofdmframesync_get_cfo(ofdmframesync _q)413 float ofdmframesync_get_cfo(ofdmframesync _q)
414 {
415     return nco_crcf_get_frequency(_q->nco_rx);
416 }
417 
418 // set receiver carrier frequency offset estimate
ofdmframesync_set_cfo(ofdmframesync _q,float _cfo)419 void ofdmframesync_set_cfo(ofdmframesync _q, float _cfo)
420 {
421     nco_crcf_set_frequency(_q->nco_rx, _cfo);
422 }
423 
424 //
425 // internal methods
426 //
427 
428 // frame detection
ofdmframesync_execute_seekplcp(ofdmframesync _q)429 void ofdmframesync_execute_seekplcp(ofdmframesync _q)
430 {
431     _q->timer++;
432 
433     if (_q->timer < _q->M)
434         return;
435 
436     // reset timer
437     _q->timer = 0;
438 
439     //
440     float complex * rc;
441     windowcf_read(_q->input_buffer, &rc);
442 
443     // estimate gain
444     unsigned int i;
445     float g = 0.0f;
446     for (i=_q->cp_len; i<_q->M + _q->cp_len; i++) {
447         // compute |rc[i]|^2 efficiently
448         g += crealf(rc[i])*crealf(rc[i]) + cimagf(rc[i])*cimagf(rc[i]);
449     }
450     g = (float)(_q->M) / g;
451 
452 #if OFDMFRAMESYNC_ENABLE_SQUELCH
453     // TODO : squelch here
454     if ( -10*log10f( sqrtf(g) ) < _q->squelch_threshold &&
455          _q->squelch_enabled)
456     {
457         printf("squelch\n");
458         return;
459     }
460 #endif
461 
462     // estimate S0 gain
463     ofdmframesync_estimate_gain_S0(_q, &rc[_q->cp_len], _q->G0a);
464 
465     float complex s_hat;
466     ofdmframesync_S0_metrics(_q, _q->G0a, &s_hat);
467     s_hat *= g;
468 
469     float tau_hat  = cargf(s_hat) * (float)(_q->M2) / (2*M_PI);
470 #if DEBUG_OFDMFRAMESYNC_PRINT
471     printf(" - gain=%12.3f, rssi=%12.8f, s_hat=%12.4f <%12.8f>, tau_hat=%8.3f\n",
472             sqrt(g),
473             -10*log10(g),
474             cabsf(s_hat), cargf(s_hat),
475             tau_hat);
476 #endif
477 
478     // save gain (permits dynamic invocation of get_rssi() method)
479     _q->g0 = g;
480 
481     //
482     if (cabsf(s_hat) > _q->plcp_detect_thresh) {
483 
484         int dt = (int)roundf(tau_hat);
485         // set timer appropriately...
486         _q->timer = (_q->M + dt) % (_q->M2);
487         _q->timer += _q->M; // add delay to help ensure good S0 estimate
488         _q->state = OFDMFRAMESYNC_STATE_PLCPSHORT0;
489 
490 #if DEBUG_OFDMFRAMESYNC_PRINT
491         printf("********** frame detected! ************\n");
492         printf("    s_hat   :   %12.8f <%12.8f>\n", cabsf(s_hat), cargf(s_hat));
493         printf("  tau_hat   :   %12.8f\n", tau_hat);
494         printf("    dt      :   %12d\n", dt);
495         printf("    timer   :   %12u\n", _q->timer);
496 #endif
497         //printf("exiting prematurely\n");
498         //ofdmframesync_destroy(_q);
499         //exit(1);
500     }
501 
502 }
503 
504 // frame detection
ofdmframesync_execute_S0a(ofdmframesync _q)505 void ofdmframesync_execute_S0a(ofdmframesync _q)
506 {
507     //printf("t : %u\n", _q->timer);
508     _q->timer++;
509 
510     if (_q->timer < _q->M2)
511         return;
512 
513     // reset timer
514     _q->timer = 0;
515 
516     //
517     float complex * rc;
518     windowcf_read(_q->input_buffer, &rc);
519 
520     // TODO : re-estimate nominal gain
521 
522     // estimate S0 gain
523     ofdmframesync_estimate_gain_S0(_q, &rc[_q->cp_len], _q->G0a);
524 
525     float complex s_hat;
526     ofdmframesync_S0_metrics(_q, _q->G0a, &s_hat);
527     s_hat *= _q->g0;
528 
529     _q->s_hat_0 = s_hat;
530 
531 #if DEBUG_OFDMFRAMESYNC_PRINT
532     float tau_hat  = cargf(s_hat) * (float)(_q->M2) / (2*M_PI);
533     printf("********** S0[0] received ************\n");
534     printf("    s_hat   :   %12.8f <%12.8f>\n", cabsf(s_hat), cargf(s_hat));
535     printf("  tau_hat   :   %12.8f\n", tau_hat);
536 #endif
537 
538 #if 0
539     // TODO : also check for phase of s_hat (should be small)
540     if (cabsf(s_hat) < 0.3f) {
541         // false alarm
542 #if DEBUG_OFDMFRAMESYNC_PRINT
543         printf("false alarm S0[0]\n");
544 #endif
545         ofdmframesync_reset(_q);
546         return;
547     }
548 #endif
549     _q->state = OFDMFRAMESYNC_STATE_PLCPSHORT1;
550 }
551 
552 // frame detection
ofdmframesync_execute_S0b(ofdmframesync _q)553 void ofdmframesync_execute_S0b(ofdmframesync _q)
554 {
555     //printf("t = %u\n", _q->timer);
556     _q->timer++;
557 
558     if (_q->timer < _q->M2)
559         return;
560 
561     // reset timer
562     _q->timer = _q->M + _q->cp_len - _q->backoff;
563 
564     //
565     float complex * rc;
566     windowcf_read(_q->input_buffer, &rc);
567 
568     // estimate S0 gain
569     ofdmframesync_estimate_gain_S0(_q, &rc[_q->cp_len], _q->G0b);
570 
571     float complex s_hat;
572     ofdmframesync_S0_metrics(_q, _q->G0b, &s_hat);
573     s_hat *= _q->g0;
574 
575     _q->s_hat_1 = s_hat;
576 
577 #if DEBUG_OFDMFRAMESYNC_PRINT
578     float tau_hat  = cargf(s_hat) * (float)(_q->M2) / (2*M_PI);
579     printf("********** S0[1] received ************\n");
580     printf("    s_hat   :   %12.8f <%12.8f>\n", cabsf(s_hat), cargf(s_hat));
581     printf("  tau_hat   :   %12.8f\n", tau_hat);
582 
583     // new timing offset estimate
584     tau_hat  = cargf(_q->s_hat_0 + _q->s_hat_1) * (float)(_q->M2) / (2*M_PI);
585     printf("  tau_hat * :   %12.8f\n", tau_hat);
586 
587     printf("**********\n");
588 #endif
589 
590     // re-adjust timer accordingly
591     float tau_prime = cargf(_q->s_hat_0 + _q->s_hat_1) * (float)(_q->M2) / (2*M_PI);
592     _q->timer -= (int)roundf(tau_prime);
593 
594 #if 0
595     if (cabsf(s_hat) < 0.3f) {
596 #if DEBUG_OFDMFRAMESYNC_PRINT
597         printf("false alarm S0[1]\n");
598 #endif
599         // false alarm
600         ofdmframesync_reset(_q);
601         return;
602     }
603 #endif
604 
605     float complex g_hat = 0.0f;
606     unsigned int i;
607     for (i=0; i<_q->M; i++)
608         g_hat += _q->G0b[i] * conjf(_q->G0a[i]);
609 
610 #if 0
611     // compute carrier frequency offset estimate using freq. domain method
612     float nu_hat = 2.0f * cargf(g_hat) / (float)(_q->M);
613 #else
614     // compute carrier frequency offset estimate using ML method
615     float complex t0 = 0.0f;
616     for (i=0; i<_q->M2; i++) {
617         t0 += conjf(rc[i])       *       _q->s0[i] *
618                     rc[i+_q->M2] * conjf(_q->s0[i+_q->M2]);
619     }
620     float nu_hat = cargf(t0) / (float)(_q->M2);
621 #endif
622 
623 #if DEBUG_OFDMFRAMESYNC_PRINT
624     printf("   nu_hat   :   %12.8f\n", nu_hat);
625 #endif
626 
627     // set NCO frequency
628     nco_crcf_set_frequency(_q->nco_rx, nu_hat);
629 
630     _q->state = OFDMFRAMESYNC_STATE_PLCPLONG;
631 }
632 
ofdmframesync_execute_S1(ofdmframesync _q)633 void ofdmframesync_execute_S1(ofdmframesync _q)
634 {
635     _q->timer--;
636 
637     if (_q->timer > 0)
638         return;
639 
640     // increment number of symbols observed
641     _q->num_symbols++;
642 
643     // run fft
644     float complex * rc;
645     windowcf_read(_q->input_buffer, &rc);
646 
647     // estimate S1 gain
648     // TODO : add backoff in gain estimation
649     ofdmframesync_estimate_gain_S1(_q, &rc[_q->cp_len], _q->G);
650 
651     // compute detector output
652     float complex g_hat = 0.0f;
653     unsigned int i;
654     for (i=0; i<_q->M; i++) {
655         //g_hat += _q->G[(i+1+_q->M)%_q->M]*conjf(_q->G[(i+_q->M)%_q->M]);
656         g_hat += _q->G[(i+1)%_q->M]*conjf(_q->G[i]);
657     }
658     g_hat /= _q->M_S1; // normalize output
659     g_hat *= _q->g0;
660 
661     // rotate by complex phasor relative to timing backoff
662     g_hat *= liquid_cexpjf((float)(_q->backoff)*2.0f*M_PI/(float)(_q->M));
663 
664 #if DEBUG_OFDMFRAMESYNC_PRINT
665     printf("    g_hat   :   %12.4f <%12.8f>\n", cabsf(g_hat), cargf(g_hat));
666 #endif
667 
668     // check conditions for g_hat:
669     //  1. magnitude should be large (near unity) when aligned
670     //  2. phase should be very near zero (time aligned)
671     if (cabsf(g_hat) > _q->plcp_sync_thresh && fabsf(cargf(g_hat)) < 0.1f*M_PI ) {
672         //printf("    acquisition\n");
673         _q->state = OFDMFRAMESYNC_STATE_RXSYMBOLS;
674         // reset timer
675         _q->timer = _q->M + _q->cp_len + _q->backoff;
676         _q->num_symbols = 0;
677 
678         // normalize gain by subcarriers, apply timing backoff correction
679         float g = (float)(_q->M) / sqrtf(_q->M_pilot + _q->M_data);
680         for (i=0; i<_q->M; i++) {
681             _q->G[i] *= g;          // gain due to relative subcarrier allocation
682             _q->G[i] *= _q->B[i];   // timing backoff correction
683         }
684 
685 #if 0
686         // TODO : choose number of taps more appropriately
687         //unsigned int ntaps = _q->M / 4;
688         unsigned int ntaps = (_q->M < 8) ? 2 : 8;
689         // FIXME : this is by far the most computationally complex part of synchronization
690         ofdmframesync_estimate_eqgain(_q, ntaps);
691 #else
692         unsigned int poly_order = 4;
693         if (poly_order >= _q->M_pilot + _q->M_data)
694             poly_order = _q->M_pilot + _q->M_data - 1;
695         ofdmframesync_estimate_eqgain_poly(_q, poly_order);
696 #endif
697 
698 #if 1
699         // compute composite gain
700         unsigned int i;
701         for (i=0; i<_q->M; i++)
702             _q->R[i] = _q->B[i] / _q->G[i];
703 #endif
704 
705         return;
706 #if 0
707         printf("exiting prematurely\n");
708         ofdmframesync_destroy(_q);
709         exit(1);
710 #endif
711     }
712 
713     // check if we are stuck searching for the S1 symbol
714     if (_q->num_symbols == 16) {
715 #if DEBUG_OFDMFRAMESYNC_PRINT
716         printf("could not find S1 symbol. bailing...\n");
717 #endif
718         ofdmframesync_reset(_q);
719     }
720 
721     // 'reset' timer (wait another half symbol)
722     _q->timer = _q->M2;
723 }
724 
ofdmframesync_execute_rxsymbols(ofdmframesync _q)725 void ofdmframesync_execute_rxsymbols(ofdmframesync _q)
726 {
727     // wait for timeout
728     _q->timer--;
729 
730     if (_q->timer == 0) {
731 
732         // run fft
733         float complex * rc;
734         windowcf_read(_q->input_buffer, &rc);
735         memmove(_q->x, &rc[_q->cp_len-_q->backoff], (_q->M)*sizeof(float complex));
736         FFT_EXECUTE(_q->fft);
737 
738         // recover symbol in internal _q->X buffer
739         ofdmframesync_rxsymbol(_q);
740 
741 #if DEBUG_OFDMFRAMESYNC
742         if (_q->debug_enabled) {
743             unsigned int i;
744             for (i=0; i<_q->M; i++) {
745                 if (_q->p[i] == OFDMFRAME_SCTYPE_DATA)
746                     windowcf_push(_q->debug_framesyms, _q->X[i]);
747             }
748         }
749 #endif
750         // invoke callback
751         if (_q->callback != NULL) {
752             int retval = _q->callback(_q->X, _q->p, _q->M, _q->userdata);
753 
754             if (retval != 0)
755                 ofdmframesync_reset(_q);
756         }
757 
758         // reset timer
759         _q->timer = _q->M + _q->cp_len;
760     }
761 
762 }
763 
764 // compute S0 metrics
ofdmframesync_S0_metrics(ofdmframesync _q,float complex * _G,float complex * _s_hat)765 void ofdmframesync_S0_metrics(ofdmframesync _q,
766                               float complex * _G,
767                               float complex * _s_hat)
768 {
769     // timing, carrier offset correction
770     unsigned int i;
771     float complex s_hat = 0.0f;
772 
773     // compute timing estimate, accumulate phase difference across
774     // gains on subsequent pilot subcarriers (note that all the odd
775     // subcarriers are NULL)
776     for (i=0; i<_q->M; i+=2) {
777         s_hat += _G[(i+2)%_q->M]*conjf(_G[i]);
778     }
779     s_hat /= _q->M_S0; // normalize output
780 
781     // set output values
782     *_s_hat = s_hat;
783 }
784 
785 // estimate short sequence gain
786 //  _q      :   ofdmframesync object
787 //  _x      :   input array (time), [size: M x 1]
788 //  _G      :   output gain (freq)
ofdmframesync_estimate_gain_S0(ofdmframesync _q,float complex * _x,float complex * _G)789 void ofdmframesync_estimate_gain_S0(ofdmframesync   _q,
790                                     float complex * _x,
791                                     float complex * _G)
792 {
793     // move input array into fft input buffer
794     memmove(_q->x, _x, (_q->M)*sizeof(float complex));
795 
796     // compute fft, storing result into _q->X
797     FFT_EXECUTE(_q->fft);
798 
799     // compute gain, ignoring NULL subcarriers
800     unsigned int i;
801     float gain = sqrtf(_q->M_S0) / (float)(_q->M);
802 
803     for (i=0; i<_q->M; i++) {
804         if (_q->p[i] != OFDMFRAME_SCTYPE_NULL && (i%2)==0) {
805             // NOTE : if cabsf(_q->S0[i]) == 0 then we can multiply by conjugate
806             //        rather than compute division
807             //_G[i] = _q->X[i] / _q->S0[i];
808             _G[i] = _q->X[i] * conjf(_q->S0[i]);
809         } else {
810             _G[i] = 0.0f;
811         }
812 
813         // normalize gain
814         _G[i] *= gain;
815     }
816 }
817 
818 // estimate long sequence gain
819 //  _q      :   ofdmframesync object
820 //  _x      :   input array (time), [size: M x 1]
821 //  _G      :   output gain (freq)
ofdmframesync_estimate_gain_S1(ofdmframesync _q,float complex * _x,float complex * _G)822 void ofdmframesync_estimate_gain_S1(ofdmframesync _q,
823                                     float complex * _x,
824                                     float complex * _G)
825 {
826     // move input array into fft input buffer
827     memmove(_q->x, _x, (_q->M)*sizeof(float complex));
828 
829     // compute fft, storing result into _q->X
830     FFT_EXECUTE(_q->fft);
831 
832     // compute gain, ignoring NULL subcarriers
833     unsigned int i;
834     float gain = sqrtf(_q->M_S1) / (float)(_q->M);
835     for (i=0; i<_q->M; i++) {
836         if (_q->p[i] != OFDMFRAME_SCTYPE_NULL) {
837             // NOTE : if cabsf(_q->S1[i]) == 0 then we can multiply by conjugate
838             //        rather than compute division
839             //_G[i] = _q->X[i] / _q->S1[i];
840             _G[i] = _q->X[i] * conjf(_q->S1[i]);
841         } else {
842             _G[i] = 0.0f;
843         }
844 
845         // normalize gain
846         _G[i] *= gain;
847     }
848 }
849 
850 // estimate complex equalizer gain from G0 and G1
851 //  _q      :   ofdmframesync object
852 //  _ntaps  :   number of time-domain taps for smoothing
ofdmframesync_estimate_eqgain(ofdmframesync _q,unsigned int _ntaps)853 void ofdmframesync_estimate_eqgain(ofdmframesync _q,
854                                    unsigned int _ntaps)
855 {
856 #if DEBUG_OFDMFRAMESYNC
857     if (_q->debug_enabled) {
858         // copy pre-smoothed gain
859         memmove(_q->G_hat, _q->G, _q->M*sizeof(float complex));
860     }
861 #endif
862 
863     // validate input
864     if (_ntaps == 0 || _ntaps > _q->M) {
865         fprintf(stderr, "error: ofdmframesync_estimate_eqgain(), ntaps must be in [1,M]\n");
866         exit(1);
867     }
868 
869     unsigned int i;
870 
871     // generate smoothing window (fft of temporal window)
872     for (i=0; i<_q->M; i++)
873         _q->x[i] = (i < _ntaps) ? 1.0f : 0.0f;
874     FFT_EXECUTE(_q->fft);
875 
876     memmove(_q->G0a, _q->G, _q->M*sizeof(float complex));
877 
878     // smooth complex equalizer gains
879     for (i=0; i<_q->M; i++) {
880         // set gain to zero for null subcarriers
881         if (_q->p[i] == OFDMFRAME_SCTYPE_NULL) {
882             _q->G[i] = 0.0f;
883             continue;
884         }
885 
886         float complex w;
887         float complex w0 = 0.0f;
888         float complex G_hat = 0.0f;
889 
890         unsigned int j;
891         for (j=0; j<_q->M; j++) {
892             if (_q->p[j] == OFDMFRAME_SCTYPE_NULL) continue;
893 
894             // select window sample from array
895             w = _q->X[(i + _q->M - j) % _q->M];
896 
897             // accumulate gain
898             //G_hat += w * 0.5f * (_q->G0a[j] + _q->G0b[j]);
899             G_hat += w * _q->G0a[j];
900             w0 += w;
901         }
902 
903         // eliminate divide-by-zero issues
904         if (cabsf(w0) < 1e-4f) {
905             fprintf(stderr,"error: ofdmframesync_estimate_eqgain(), weighting factor is zero\n");
906             w0 = 1.0f;
907         }
908         _q->G[i] = G_hat / w0;
909     }
910 }
911 
912 // estimate complex equalizer gain from G0 and G1 using polynomial fit
913 //  _q      :   ofdmframesync object
914 //  _order  :   polynomial order
ofdmframesync_estimate_eqgain_poly(ofdmframesync _q,unsigned int _order)915 void ofdmframesync_estimate_eqgain_poly(ofdmframesync _q,
916                                         unsigned int _order)
917 {
918 #if DEBUG_OFDMFRAMESYNC
919     if (_q->debug_enabled) {
920         // copy pre-smoothed gain
921         memmove(_q->G_hat, _q->G, _q->M*sizeof(float complex));
922     }
923 #endif
924 
925     // polynomial interpolation
926     unsigned int i;
927     unsigned int N = _q->M_pilot + _q->M_data;
928     if (_order > N-1) _order = N-1;
929     if (_order > 10)  _order = 10;
930     float x_freq[N];
931     float y_abs[N];
932     float y_arg[N];
933     float p_abs[_order+1];
934     float p_arg[_order+1];
935 
936     unsigned int n=0;
937     unsigned int k;
938     for (i=0; i<_q->M; i++) {
939 
940         // start at mid-point (effective fftshift)
941         k = (i + _q->M2) % _q->M;
942 
943         if (_q->p[k] != OFDMFRAME_SCTYPE_NULL) {
944             if (n == N) {
945                 fprintf(stderr, "error: ofdmframesync_estimate_eqgain_poly(), pilot subcarrier mismatch\n");
946                 exit(1);
947             }
948             // store resulting...
949             x_freq[n] = (k > _q->M2) ? (float)k - (float)(_q->M) : (float)k;
950             x_freq[n] = x_freq[n] / (float)(_q->M);
951             y_abs[n] = cabsf(_q->G[k]);
952             y_arg[n] = cargf(_q->G[k]);
953 
954             // update counter
955             n++;
956         }
957     }
958 
959     if (n != N) {
960         fprintf(stderr, "error: ofdmframesync_estimate_eqgain_poly(), pilot subcarrier mismatch\n");
961         exit(1);
962     }
963 
964     // try to unwrap phase
965     for (i=1; i<N; i++) {
966         while ((y_arg[i] - y_arg[i-1]) >  M_PI)
967             y_arg[i] -= 2*M_PI;
968         while ((y_arg[i] - y_arg[i-1]) < -M_PI)
969             y_arg[i] += 2*M_PI;
970     }
971 
972     // fit to polynomial
973     polyf_fit(x_freq, y_abs, N, p_abs, _order+1);
974     polyf_fit(x_freq, y_arg, N, p_arg, _order+1);
975 
976     // compute subcarrier gain
977     for (i=0; i<_q->M; i++) {
978         float freq = (i > _q->M2) ? (float)i - (float)(_q->M) : (float)i;
979         freq = freq / (float)(_q->M);
980         float A     = polyf_val(p_abs, _order+1, freq);
981         float theta = polyf_val(p_arg, _order+1, freq);
982         _q->G[i] = (_q->p[i] == OFDMFRAME_SCTYPE_NULL) ? 0.0f : A * liquid_cexpjf(theta);
983     }
984 
985 #if 0
986     for (i=0; i<N; i++)
987         printf("x(%3u) = %12.8f; y_abs(%3u) = %12.8f; y_arg(%3u) = %12.8f;\n",
988                 i+1, x_freq[i],
989                 i+1, y_abs[i],
990                 i+1, y_arg[i]);
991 
992     for (i=0; i<=_order; i++)
993         printf("p_abs(%3u) = %12.8f;\n", i+1, p_abs[i]);
994     for (i=0; i<=_order; i++)
995         printf("p_arg(%3u) = %12.8f;\n", i+1, p_arg[i]);
996 #endif
997 }
998 
999 // recover symbol, correcting for gain, pilot phase, etc.
ofdmframesync_rxsymbol(ofdmframesync _q)1000 void ofdmframesync_rxsymbol(ofdmframesync _q)
1001 {
1002     // apply gain
1003     unsigned int i;
1004     for (i=0; i<_q->M; i++)
1005         _q->X[i] *= _q->R[i];
1006 
1007     // polynomial curve-fit
1008     float x_phase[_q->M_pilot];
1009     float y_phase[_q->M_pilot];
1010     float p_phase[2];
1011 
1012     unsigned int n=0;
1013     unsigned int k;
1014     float complex pilot = 1.0f;
1015     for (i=0; i<_q->M; i++) {
1016 
1017         // start at mid-point (effective fftshift)
1018         k = (i + _q->M2) % _q->M;
1019 
1020         if (_q->p[k]==OFDMFRAME_SCTYPE_PILOT) {
1021             if (n == _q->M_pilot) {
1022                 fprintf(stderr,"warning: ofdmframesync_rxsymbol(), pilot subcarrier mismatch\n");
1023                 return;
1024             }
1025             pilot = (msequence_advance(_q->ms_pilot) ? 1.0f : -1.0f);
1026 #if 0
1027             printf("pilot[%3u] = %12.4e + j*%12.4e (expected %12.4e + j*%12.4e)\n",
1028                     k,
1029                     crealf(_q->X[k]), cimagf(_q->X[k]),
1030                     crealf(pilot),    cimagf(pilot));
1031 #endif
1032             // store resulting...
1033             x_phase[n] = (k > _q->M2) ? (float)k - (float)(_q->M) : (float)k;
1034             y_phase[n] = cargf(_q->X[k]*conjf(pilot));
1035 
1036             // update counter
1037             n++;
1038 
1039         }
1040     }
1041 
1042     if (n != _q->M_pilot) {
1043         fprintf(stderr,"warning: ofdmframesync_rxsymbol(), pilot subcarrier mismatch\n");
1044         return;
1045     }
1046 
1047     // try to unwrap phase
1048     for (i=1; i<_q->M_pilot; i++) {
1049         while ((y_phase[i] - y_phase[i-1]) >  M_PI)
1050             y_phase[i] -= 2*M_PI;
1051         while ((y_phase[i] - y_phase[i-1]) < -M_PI)
1052             y_phase[i] += 2*M_PI;
1053     }
1054 
1055     // fit phase to 1st-order polynomial (2 coefficients)
1056     polyf_fit(x_phase, y_phase, _q->M_pilot, p_phase, 2);
1057 
1058     // filter slope estimate (timing offset)
1059     float alpha = 0.3f;
1060     p_phase[1] = alpha*p_phase[1] + (1-alpha)*_q->p1_prime;
1061     _q->p1_prime = p_phase[1];
1062 
1063 #if DEBUG_OFDMFRAMESYNC
1064     if (_q->debug_enabled) {
1065         // save pilots
1066         memmove(_q->px, x_phase, _q->M_pilot*sizeof(float));
1067         memmove(_q->py, y_phase, _q->M_pilot*sizeof(float));
1068 
1069         // NOTE : swapping values for octave
1070         _q->p_phase[0] = p_phase[1];
1071         _q->p_phase[1] = p_phase[0];
1072 
1073         windowf_push(_q->debug_pilot_0, p_phase[0]);
1074         windowf_push(_q->debug_pilot_1, p_phase[1]);
1075     }
1076 #endif
1077 
1078     // compensate for phase offset
1079     // TODO : find more computationally efficient way to do this
1080     for (i=0; i<_q->M; i++) {
1081         // only apply to data/pilot subcarriers
1082         if (_q->p[i] == OFDMFRAME_SCTYPE_NULL) {
1083             _q->X[i] = 0.0f;
1084         } else {
1085             float fx    = (i > _q->M2) ? (float)i - (float)(_q->M) : (float)i;
1086             float theta = polyf_val(p_phase, 2, fx);
1087             _q->X[i] *= liquid_cexpjf(-theta);
1088         }
1089     }
1090 
1091     // adjust NCO frequency based on differential phase
1092     if (_q->num_symbols > 0) {
1093         // compute phase error (unwrapped)
1094         float dphi_prime = p_phase[0] - _q->phi_prime;
1095         while (dphi_prime >  M_PI) dphi_prime -= M_2_PI;
1096         while (dphi_prime < -M_PI) dphi_prime += M_2_PI;
1097 
1098         // adjust NCO proportionally to phase error
1099         nco_crcf_adjust_frequency(_q->nco_rx, 1e-3f*dphi_prime);
1100     }
1101     // set internal phase state
1102     _q->phi_prime = p_phase[0];
1103     //printf("%3u : theta : %12.8f, nco freq: %12.8f\n", _q->num_symbols, p_phase[0], nco_crcf_get_frequency(_q->nco_rx));
1104 
1105     // increment symbol counter
1106     _q->num_symbols++;
1107 
1108 #if 0
1109     for (i=0; i<_q->M_pilot; i++)
1110         printf("x_phase(%3u) = %12.8f; y_phase(%3u) = %12.8f;\n", i+1, x_phase[i], i+1, y_phase[i]);
1111     printf("poly : p0=%12.8f, p1=%12.8f\n", p_phase[0], p_phase[1]);
1112 #endif
1113 }
1114 
1115 // enable debugging
ofdmframesync_debug_enable(ofdmframesync _q)1116 void ofdmframesync_debug_enable(ofdmframesync _q)
1117 {
1118     // create debugging objects if necessary
1119 #if DEBUG_OFDMFRAMESYNC
1120     if (_q->debug_objects_created)
1121         return;
1122 
1123     _q->debug_x         = windowcf_create(DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1124     _q->debug_rssi      = windowf_create(DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1125     _q->debug_framesyms = windowcf_create(DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1126     _q->G_hat           = (float complex*) malloc((_q->M)*sizeof(float complex));
1127 
1128     _q->px = (float*) malloc((_q->M_pilot)*sizeof(float));
1129     _q->py = (float*) malloc((_q->M_pilot)*sizeof(float));
1130 
1131     _q->debug_pilot_0 = windowf_create(DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1132     _q->debug_pilot_1 = windowf_create(DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1133 
1134     _q->debug_enabled   = 1;
1135     _q->debug_objects_created = 1;
1136 #else
1137     fprintf(stderr,"ofdmframesync_debug_enable(): compile-time debugging disabled\n");
1138 #endif
1139 }
1140 
ofdmframesync_debug_disable(ofdmframesync _q)1141 void ofdmframesync_debug_disable(ofdmframesync _q)
1142 {
1143     // disable debugging
1144 #if DEBUG_OFDMFRAMESYNC
1145     _q->debug_enabled = 0;
1146 #else
1147     fprintf(stderr,"ofdmframesync_debug_disable(): compile-time debugging disabled\n");
1148 #endif
1149 }
1150 
ofdmframesync_debug_print(ofdmframesync _q,const char * _filename)1151 void ofdmframesync_debug_print(ofdmframesync _q,
1152                                const char * _filename)
1153 {
1154 #if DEBUG_OFDMFRAMESYNC
1155     if (!_q->debug_objects_created) {
1156         fprintf(stderr,"error: ofdmframe_debug_print(), debugging objects don't exist; enable debugging first\n");
1157         return;
1158     }
1159 
1160     FILE * fid = fopen(_filename,"w");
1161     if (!fid) {
1162         fprintf(stderr,"error: ofdmframe_debug_print(), could not open '%s' for writing\n", _filename);
1163         return;
1164     }
1165     fprintf(fid,"%% %s : auto-generated file\n", DEBUG_OFDMFRAMESYNC_FILENAME);
1166     fprintf(fid,"close all;\n");
1167     fprintf(fid,"clear all;\n");
1168     fprintf(fid,"n = %u;\n", DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1169     fprintf(fid,"M = %u;\n", _q->M);
1170     fprintf(fid,"M_null  = %u;\n", _q->M_null);
1171     fprintf(fid,"M_pilot = %u;\n", _q->M_pilot);
1172     fprintf(fid,"M_data  = %u;\n", _q->M_data);
1173     unsigned int i;
1174     float complex * rc;
1175     float * r;
1176 
1177     // save subcarrier allocation
1178     fprintf(fid,"p = zeros(1,M);\n");
1179     for (i=0; i<_q->M; i++)
1180         fprintf(fid,"p(%4u) = %d;\n", i+1, _q->p[i]);
1181     fprintf(fid,"i_null  = find(p==%d);\n", OFDMFRAME_SCTYPE_NULL);
1182     fprintf(fid,"i_pilot = find(p==%d);\n", OFDMFRAME_SCTYPE_PILOT);
1183     fprintf(fid,"i_data  = find(p==%d);\n", OFDMFRAME_SCTYPE_DATA);
1184 
1185     // short, long, training sequences
1186     for (i=0; i<_q->M; i++) {
1187         fprintf(fid,"S0(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(_q->S0[i]), cimagf(_q->S0[i]));
1188         fprintf(fid,"S1(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(_q->S1[i]), cimagf(_q->S1[i]));
1189     }
1190 
1191     fprintf(fid,"x = zeros(1,n);\n");
1192     windowcf_read(_q->debug_x, &rc);
1193     for (i=0; i<DEBUG_OFDMFRAMESYNC_BUFFER_LEN; i++)
1194         fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(rc[i]), cimagf(rc[i]));
1195     fprintf(fid,"figure;\n");
1196     fprintf(fid,"plot(0:(n-1),real(x),0:(n-1),imag(x));\n");
1197     fprintf(fid,"xlabel('sample index');\n");
1198     fprintf(fid,"ylabel('received signal, x');\n");
1199 
1200 
1201     fprintf(fid,"s1 = [];\n");
1202     for (i=0; i<_q->M; i++)
1203         fprintf(fid,"s1(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(_q->s1[i]), cimagf(_q->s1[i]));
1204 
1205 
1206     // write agc_rssi
1207     fprintf(fid,"\n\n");
1208     fprintf(fid,"agc_rssi = zeros(1,%u);\n", DEBUG_OFDMFRAMESYNC_BUFFER_LEN);
1209     windowf_read(_q->debug_rssi, &r);
1210     for (i=0; i<DEBUG_OFDMFRAMESYNC_BUFFER_LEN; i++)
1211         fprintf(fid,"agc_rssi(%4u) = %12.4e;\n", i+1, r[i]);
1212     fprintf(fid,"agc_rssi = filter([0.00362168 0.00724336 0.00362168],[1 -1.82269490 0.83718163],agc_rssi);\n");
1213     fprintf(fid,"agc_rssi = 10*log10( agc_rssi );\n");
1214     fprintf(fid,"figure;\n");
1215     fprintf(fid,"plot(agc_rssi)\n");
1216     fprintf(fid,"ylabel('RSSI [dB]');\n");
1217 
1218     // write short, long symbols
1219     fprintf(fid,"\n\n");
1220     fprintf(fid,"S0 = zeros(1,M);\n");
1221     fprintf(fid,"S1 = zeros(1,M);\n");
1222     for (i=0; i<_q->M; i++) {
1223         fprintf(fid,"S0(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(_q->S0[i]), cimagf(_q->S0[i]));
1224         fprintf(fid,"S1(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(_q->S1[i]), cimagf(_q->S1[i]));
1225     }
1226 
1227 
1228     // write gain arrays
1229     fprintf(fid,"\n\n");
1230     fprintf(fid,"G0     = zeros(1,M);\n");
1231     fprintf(fid,"G1     = zeros(1,M);\n");
1232     fprintf(fid,"G_hat  = zeros(1,M);\n");
1233     fprintf(fid,"G      = zeros(1,M);\n");
1234     for (i=0; i<_q->M; i++) {
1235         fprintf(fid,"G0(%3u)    = %12.8f + j*%12.8f;\n", i+1, crealf(_q->G0a[i]),  cimagf(_q->G0a[i]));
1236         fprintf(fid,"G1(%3u)    = %12.8f + j*%12.8f;\n", i+1, crealf(_q->G0b[i]),  cimagf(_q->G0b[i]));
1237         fprintf(fid,"G_hat(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(_q->G_hat[i]),cimagf(_q->G_hat[i]));
1238         fprintf(fid,"G(%3u)     = %12.8f + j*%12.8f;\n", i+1, crealf(_q->G[i]),    cimagf(_q->G[i]));
1239     }
1240     fprintf(fid,"f = [0:(M-1)];\n");
1241     fprintf(fid,"figure;\n");
1242     fprintf(fid,"subplot(2,1,1);\n");
1243     fprintf(fid,"  plot(f, fftshift(abs(G_hat)),'sb',...\n");
1244     fprintf(fid,"       f, fftshift(abs(G)),'-k','LineWidth',2);\n");
1245     fprintf(fid,"  grid on;\n");
1246     fprintf(fid,"  xlabel('subcarrier index');\n");
1247     fprintf(fid,"  ylabel('gain estimate (mag)');\n");
1248     fprintf(fid,"subplot(2,1,2);\n");
1249     fprintf(fid,"  plot(f, fftshift(arg(G_hat).*[abs(G0) > 1e-3]),'sb',...\n");
1250     fprintf(fid,"       f, fftshift(arg(G)),'-k','LineWidth',2);\n");
1251     fprintf(fid,"  grid on;\n");
1252     fprintf(fid,"  xlabel('subcarrier index');\n");
1253     fprintf(fid,"  ylabel('gain estimate (phase)');\n");
1254 
1255     // write pilot response
1256     fprintf(fid,"\n\n");
1257     fprintf(fid,"px = zeros(1,M_pilot);\n");
1258     fprintf(fid,"py = zeros(1,M_pilot);\n");
1259     for (i=0; i<_q->M_pilot; i++) {
1260         fprintf(fid,"px(%3u) = %12.8f;\n", i+1, _q->px[i]);
1261         fprintf(fid,"py(%3u) = %12.8f;\n", i+1, _q->py[i]);
1262     }
1263     fprintf(fid,"p_phase(1) = %12.8f;\n", _q->p_phase[0]);
1264     fprintf(fid,"p_phase(2) = %12.8f;\n", _q->p_phase[1]);
1265 
1266     // save pilot history
1267     fprintf(fid,"p0 = zeros(1,M);\n");
1268     windowf_read(_q->debug_pilot_0, &r);
1269     for (i=0; i<DEBUG_OFDMFRAMESYNC_BUFFER_LEN; i++)
1270         fprintf(fid,"p0(%4u) = %12.4e;\n", i+1, r[i]);
1271 
1272     fprintf(fid,"p1 = zeros(1,M);\n");
1273     windowf_read(_q->debug_pilot_1, &r);
1274     for (i=0; i<DEBUG_OFDMFRAMESYNC_BUFFER_LEN; i++)
1275         fprintf(fid,"p1(%4u) = %12.4e;\n", i+1, r[i]);
1276 
1277     fprintf(fid,"figure;\n");
1278     fprintf(fid,"fp = (-M/2):(M/2);\n");
1279     fprintf(fid,"subplot(3,1,1);\n");
1280     fprintf(fid,"  plot(px, py, 'sb',...\n");
1281     fprintf(fid,"       fp, polyval(p_phase, fp), '-k');\n");
1282     fprintf(fid,"  grid on;\n");
1283     fprintf(fid,"  legend('pilots','polyfit',0);\n");
1284     fprintf(fid,"  xlabel('subcarrier');\n");
1285     fprintf(fid,"  ylabel('phase');\n");
1286     fprintf(fid,"subplot(3,1,2);\n");
1287     fprintf(fid,"  plot(1:length(p0), p0);\n");
1288     fprintf(fid,"  grid on;\n");
1289     fprintf(fid,"  ylabel('p0 (phase offset)');\n");
1290     fprintf(fid,"subplot(3,1,3);\n");
1291     fprintf(fid,"  plot(1:length(p1), p1);\n");
1292     fprintf(fid,"  grid on;\n");
1293     fprintf(fid,"  ylabel('p1 (phase slope)');\n");
1294 
1295     // write frame symbols
1296     fprintf(fid,"framesyms = zeros(1,n);\n");
1297     windowcf_read(_q->debug_framesyms, &rc);
1298     for (i=0; i<DEBUG_OFDMFRAMESYNC_BUFFER_LEN; i++)
1299         fprintf(fid,"framesyms(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(rc[i]), cimagf(rc[i]));
1300     fprintf(fid,"figure;\n");
1301     fprintf(fid,"plot(real(framesyms), imag(framesyms), 'x');\n");
1302     fprintf(fid,"xlabel('I');\n");
1303     fprintf(fid,"ylabel('Q');\n");
1304     fprintf(fid,"axis([-1 1 -1 1]*1.6);\n");
1305     fprintf(fid,"axis square;\n");
1306     fprintf(fid,"grid on;\n");
1307 
1308     fclose(fid);
1309     printf("ofdmframesync/debug: results written to '%s'\n", _filename);
1310 #else
1311     fprintf(stderr,"ofdmframesync_debug_print(): compile-time debugging disabled\n");
1312 #endif
1313 }
1314 
1315