1 //
2 // modem_demodulate_soft_gentab.c
3 //
4 // Generates table of nearest symbols and tests soft demodulation
5 //
6
7 #include <stdlib.h>
8 #include <stdio.h>
9 #include <string.h>
10 #include <math.h>
11 #include <time.h>
12
13 #include "liquid.h"
14
15 #define OUTPUT_FILENAME "modem_demodulate_soft_gentab.m"
16 #define DEBUG 0
17
18 // print a string of bits to the standard output
print_bitstring(unsigned int _x,unsigned int _n)19 void print_bitstring(unsigned int _x,
20 unsigned int _n)
21 {
22 unsigned int i;
23 for (i=0; i<_n; i++)
24 printf("%1u", (_x >> (_n-i-1)) & 1);
25 }
26
27 // soft demodulation
28 // _q : demodulator object
29 // _c : constellation table
30 // _cp : nearest neighbor table [size: _p x 1]
31 // _p : size of table
32 // _r : received sample
33 // _s : hard demodulator output
34 // _soft_bits : soft bit ouput (approximate log-likelihood ratio)
35 void modem_demodulate_soft_tab(modem _q,
36 float complex * _c,
37 unsigned int * _cp,
38 unsigned int _p,
39 float complex _r,
40 unsigned int * _s,
41 float * _soft_bits);
42
main(int argc,char * argv[])43 int main(int argc, char*argv[])
44 {
45 srand(time(NULL));
46
47 // options
48 modulation_scheme ms = LIQUID_MODEM_QAM16; // modulation scheme
49 unsigned int p = 4; // number of 'nearest symbols'
50 float complex e = 0.1f + _Complex_I*0.15f; // error
51
52 // validate input
53 if (ms == LIQUID_MODEM_UNKNOWN || ms >= LIQUID_MODEM_NUM_SCHEMES) {
54 fprintf(stderr,"error: %s, invalid modulation scheme\n", argv[0]);
55 exit(1);
56 }
57
58 unsigned int i;
59 unsigned int j;
60 unsigned int k;
61
62 // derived values
63 unsigned int bps = modulation_types[ms].bps;
64 unsigned int M = 1 << bps; // constellation size
65 // ensure number of nearest symbols is not too large
66 if (p > (M-1)) {
67 fprintf(stderr,"error: requesting too many neighbors\n");
68 exit(1);
69 }
70 float sig = 0.2f; // noise standard deviation
71
72 // generate constellation
73 modem q = modem_create(ms);
74 float complex c[M]; // constellation
75 for (i=0; i<M; i++)
76 modem_modulate(q, i, &c[i]);
77 modem_destroy(q);
78
79 //
80 // find nearest symbols
81 //
82 unsigned int cp[M*p];
83
84 // initialize table (with 'M' as invalid symbol)
85 for (i=0; i<M; i++) {
86 for (k=0; k<p; k++)
87 cp[i*p + k] = M;
88 }
89
90 int symbol_valid;
91 for (i=0; i<M; i++) {
92 #if DEBUG
93 printf("\n******** %3u ********\n", i);
94 #endif
95 for (k=0; k<p; k++) {
96 #if DEBUG
97 printf("iteration %u\n", k);
98 #endif
99 float dmin=1e9f;
100 for (j=0; j<M; j++) {
101 symbol_valid = 1;
102 // ignore this symbol
103 if (i==j) symbol_valid = 0;
104
105 // also ignore symbol if it is already in the index
106 unsigned int l;
107 for (l=0; l<p; l++) {
108 if ( cp[i*p + l] == j )
109 symbol_valid = 0;
110 }
111
112 // compute distance
113 float d = cabsf( c[i] - c[j] );
114 #if DEBUG
115 printf(" testing %3u... ", j);
116 if (!symbol_valid)
117 printf("invalid\n");
118 else
119 printf(" d = %8.6f ", d);
120 #endif
121
122 if ( d < dmin && symbol_valid ) {
123 dmin = d;
124 cp[i*p + k] = j;
125 #if DEBUG
126 printf("*\n");
127 } else {
128 if (symbol_valid)
129 printf("\n");
130 #endif
131 }
132
133 // ...
134 //cp[i*p + k] = rand() % (M+4);
135 }
136 #if DEBUG
137 printf("new symbol : %3u, dmin=%8.6f\n", cp[i*p + k], dmin);
138 #endif
139 }
140 }
141
142 //
143 // print results
144 //
145 for (i=0; i<M; i++) {
146 printf(" %4u [", i);
147 print_bitstring(i,bps);
148 printf("] : ");
149 for (k=0; k<p; k++) {
150 printf(" ");
151 print_bitstring(cp[i*p + k],bps);
152 if (cp[i*p + k] < M)
153 printf("(%6.4f)", cabsf( c[i]-c[cp[i*p+k]] ));
154 else
155 printf(" ");
156 }
157 printf("\n");
158 }
159
160 // print c-type array
161 printf("\n");
162 printf("// %s soft demodulation nearest neighbors (p=%u)\n", modulation_types[ms].name, p);
163 printf("const unsigned char %s_demod_soft_neighbors[%u] = {\n", modulation_types[ms].name, p*M);
164 printf(" ");
165 for (i=0; i<M; i++) {
166 for (k=0; k<p; k++) {
167 printf("%4u", cp[i*p+k]);
168 if (k != p-1) printf(",");
169 }
170 if (i != M-1) {
171 printf(", // ");
172 print_bitstring(i,bps);
173 printf("\n ");
174 } else {
175 printf("}; // ");
176 print_bitstring(i,bps);
177 printf("\n\n");
178 }
179 }
180
181 // select input symbol and compute received symbol
182 unsigned int sym_in = rand() % M;
183 if (ms == LIQUID_MODEM_QAM16)
184 sym_in = 13;
185 float complex r = c[sym_in] + e;
186
187 // run soft demodulation for each bit
188 float soft_bits[bps];
189 for (k=0; k<bps; k++) {
190 #if DEBUG
191 printf("\n");
192 printf("********* bit index %u ************\n", k);
193 #endif
194 // reset soft bit value
195 soft_bits[k] = 0.0f;
196 float bit_0 = 0.0f;
197 float bit_1 = 0.0f;
198
199 // compute LLR for this bit
200 for (i=0; i<M; i++) {
201 // compute distance between received point and symbol
202 float d = crealf( (r-c[i])*conjf(r-c[i]) );
203 float t = expf( -d / (2.0f*sig*sig) );
204
205 // check if this symbol has a '0' or '1' at this bit index
206 unsigned int bit = (i >> (bps-k-1)) & 0x01;
207 //printf("%c", bit ? '1' : '0');
208
209 if (bit) bit_1 += t;
210 else bit_0 += t;
211
212 #if DEBUG
213 printf(" symbol : ");
214 print_bitstring(i,bps);
215 printf(" [%c]", bit ? '1' : '0');
216 printf(" { %7.4f %7.4f}, d=%7.4f, t=%12.8f\n", crealf(c[i]), cimagf(c[i]), d, t);
217 #endif
218 }
219
220 soft_bits[k] = logf(bit_1) - logf(bit_0);
221 #if DEBUG
222 printf(" {0 : %12.8f, 1 : %12.8f}\n", bit_0, bit_1);
223 #endif
224 }
225
226 //
227 // demodulate using internal method
228 //
229 q = modem_create(ms);
230 unsigned int sym_out_tab;
231 float soft_bits_tab[bps];
232 modem_demodulate_soft_tab(q,c,cp,p,r,&sym_out_tab,soft_bits_tab);
233 modem_destroy(q);
234
235 // print results
236 printf("\n");
237 printf(" input symbol : ");
238 print_bitstring(sym_in,bps);
239 printf(" {%12.8f, %12.8f}\n", crealf(c[sym_in]), cimagf(c[sym_in]));
240
241 printf(" soft bits :\n");
242 printf(" /----------- LLR -----------\\ /------- min dist --------\\\n");
243 for (k=0; k<bps; k++) {
244 printf(" %1u : ", (sym_in >> (bps-k-1)) & 0x01);
245 printf("%12.8f > ", soft_bits[k]);
246 int soft_bit = (soft_bits[k]*16 + 127);
247 if (soft_bit > 255) soft_bit = 255;
248 if (soft_bit < 0) soft_bit = 0;
249 printf("%5d > %1u", soft_bit, soft_bit & 0x80 ? 1 : 0);
250
251 //
252 // print tabular method
253 //
254 printf(" ");
255 printf("%12.8f > ", soft_bits_tab[k]);
256 soft_bit = (soft_bits_tab[k]*16 + 127);
257 if (soft_bit > 255) soft_bit = 255;
258 if (soft_bit < 0) soft_bit = 0;
259 printf("%5d > %1u", soft_bit, soft_bit & 0x80 ? 1 : 0);
260
261 printf("\n");
262 }
263
264 //
265 // export results to file
266 //
267 FILE * fid = fopen(OUTPUT_FILENAME,"w");
268 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
269 fprintf(fid,"clear all;\n");
270 fprintf(fid,"close all;\n\n");
271 fprintf(fid,"m = %u;\n", bps);
272 fprintf(fid,"M = %u;\n", 1<<bps);
273 fprintf(fid,"c = zeros(1,M);\n");
274 fprintf(fid,"i_str = cell(1,M);\n");
275
276 for (i=0; i<M; i++) {
277 // write symbol to output file
278 fprintf(fid,"c(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(c[i]), cimagf(c[i]));
279 fprintf(fid,"i_str{%3u} = '", i+1);
280 for (j=0; j<bps; j++)
281 fprintf(fid,"%c", (i >> (bps-j-1)) & 0x01 ? '1' : '0');
282 fprintf(fid,"';\n");
283 }
284 fprintf(fid,"x = %12.8f + j*%12.8f;\n", crealf(c[sym_in]), cimagf(c[sym_in]));
285 fprintf(fid,"r = %12.8f + j*%12.8f;\n", crealf(r), cimagf(r));
286
287 // plot results
288 fprintf(fid,"\n\n");
289 fprintf(fid,"figure;\n");
290 fprintf(fid,"plot(c,'o','MarkerSize',4,r,'rx',[x r]);\n");
291 fprintf(fid,"hold on;\n");
292 // print lines to 'nearest neighbor' points
293 fprintf(fid,"plot(");
294 for (i=0; i<p; i++) {
295 float complex x_hat = c[ cp[sym_out_tab*p+i] ];
296 fprintf(fid,"[%12.8f %12.8f],[%12.8f %12.8f],'Color',[1 1 1]*0.8",
297 crealf(r), crealf(x_hat),
298 cimagf(r), cimagf(x_hat));
299 if (i == p-1) fprintf(fid,");\n");
300 else fprintf(fid,",...\n ");
301 }
302 fprintf(fid,"text(real(c)+0.02, imag(c)+0.02, i_str);\n");
303 fprintf(fid,"hold off;\n");
304 fprintf(fid,"axis([-1 1 -1 1]*1.6);\n");
305 fprintf(fid,"axis square;\n");
306 fprintf(fid,"grid on;\n");
307 fprintf(fid,"xlabel('in phase');\n");
308 fprintf(fid,"ylabel('quadrature phase');\n");
309
310 fclose(fid);
311 printf("results written to %s.\n", OUTPUT_FILENAME);
312
313 printf("done.\n");
314 return 0;
315 }
316
317 // soft demodulation
318 // _q : demodulator object
319 // _c : constellation table
320 // _cp : nearest neighbor table [size: M*_p x 1]
321 // _p : size of table
322 // _r : received sample
323 // _s : hard demodulator output
324 // _soft_bits : soft bit ouput (approximate log-likelihood ratio)
modem_demodulate_soft_tab(modem _q,float complex * _c,unsigned int * _cp,unsigned int _p,float complex _r,unsigned int * _s,float * _soft_bits)325 void modem_demodulate_soft_tab(modem _q,
326 float complex * _c,
327 unsigned int * _cp,
328 unsigned int _p,
329 float complex _r,
330 unsigned int * _s,
331 float * _soft_bits)
332 {
333 #if DEBUG
334 printf("\nmodem_demodulate_soft_tab() invoked\n");
335 #endif
336 // run hard demodulation
337 modem_demodulate(_q, _r, _s);
338 #if DEBUG
339 printf(" hard demod : %3u\n", *_s);
340 #endif
341
342 unsigned int bps = modem_get_bps(_q);
343
344 float sig = 0.2f;
345
346 unsigned int k;
347 unsigned int i;
348 unsigned int bit;
349 float d;
350 for (k=0; k<bps; k++) {
351 // initialize soft bit value
352 _soft_bits[k] = 0.0f;
353
354 // find nearest 0 and nearest 1
355 float dmin_0 = 1.0f;
356 float dmin_1 = 1.0f;
357
358 // check bit of hard demodulation
359 d = crealf( (_r-_c[*_s])*conjf(_r-_c[*_s]) );
360 bit = (*_s >> (bps-k-1)) & 0x01;
361 if (bit) dmin_1 = d;
362 else dmin_0 = d;
363
364 // check symbols in table
365 #if DEBUG
366 printf(" index %2u : ", k);
367 #endif
368 for (i=0; i<_p; i++) {
369 bit = (_cp[(*_s)*_p+i] >> (bps-k-1)) & 0x01;
370
371 #if DEBUG
372 print_bitstring(_cp[(*_s)*_p+i],bps);
373 printf("[%1u]", bit);
374 #endif
375
376 // compute distance
377 float complex x_hat = _c[ _cp[(*_s)*_p + i] ];
378 d = crealf( (_r-x_hat)*conjf(_r-x_hat) );
379 #if DEBUG
380 printf("(%8.6f) ", d);
381 #endif
382 if (bit) {
383 if (d < dmin_1) dmin_1 = d;
384 } else {
385 if (d < dmin_0) dmin_0 = d;
386 }
387 }
388 #if DEBUG
389 printf("\n");
390 printf(" dmin_0 : %12.8f\n", dmin_0);
391 printf(" dmin_1 : %12.8f\n", dmin_1);
392 #endif
393
394 // make assignments
395 #if 0
396 _soft_bits[k] = logf(expf(-dmin_1/(2.0f*sig*sig)))
397 - logf(expf(-dmin_0/(2.0f*sig*sig)));
398 #else
399 _soft_bits[k] = (-dmin_1/(2.0f*sig*sig))
400 - (-dmin_0/(2.0f*sig*sig));
401 #endif
402 }
403 }
404
405