1 /*
2  * Copyright (c) 2007 - 2019 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 // modem_common.c : common utilities specific to precision
25 //
26 
27 #include <stdlib.h>
28 #include <stdio.h>
29 #include <string.h>
30 
31 #include "liquid.internal.h"
32 
33 #define DEBUG_DEMODULATE_SOFT 0
34 
35 // modem structure used for both modulation and demodulation
36 //
37 // The modem structure implements a variety of common modulation schemes,
38 // including (differential) phase-shift keying, and (quadrature) amplitude
39 // modulation.
40 //
41 // While the same modem structure may be used for both modulation and
42 // demodulation for most schemes, it is important to use separate objects
43 // for differential-mode modems (e.g. LIQUID_MODEM_DPSK) as the internal state
44 // will change after each symbol.  It is usually good practice to keep
45 // separate instances of modulators and demodulators.
MODEM(_s)46 struct MODEM(_s)
47 {
48     // common data
49     modulation_scheme scheme;   // modulation scheme
50     unsigned int m;             // bits per symbol (modulation depth)
51     unsigned int M;             // constellation size, M=2^m
52 
53     // Reference vector for demodulating linear arrays
54     //
55     // By storing these values in an array they do not need to be
56     // calculated during run-time.  This speeds up the demodulation by
57     // approximately 8%.
58     T ref[MAX_MOD_BITS_PER_SYMBOL];
59 
60     // modulation
61     TC * symbol_map;     // complete symbol map
62     int modulate_using_map;         // modulate using map (look-up table) flag
63 
64     // demodulation
65     TC r;                // received state vector
66     TC x_hat;            // estimated symbol (demodulator)
67 
68     // common data structure shared between specific modem types
69     union {
70         // PSK modem
71         struct {
72             T d_phi;            // half of phase between symbols
73             T alpha;            // scaling factor for phase symbols
74         } psk;
75 
76         // DPSK modem
77         struct {
78             T d_phi;            // half of phase between symbols
79             T phi;              // angle state for differential PSK
80             T alpha;            // scaling factor for phase symbols
81         } dpsk;
82 
83         // ASK modem
84         struct {
85             T alpha;            // scaling factor to ensure unity energy
86         } ask;
87 
88         // QAM modem
89         struct {
90             unsigned int m_i;   // bits per symbol, in-phase
91             unsigned int m_q;   // bits per symbol, quadrature
92             unsigned int M_i;   // in-phase dimension, M_i=2^{m_i}
93             unsigned int M_q;   // quadrature dimension, M_q=2^{m_q}
94             T alpha;            // scaling factor to ensure unity energy
95         } qam;
96 
97         // APSK modem
98         struct {
99             unsigned int num_levels;    // number of levels
100             unsigned int p[8];          // number of symbols per level
101             T r[8];                     // radii of levels
102             T r_slicer[8];              // slicer radii of levels
103             T phi[8];                   // phase offset of levels
104             unsigned char * map;        // symbol mapping (allocated)
105         } apsk;
106 
107         // 'square' 32-QAM
108         struct {
109             TC * map;           // 8-sample sub-map (first quadrant)
110         } sqam32;
111 
112         // 'square' 128-QAM
113         struct {
114             TC * map;           // 32-sample sub-map (first quadrant)
115         } sqam128;
116     } data;
117 
118     // modulate function pointer
119     void (*modulate_func)(MODEM() _q,
120                           unsigned int _symbol_in,
121                           TC * _y);
122 
123     // demodulate function pointer
124     void (*demodulate_func)(MODEM() _q,
125                             TC _x,
126                             unsigned int * _symbol_out);
127 
128     // soft demodulation
129     //int demodulate_soft;    // soft demodulation flag
130     // neighbors array
131     unsigned char * demod_soft_neighbors;   // array of nearest neighbors
132     unsigned int demod_soft_p;              // number of neighbors in array
133 };
134 
135 // create digital modem of a specific scheme and bits/symbol
MODEM(_create)136 MODEM() MODEM(_create)(modulation_scheme _scheme)
137 {
138     switch (_scheme) {
139 
140     // Phase-shift keying (PSK)
141     case LIQUID_MODEM_PSK2:     return MODEM(_create_psk)(1);
142     case LIQUID_MODEM_PSK4:     return MODEM(_create_psk)(2);
143     case LIQUID_MODEM_PSK8:     return MODEM(_create_psk)(3);
144     case LIQUID_MODEM_PSK16:    return MODEM(_create_psk)(4);
145     case LIQUID_MODEM_PSK32:    return MODEM(_create_psk)(5);
146     case LIQUID_MODEM_PSK64:    return MODEM(_create_psk)(6);
147     case LIQUID_MODEM_PSK128:   return MODEM(_create_psk)(7);
148     case LIQUID_MODEM_PSK256:   return MODEM(_create_psk)(8);
149 
150     // Differential phase-shift keying (DPSK)
151     case LIQUID_MODEM_DPSK2:    return MODEM(_create_dpsk)(1);
152     case LIQUID_MODEM_DPSK4:    return MODEM(_create_dpsk)(2);
153     case LIQUID_MODEM_DPSK8:    return MODEM(_create_dpsk)(3);
154     case LIQUID_MODEM_DPSK16:   return MODEM(_create_dpsk)(4);
155     case LIQUID_MODEM_DPSK32:   return MODEM(_create_dpsk)(5);
156     case LIQUID_MODEM_DPSK64:   return MODEM(_create_dpsk)(6);
157     case LIQUID_MODEM_DPSK128:  return MODEM(_create_dpsk)(7);
158     case LIQUID_MODEM_DPSK256:  return MODEM(_create_dpsk)(8);
159 
160     // amplitude-shift keying (ASK)
161     case LIQUID_MODEM_ASK2:     return MODEM(_create_ask)(1);
162     case LIQUID_MODEM_ASK4:     return MODEM(_create_ask)(2);
163     case LIQUID_MODEM_ASK8:     return MODEM(_create_ask)(3);
164     case LIQUID_MODEM_ASK16:    return MODEM(_create_ask)(4);
165     case LIQUID_MODEM_ASK32:    return MODEM(_create_ask)(5);
166     case LIQUID_MODEM_ASK64:    return MODEM(_create_ask)(6);
167     case LIQUID_MODEM_ASK128:   return MODEM(_create_ask)(7);
168     case LIQUID_MODEM_ASK256:   return MODEM(_create_ask)(8);
169 
170     // rectangular quadrature amplitude-shift keying (QAM)
171     case LIQUID_MODEM_QAM4:     return MODEM(_create_qam)(2);
172     case LIQUID_MODEM_QAM8:     return MODEM(_create_qam)(3);
173     case LIQUID_MODEM_QAM16:    return MODEM(_create_qam)(4);
174     case LIQUID_MODEM_QAM32:    return MODEM(_create_qam)(5);
175     case LIQUID_MODEM_QAM64:    return MODEM(_create_qam)(6);
176     case LIQUID_MODEM_QAM128:   return MODEM(_create_qam)(7);
177     case LIQUID_MODEM_QAM256:   return MODEM(_create_qam)(8);
178 
179     // amplitude phase-shift keying (APSK)
180     case LIQUID_MODEM_APSK4:    return MODEM(_create_apsk)(2);
181     case LIQUID_MODEM_APSK8:    return MODEM(_create_apsk)(3);
182     case LIQUID_MODEM_APSK16:   return MODEM(_create_apsk)(4);
183     case LIQUID_MODEM_APSK32:   return MODEM(_create_apsk)(5);
184     case LIQUID_MODEM_APSK64:   return MODEM(_create_apsk)(6);
185     case LIQUID_MODEM_APSK128:  return MODEM(_create_apsk)(7);
186     case LIQUID_MODEM_APSK256:  return MODEM(_create_apsk)(8);
187 
188     // specific modems
189     case LIQUID_MODEM_BPSK:      return MODEM(_create_bpsk)();
190     case LIQUID_MODEM_QPSK:      return MODEM(_create_qpsk)();
191     case LIQUID_MODEM_OOK:       return MODEM(_create_ook)();
192     case LIQUID_MODEM_SQAM32:    return MODEM(_create_sqam32)();
193     case LIQUID_MODEM_SQAM128:   return MODEM(_create_sqam128)();
194     case LIQUID_MODEM_V29:       return MODEM(_create_V29)();
195     case LIQUID_MODEM_ARB16OPT:  return MODEM(_create_arb16opt)();
196     case LIQUID_MODEM_ARB32OPT:  return MODEM(_create_arb32opt)();
197     case LIQUID_MODEM_ARB64OPT:  return MODEM(_create_arb64opt)();
198     case LIQUID_MODEM_ARB128OPT: return MODEM(_create_arb128opt)();
199     case LIQUID_MODEM_ARB256OPT: return MODEM(_create_arb256opt)();
200     case LIQUID_MODEM_ARB64VT:   return MODEM(_create_arb64vt)();
201 
202     // arbitrary modem
203     case LIQUID_MODEM_ARB:
204         fprintf(stderr,"error: modem_create(), cannot create arbitrary modem (LIQUID_MODEM_ARB) without specifying constellation\n");
205         exit(1);
206 
207     // unknown modulation scheme
208     default:
209         fprintf(stderr,"error: modem_create(), unknown/unsupported modulation scheme : %u\n", _scheme);
210         exit(1);
211     }
212 
213     // should never get to this point, but adding return statment
214     // to keep compiler happy
215     return NULL;
216 }
217 
218 // recreate modulation scheme, re-allocating memory as necessary
MODEM(_recreate)219 MODEM() MODEM(_recreate)(MODEM() _q,
220                          modulation_scheme _scheme)
221 {
222     // TODO : regenerate modem only when truly necessary
223     if (_q->scheme != _scheme) {
224         // destroy and re-create modem
225         MODEM(_destroy)(_q);
226         _q = MODEM(_create)(_scheme);
227     }
228 
229     // return object
230     return _q;
231 }
232 
233 // destroy a modem object
MODEM(_destroy)234 void MODEM(_destroy)(MODEM() _q)
235 {
236     // free symbol map
237     if (_q->symbol_map != NULL)
238         free(_q->symbol_map);
239 
240     // free soft-demodulation neighbors table
241     if (_q->demod_soft_neighbors != NULL)
242         free(_q->demod_soft_neighbors);
243 
244     // free memory in specific data types
245     if (_q->scheme == LIQUID_MODEM_SQAM32) {
246         free(_q->data.sqam32.map);
247     } else if (_q->scheme == LIQUID_MODEM_SQAM128) {
248         free(_q->data.sqam128.map);
249     } else if (liquid_modem_is_apsk(_q->scheme)) {
250         free(_q->data.apsk.map);
251     }
252 
253     // free main object memory
254     free(_q);
255 }
256 
257 // print a modem object
MODEM(_print)258 void MODEM(_print)(MODEM() _q)
259 {
260     printf("linear modem:\n");
261     printf("    scheme:         %s\n", modulation_types[_q->scheme].name);
262     printf("    bits/symbol:    %u\n", _q->m);
263 }
264 
265 // reset a modem object (only an issue with dpsk)
MODEM(_reset)266 void MODEM(_reset)(MODEM() _q)
267 {
268     _q->r = 1.0f;         // received sample
269     _q->x_hat = _q->r;  // estimated symbol
270 
271     if ( liquid_modem_is_dpsk(_q->scheme) )
272         _q->data.dpsk.phi = 0.0f;  // reset differential PSK phase state
273 }
274 
275 // initialize a generic modem object
MODEM(_init)276 void MODEM(_init)(MODEM() _q,
277                   unsigned int _bits_per_symbol)
278 {
279     if (_bits_per_symbol < 1 ) {
280         fprintf(stderr,"error: modem_init(), modem must have at least 1 bit/symbol\n");
281         exit(1);
282     } else if (_bits_per_symbol > MAX_MOD_BITS_PER_SYMBOL) {
283         fprintf(stderr,"error: modem_init(), maximum number of bits per symbol exceeded\n");
284         exit(1);
285     }
286 
287     // initialize common elements
288     _q->symbol_map = NULL;    // symbol map (LIQUID_MODEM_ARB only)
289     _q->modulate_using_map=0; // modulate using map flag
290 
291     // common data
292     _q->m = _bits_per_symbol; // bits/symbol
293     _q->M = 1 << (_q->m);   // constellation size (2^m)
294 
295     // set function pointers initially to NULL
296     _q->modulate_func = NULL;
297     _q->demodulate_func = NULL;
298 
299     // soft demodulation
300     _q->demod_soft_neighbors = NULL;
301     _q->demod_soft_p = 0;
302 }
303 
304 // initialize symbol map for fast modulation
MODEM(_init_map)305 void MODEM(_init_map)(MODEM() _q)
306 {
307     // validate input
308     if (_q->symbol_map == NULL) {
309         fprintf(stderr,"error: modem_init_map(), symbol map array has not been allocated\n");
310         exit(1);
311     } else if (_q->M == 0 || _q->M > (1<<MAX_MOD_BITS_PER_SYMBOL)) {
312         fprintf(stderr,"error: modem_init_map(), constellation size is out of range\n");
313         exit(1);
314     } else if (_q->modulate_func == NULL) {
315         fprintf(stderr,"error: modem_init_map(), modulation function has not been initialized\n");
316         exit(1);
317     }
318 
319     unsigned int i;
320     for (i=0; i<_q->M; i++)
321         _q->modulate_func(_q, i, &_q->symbol_map[i]);
322 }
323 
324 // Generate random symbol
MODEM(_gen_rand_sym)325 unsigned int MODEM(_gen_rand_sym)(MODEM() _q)
326 {
327     return rand() % (_q->M);
328 }
329 
330 // Get modem depth (bits/symbol)
MODEM(_get_bps)331 unsigned int MODEM(_get_bps)(MODEM() _q)
332 {
333     return _q->m;
334 }
335 
336 // get modulation scheme
MODEM(_get_scheme)337 modulation_scheme MODEM(_get_scheme)(MODEM() _q)
338 {
339     return _q->scheme;
340 }
341 
342 // generic modulatio function
343 //  _q          :   modem object
344 //  _symbol_in  :   input symbol
345 //  _y          :   output sample
MODEM(_modulate)346 void MODEM(_modulate)(MODEM() _q,
347                       unsigned int _symbol_in,
348                       TC * _y)
349 {
350     // validate input
351     if (_symbol_in >= _q->M) {
352         fprintf(stderr,"error: modem_modulate(), input symbol exceeds constellation size\n");
353         exit(1);
354     }
355 
356     if (_q->modulate_using_map) {
357         // modulate simply using map (look-up table)
358         MODEM(_modulate_map)(_q, _symbol_in, _y);
359     } else {
360         // invoke method specific to scheme (calculate symbol on the fly)
361         _q->modulate_func(_q, _symbol_in, _y);
362     }
363 }
364 
365 // modulate using symbol map (look-up table)
MODEM(_modulate_map)366 void MODEM(_modulate_map)(MODEM() _q,
367                           unsigned int _symbol_in,
368                           TC * _y)
369 {
370     if (_symbol_in >= _q->M) {
371         fprintf(stderr,"error: modem_modulate_table(), input symbol exceeds maximum\n");
372         exit(1);
373     } else if (_q->symbol_map == NULL) {
374         fprintf(stderr,"error: modem_modulate_table(), symbol table not initialized\n");
375         exit(1);
376     }
377 
378     // map sample directly to output
379     *_y = _q->symbol_map[_symbol_in];
380 }
381 
382 // generic demodulation
MODEM(_demodulate)383 void MODEM(_demodulate)(MODEM() _q,
384                         TC x,
385                         unsigned int *symbol_out)
386 {
387     // invoke method specific to scheme (calculate symbol on the fly)
388     _q->demodulate_func(_q, x, symbol_out);
389 }
390 
391 // generic soft demodulation
MODEM(_demodulate_soft)392 void MODEM(_demodulate_soft)(MODEM() _q,
393                              TC _x,
394                              unsigned int  * _s,
395                              unsigned char * _soft_bits)
396 {
397     // switch scheme
398     switch (_q->scheme) {
399     case LIQUID_MODEM_ARB:  MODEM(_demodulate_soft_arb)( _q,_x,_s,_soft_bits); return;
400     case LIQUID_MODEM_BPSK: MODEM(_demodulate_soft_bpsk)(_q,_x,_s,_soft_bits); return;
401     case LIQUID_MODEM_QPSK: MODEM(_demodulate_soft_qpsk)(_q,_x,_s,_soft_bits); return;
402     default:;
403     }
404 
405     // check if...
406     if (_q->demod_soft_neighbors != NULL && _q->demod_soft_p != 0) {
407         // demodulate using approximate log-likelihood method with
408         // look-up table for nearest neighbors
409         MODEM(_demodulate_soft_table)(_q, _x, _s, _soft_bits);
410 
411         return;
412     }
413 
414     // for now demodulate normally and simply copy the
415     // hard-demodulated bits
416     unsigned int symbol_out;
417     _q->demodulate_func(_q, _x, &symbol_out);
418     *_s = symbol_out;
419 
420     // unpack soft bits
421     liquid_unpack_soft_bits(symbol_out, _q->m, _soft_bits);
422 }
423 
424 #if DEBUG_DEMODULATE_SOFT
425 // print a string of bits to the standard output
print_bitstring_demod_soft(unsigned int _x,unsigned int _n)426 void print_bitstring_demod_soft(unsigned int _x,
427                                 unsigned int _n)
428 {
429     unsigned int i;
430     for (i=0; i<_n; i++)
431         printf("%1u", (_x >> (_n-i-1)) & 1);
432 }
433 #endif
434 
435 // generic soft demodulation using look-up table...
436 //  _q          :   demodulator object
437 //  _r          :   received sample
438 //  _s          :   hard demodulator output
439 //  _soft_bits  :   soft bit ouput (approximate log-likelihood ratio)
MODEM(_demodulate_soft_table)440 void MODEM(_demodulate_soft_table)(MODEM() _q,
441                                    TC _r,
442                                    unsigned int * _s,
443                                    unsigned char * _soft_bits)
444 {
445     // run hard demodulation; this will store re-modulated sample
446     // as internal variable x_hat
447     unsigned int s;
448     MODEM(_demodulate)(_q, _r, &s);
449 
450     unsigned int bps = MODEM(_get_bps)(_q);
451 
452     // gamma = 1/(2*sigma^2), approximate for constellation size
453     T gamma = 1.2f*_q->M;
454 
455     // set and initialize minimum bit values
456     unsigned int i;
457     unsigned int k;
458     T dmin_0[bps];
459     T dmin_1[bps];
460     for (k=0; k<bps; k++) {
461         dmin_0[k] = 8.0f;
462         dmin_1[k] = 8.0f;
463     }
464 
465     unsigned int bit;
466     T d;
467     TC x_hat;    // re-modulated symbol
468     unsigned char * softab = _q->demod_soft_neighbors;
469     unsigned int p = _q->demod_soft_p;
470 
471     // check hard demodulation
472     d = crealf( (_r-_q->x_hat)*conjf(_r-_q->x_hat) );
473     for (k=0; k<bps; k++) {
474         bit = (s >> (bps-k-1)) & 0x01;
475         if (bit) dmin_1[k] = d;
476         else     dmin_0[k] = d;
477     }
478 
479     // parse all 'nearest neighbors' and find minimum distance for each bit
480     for (i=0; i<p; i++) {
481         // remodulate symbol
482         if (_q->modulate_using_map)
483             x_hat = _q->symbol_map[ softab[s*p + i] ];
484         else
485             MODEM(_modulate)(_q, softab[s*p+i], &x_hat);
486 
487         // compute magnitude squared of Euclidean distance
488         //d = crealf( (_r-x_hat)*conjf(_r-x_hat) );
489         // (same as above, but faster)
490         TC e = _r - x_hat;
491         d = crealf(e)*crealf(e) + cimagf(e)*cimagf(e);
492 
493         // look at each bit in 'nearest neighbor' and update minimum
494         for (k=0; k<bps; k++) {
495             // strip bit
496             unsigned int bit = (softab[s*p+i] >> (bps-k-1)) & 0x01;
497             if ( bit ) {
498                 if (d < dmin_1[k]) dmin_1[k] = d;
499             } else {
500                 if (d < dmin_0[k]) dmin_0[k] = d;
501             }
502         }
503     }
504 
505     // make soft bit assignments
506     for (k=0; k<bps; k++) {
507         int soft_bit = ((dmin_0[k] - dmin_1[k])*gamma)*16 + 127;
508         if (soft_bit > 255) soft_bit = 255;
509         if (soft_bit <   0) soft_bit = 0;
510         _soft_bits[k] = (unsigned char)soft_bit;
511     }
512 
513     // set hard output symbol
514     *_s = s;
515 }
516 
517 
518 
519 // get demodulator's estimated transmit sample
MODEM(_get_demodulator_sample)520 void MODEM(_get_demodulator_sample)(MODEM() _q,
521                                     TC * _x_hat)
522 {
523     *_x_hat = _q->x_hat;
524 }
525 
526 // get demodulator phase error
MODEM(_get_demodulator_phase_error)527 T MODEM(_get_demodulator_phase_error)(MODEM() _q)
528 {
529     return cimagf(_q->r*conjf(_q->x_hat));
530 }
531 
532 // get error vector magnitude
MODEM(_get_demodulator_evm)533 T MODEM(_get_demodulator_evm)(MODEM() _q)
534 {
535     return cabsf(_q->x_hat - _q->r);
536 }
537 
538 // Demodulate a linear symbol constellation using dynamic threshold calculation
539 //  _v      :   input value
540 //  _m      :   bits per symbol
541 //  _alpha  :   scaling factor
542 //  _s      :   demodulated symbol
543 //  _res    :   residual
MODEM(_demodulate_linear_array)544 void MODEM(_demodulate_linear_array)(T              _v,
545                                      unsigned int   _m,
546                                      T              _alpha,
547                                      unsigned int * _s,
548                                      T *            _res)
549 {
550     unsigned int s=0;
551     unsigned int i, k = _m;
552     T ref=0.0f;
553     for (i=0; i<_m; i++) {
554         s <<= 1;
555         s |= (_v > 0) ? 1 : 0;
556         ref = _alpha * (1<<(k-1));
557         _v += (_v < 0) ? ref : -ref;
558         k--;
559     }
560     *_s = s;
561     *_res = _v;
562 }
563 
564 // Demodulate a linear symbol constellation using refereneced lookup table
565 //  _v      :   input value
566 //  _m      :   bits per symbol
567 //  _ref    :   array of thresholds
568 //  _s      :   demodulated symbol
569 //  _res    :   residual
MODEM(_demodulate_linear_array_ref)570 void MODEM(_demodulate_linear_array_ref)(T              _v,
571                                          unsigned int   _m,
572                                          T *            _ref,
573                                          unsigned int * _s,
574                                          T *            _res)
575 {
576     // initialize loop counter
577     register unsigned int i;
578 
579     // initialize demodulated symbol
580     register unsigned int s=0;
581 
582     for (i=0; i<_m; i++) {
583         // prepare symbol for next demodulated bit
584         s <<= 1;
585 
586         // compare received value to zero
587         if ( _v > 0 ) {
588             // shift '1' into symbol, subtract reference
589             s |= 1;
590             _v -= _ref[_m-i-1];
591         } else {
592             // shift '0' into symbol, add reference
593             s |= 0;
594             _v += _ref[_m-i-1];
595         }
596     }
597     // return demodulated symbol
598     *_s = s;
599 
600     // return residual
601     *_res = _v;
602 }
603 
604 
605 // generate soft demodulation look-up table
MODEM(_demodsoft_gentab)606 void MODEM(_demodsoft_gentab)(MODEM()      _q,
607                               unsigned int _p)
608 {
609     // validate input: ensure number of nearest symbols is not too large
610     if (_p > (_q->M-1)) {
611         fprintf(stderr,"error: modem_demodsoft_gentab(), requesting too many neighbors\n");
612         exit(1);
613     }
614 
615     // allocate internal memory
616     _q->demod_soft_p = _p;
617     _q->demod_soft_neighbors = (unsigned char*)malloc(_q->M*_p*sizeof(unsigned char));
618 
619     unsigned int i;
620     unsigned int j;
621     unsigned int k;
622 
623     // generate constellation
624     // TODO : enforce full constellation for modulation
625     unsigned int M = _q->M;  // constellation size
626     TC c[M];         // constellation
627     for (i=0; i<M; i++)
628         modem_modulate(_q, i, &c[i]);
629 
630     //
631     // find nearest symbols (see algorithm in sandbox/modem_demodulate_soft_gentab.c)
632     //
633 
634     // initialize table (with 'M' as invalid symbol)
635     for (i=0; i<M; i++) {
636         for (k=0; k<_p; k++)
637             _q->demod_soft_neighbors[i*_p + k] = M;
638     }
639 
640     int symbol_valid;
641     for (i=0; i<M; i++) {
642         for (k=0; k<_p; k++) {
643             T dmin=1e9f;
644             for (j=0; j<M; j++) {
645                 symbol_valid = 1;
646                 // ignore this symbol
647                 if (i==j) symbol_valid = 0;
648 
649                 // also ignore symbol if it is already in the index
650                 unsigned int l;
651                 for (l=0; l<_p; l++) {
652                     if ( _q->demod_soft_neighbors[i*_p + l] == j )
653                         symbol_valid = 0;
654                 }
655 
656                 // compute distance
657                 T d = cabsf( c[i] - c[j] );
658 
659                 if ( d < dmin && symbol_valid ) {
660                     dmin = d;
661                     _q->demod_soft_neighbors[i*_p + k] = j;
662                 }
663 
664             }
665         }
666     }
667 
668 #if DEBUG_DEMODULATE_SOFT
669     //
670     // print results
671     //
672     unsigned int bps = _q->m;
673     for (i=0; i<M; i++) {
674         printf(" %4u [", i);
675         print_bitstring_demod_soft(i,bps);
676         printf("] : ");
677         for (k=0; k<_p; k++) {
678             printf(" ");
679             print_bitstring_demod_soft(_q->demod_soft_neighbors[i*_p + k],bps);
680             if (_q->demod_soft_neighbors[i*_p + k] < M)
681                 printf("(%6.4f)", cabsf( c[i]-c[_q->demod_soft_neighbors[i*_p+k]] ));
682             else
683                 printf("        ");
684         }
685         printf("\n");
686     }
687 
688     // print c-type array
689     printf("\n");
690     printf("// %s%u soft demodulation nearest neighbors (p=%u)\n", modulation_types[_q->scheme].name, M, _p);
691     printf("const unsigned char %s%u_demod_soft_neighbors[%u] = {\n", modulation_types[_q->scheme].name, M, _p*M);
692     printf("    ");
693     for (i=0; i<M; i++) {
694         for (k=0; k<_p; k++) {
695             printf("%4u", _q->demod_soft_neighbors[i*_p+k]);
696             if (k != _p-1) printf(",");
697         }
698         if (i != M-1) {
699             printf(",   // ");
700             print_bitstring_demod_soft(i,bps);
701             printf("\n    ");
702         } else {
703             printf("};  // ");
704             print_bitstring_demod_soft(i,bps);
705             printf("\n\n");
706         }
707     }
708 #endif
709 }
710 
711 
712