SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
sdr.h
Go to the documentation of this file.
1 // This file is part of LeanSDR Copyright (C) 2016-2019 <pabr@pabr.org>.
2 // See the toplevel README for more information.
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16 
17 #ifndef LEANSDR_SDR_H
18 #define LEANSDR_SDR_H
19 
20 #include "leansdr/dsp.h"
21 #include "leansdr/math.h"
22 
23 namespace leansdr
24 {
25 
26 // Abbreviations for floating-point types
27 
28 typedef float f32;
29 
30 typedef complex<u8> cu8;
31 typedef complex<s8> cs8;
35 
37 // SDR blocks
39 
40 // AUTO-NOTCH FILTER
41 
42 // Periodically detects the [nslots] strongest peaks with a FFT,
43 // removes them with a first-order filter.
44 
45 template <typename T>
47 {
49  float k;
50 
52  pipebuf<complex<T>> &_in,
53  pipebuf<complex<T>> &_out,
54  int _nslots,
55  T _agc_rms_setpoint) : runnable(sch, "auto_notch"),
56  decimation(1024 * 4096),
57  k(0.002), // k(0.01)
58  fft(4096),
59  in(_in),
60  out(_out, fft.n),
61  nslots(_nslots),
62  phase(0),
63  gain(1),
64  agc_rms_setpoint(_agc_rms_setpoint)
65  {
66  __slots = new slot[nslots];
67 
68  for (int s = 0; s < nslots; ++s)
69  {
70  __slots[s].i = -1;
71  __slots[s].expj = new complex<float>[fft.n];
72  }
73  }
74 
76  {
77  delete[] __slots;
78  }
79 
80  void run()
81  {
82  while (in.readable() >= fft.n && out.writable() >= fft.n)
83  {
84  phase += fft.n;
85 
86  if (phase >= decimation)
87  {
88  phase -= decimation;
89  detect();
90  }
91 
92  process();
93  in.read(fft.n);
94  out.written(fft.n);
95  }
96  }
97 
98  void detect()
99  {
100  complex<T> *pin = in.rd();
101  complex<float> *data = new complex<float>[fft.n];
102  float m0 = 0, m2 = 0;
103 
104  for (int i = 0; i < fft.n; ++i)
105  {
106  data[i].re = pin[i].re;
107  data[i].im = pin[i].im;
108  m2 += (float)pin[i].re * pin[i].re + (float)pin[i].im * pin[i].im;
109  if (gen_abs(pin[i].re) > m0)
110  m0 = gen_abs(pin[i].re);
111  if (gen_abs(pin[i].im) > m0)
112  m0 = gen_abs(pin[i].im);
113  }
114 
115  if (agc_rms_setpoint && m2)
116  {
117  float rms = gen_sqrt(m2 / fft.n);
118  if (sch->debug)
119  fprintf(stderr, "(pow %f max %f)", rms, m0);
120  float new_gain = agc_rms_setpoint / rms;
121  gain = gain * 0.9 + new_gain * 0.1;
122  }
123 
124  fft.inplace(data, true);
125  float *amp = new float[fft.n];
126 
127  for (int i = 0; i < fft.n; ++i)
128  amp[i] = hypotf(data[i].re, data[i].im);
129 
130  for (slot *s = __slots; s < __slots + nslots; ++s)
131  {
132  int iamax = 0;
133 
134  for (int i = 0; i < fft.n; ++i)
135  if (amp[i] > amp[iamax])
136  iamax = i;
137 
138  if (iamax != s->i)
139  {
140  if (sch->debug)
141  fprintf(stderr, "%s: slot %d new peak %d -> %d\n", name, (int)(s - __slots), s->i, iamax);
142 
143  s->i = iamax;
144  s->estim.re = 0;
145  s->estim.im = 0;
146  s->estt = 0;
147 
148  for (int i = 0; i < fft.n; ++i)
149  {
150  float a = 2 * M_PI * s->i * i / fft.n;
151  s->expj[i].re = cosf(a);
152  s->expj[i].im = sinf(a);
153  }
154  }
155 
156  amp[iamax] = 0;
157 
158  if (iamax - 1 >= 0)
159  amp[iamax - 1] = 0;
160 
161  if (iamax + 1 < fft.n)
162  amp[iamax + 1] = 0;
163  }
164 
165  delete[] amp;
166  delete[] data;
167  }
168 
169  void process()
170  {
171  complex<T> *pin = in.rd(), *pend = pin + fft.n, *pout = out.wr();
172 
173  for (slot *s = __slots; s < __slots + nslots; ++s)
174  s->ej = s->expj;
175 
176  for (; pin < pend; ++pin, ++pout)
177  {
178  complex<float> out = *pin;
179  // TODO Optimize for nslots==1 ?
180 
181  for (slot *s = __slots; s < __slots + nslots; ++s->ej, ++s)
182  {
183  complex<float> bb(pin->re * s->ej->re + pin->im * s->ej->im,
184  -pin->re * s->ej->im + pin->im * s->ej->re);
185  s->estim.re = bb.re * k + s->estim.re * (1 - k);
186  s->estim.im = bb.im * k + s->estim.im * (1 - k);
187  complex<float> sub(
188  s->estim.re * s->ej->re - s->estim.im * s->ej->im,
189  s->estim.re * s->ej->im + s->estim.im * s->ej->re);
190  out.re -= sub.re;
191  out.im -= sub.im;
192  }
193 
194  pout->re = gain * out.re;
195  pout->im = gain * out.im;
196  }
197  }
198 
199  private:
203  int nslots;
204 
205  struct slot
206  {
207  int i;
210  int estt;
211  };
212 
213  slot *__slots;
214  int phase;
215  float gain;
217 };
218 
219 // SIGNAL STRENGTH ESTIMATOR
220 
221 // Outputs RMS values.
222 
223 template <typename T>
225 {
226  unsigned long window_size; // Samples per estimation
227  unsigned long decimation; // Output rate
228 
229  ss_estimator(scheduler *sch, pipebuf<complex<T>> &_in, pipebuf<T> &_out) : runnable(sch, "SS estimator"),
230  window_size(1024),
231  decimation(1024),
232  in(_in),
233  out(_out),
234  phase(0)
235  {
236  }
237 
238  void run()
239  {
240  while (in.readable() >= window_size && out.writable() >= 1)
241  {
242  phase += window_size;
243 
244  if (phase >= decimation)
245  {
246  phase -= decimation;
247  complex<T> *p = in.rd(), *pend = p + window_size;
248  float s = 0;
249 
250  for (; p < pend; ++p)
251  s += (float)p->re * p->re + (float)p->im * p->im;
252 
253  out.write(sqrtf(s / window_size));
254  }
255  in.read(window_size);
256  }
257  }
258 
259  private:
262  unsigned long phase;
263 };
264 
265 template <typename T>
267 {
268  unsigned long window_size; // Samples per estimation
269  unsigned long decimation; // Output rate
270 
272  pipebuf<complex<T>> &_in,
273  pipebuf<T> &_out_ss,
274  pipebuf<T> &_out_ampmin,
275  pipebuf<T> &_out_ampmax) : runnable(sch, "SS estimator"),
276  window_size(1024),
277  decimation(1024),
278  in(_in),
279  out_ss(_out_ss),
280  out_ampmin(_out_ampmin),
281  out_ampmax(_out_ampmax),
282  phase(0)
283  {
284  }
285 
286  void run()
287  {
288  while (in.readable() >= window_size && out_ss.writable() >= 1 && out_ampmin.writable() >= 1 && out_ampmax.writable() >= 1)
289  {
290  phase += window_size;
291 
292  if (phase >= decimation)
293  {
294  phase -= decimation;
295  complex<T> *p = in.rd(), *pend = p + window_size;
296  float s2 = 0;
297  float amin = 1e38, amax = 0;
298 
299  for (; p < pend; ++p)
300  {
301  float mag2 = (float)p->re * p->re + (float)p->im * p->im;
302  s2 += mag2;
303  float mag = sqrtf(mag2);
304  if (mag < amin)
305  amin = mag;
306  if (mag > amax)
307  amax = mag;
308  }
309 
310  out_ss.write(sqrtf(s2 / window_size));
311  out_ampmin.write(amin);
312  out_ampmax.write(amax);
313  }
314 
315  in.read(window_size);
316  }
317  }
318 
319  private:
321  pipewriter<T> out_ss, out_ampmin, out_ampmax;
322  unsigned long phase;
323 };
324 
325 // AGC
326 
327 template <typename T>
329 {
330  float out_rms; // Desired RMS output power
331  float bw; // Bandwidth
332  float estimated; // Input power
333  static const int chunk_size = 128;
334 
336  pipebuf<complex<T>> &_in,
337  pipebuf<complex<T>> &_out) : runnable(sch, "AGC"),
338  out_rms(1),
339  bw(0.001),
340  estimated(0),
341  in(_in),
342  out(_out, chunk_size)
343  {
344  }
345 
346  private:
349 
350  void run()
351  {
352  while (in.readable() >= chunk_size && out.writable() >= chunk_size)
353  {
354  complex<T> *pin = in.rd(), *pend = pin + chunk_size;
355  float amp2 = 0;
356 
357  for (; pin < pend; ++pin)
358  amp2 += pin->re * pin->re + pin->im * pin->im;
359 
360  amp2 /= chunk_size;
361 
362  if (!estimated)
363  estimated = amp2;
364 
365  estimated = estimated * (1 - bw) + amp2 * bw;
366  float gain = estimated ? out_rms / sqrtf(estimated) : 0;
367  pin = in.rd();
368  complex<T> *pout = out.wr();
369  float bwcomp = 1 - bw;
370 
371  for (; pin < pend; ++pin, ++pout)
372  {
373  pout->re = pin->re * gain;
374  pout->im = pin->im * gain;
375  }
376 
377  in.read(chunk_size);
378  out.written(chunk_size);
379  }
380  }
381 };
382 // simple_agc
383 
384 typedef uint16_t u_angle; // [0,2PI[ in 65536 steps
385 typedef int16_t s_angle; // [-PI,PI[ in 65536 steps
386 
387 // GENERIC CONSTELLATION DECODING BY LOOK-UP TABLE.
388 
389 // Metrics and phase errors are pre-computed on a RxR grid.
390 // R must be a power of 2.
391 // Up to 256 symbols.
392 
394 { // TBD obsolete
395  int16_t cost; // For Viterbi with TBM=int16_t
397 };
398 
399 // Target RMS amplitude for AGC
400 //const float cstln_amp = 73; // Best for 32APSK 9/10
401 //const float cstln_amp = 90; // Best for QPSK
402 //const float cstln_amp = 64; // Best for BPSK
403 //const float cstln_amp = 75; // Best for BPSK at 45°
404 const float cstln_amp = 75; // Trade-off
405 
406 // A struct that temporarily holds all the info we precompute for the LUT.
407 struct full_ss
408 {
409  uint8_t nearest; // Index of nearest in constellation
410  uint16_t dists2[256]; // Squared distances
411  float p[8]; // 0..1 probability of bits being 1
412 };
413 
414 // Options for soft-symbols.
415 // These functions are overloaded to keep cstln_lut<SOFTSYMB> generic:
416 // to_softsymb(const full_ss *fss, SOFTSYMB *ss)
417 // softsymb_harden(SOFTSYMB *ss) {
418 // softsymb_to_dump(const SOFTSYMB &ss, int bit) To grey 0..255
419 // For LUT initialization only. Performance is not critical.
420 
421 // Hard decision soft-symbols.
422 // Value is the symbol index, 0..255.
423 typedef uint8_t hard_ss;
424 void to_softsymb(const full_ss *fss, hard_ss *ss);
425 void softsymb_harden(hard_ss *ss);
426 uint8_t softsymb_to_dump(const hard_ss &ss, int bit);
427 
428 // Euclidian QPSK soft-symbols.
429 // Additive metric suitable for Viterbi.
430 // Backward-compatible with simplified Viterbi (TBD remove)
431 struct eucl_ss
432 {
433  static const int MAX_SYMBOLS = 4;
434  uint16_t dists2[MAX_SYMBOLS];
435  uint16_t discr2; // 2nd_nearest - nearest
437 };
438 
439 void to_softsymb(const full_ss *fss, eucl_ss *ss);
440 void softsymb_harden(eucl_ss *ss);
441 uint8_t softsymb_to_dump(const eucl_ss &ss, int bit);
442 
443 // Log-Likelihood Ratios soft-symbols
444 typedef int8_t llr_t; // log(p(0)/p(1)), clip -127=1 +127=0
445 
446 inline bool llr_harden(llr_t v)
447 {
448  return v & 128;
449 }
450 
451 struct llr_ss
452 {
453  llr_t bits[8]; // Up to 8 bit considered independent
454 };
455 
456 void to_softsymb(const full_ss *fss, llr_ss *ss);
457 void softsymb_harden(llr_ss *ss);
458 uint8_t softsymb_to_dump(const llr_ss &ss, int bit);
459 
461 {
462  enum predef
463  {
464  BPSK, // DVB-S2 (and DVB-S variant)
465  QPSK, // DVB-S
468  APSK32, // DVB-S2
469  APSK64E, // DVB-S2X
472  QAM256, // For experimentation only
473  COUNT
474  };
475 
476  static const char *names[];
477  float amp_max; // Max amplitude. 1 for PSK, 0 if not applicable.
479  int nsymbols;
481 };
482 // cstln_base
483 
484 template <typename SOFTSYMB, int R>
486 {
488  float mer = 10,
489  float gamma1 = 0,
490  float gamma2 = 0,
491  float gamma3 = 0)
492  {
493  switch (type)
494  {
495  case BPSK:
496  amp_max = 1;
497  nrotations = 2;
498  nsymbols = 2;
499  symbols = new complex<signed char>[nsymbols];
500 #if 0 // BPSK at 0°
501  symbols[0] = polar(1, 2, 0);
502  symbols[1] = polar(1, 2, 1);
503  printf("cstln_lut::cstln_lut: BPSK at 0 degrees\n");
504 #else // BPSK at 45°
505  symbols[0] = polar(1, 8, 1);
506  symbols[1] = polar(1, 8, 5);
507  printf("cstln_lut::cstln_lut: BPSK at 45 degrees\n");
508 #endif
509  make_lut_from_symbols(mer);
510  break;
511  case QPSK:
512  amp_max = 1;
513  // EN 300 421, section 4.5 Baseband shaping and modulation
514  // EN 302 307, section 5.4.1
515  nrotations = 4;
516  nsymbols = 4;
517  symbols = new complex<signed char>[nsymbols];
518  symbols[0] = polar(1, 4, 0.5);
519  symbols[1] = polar(1, 4, 3.5);
520  symbols[2] = polar(1, 4, 1.5);
521  symbols[3] = polar(1, 4, 2.5);
522  make_lut_from_symbols(mer);
523  printf("cstln_lut::cstln_lut: QPSK\n");
524  break;
525  case PSK8:
526  amp_max = 1;
527  // EN 302 307, section 5.4.2
528  nrotations = 8;
529  nsymbols = 8;
530  symbols = new complex<signed char>[nsymbols];
531  symbols[0] = polar(1, 8, 1);
532  symbols[1] = polar(1, 8, 0);
533  symbols[2] = polar(1, 8, 4);
534  symbols[3] = polar(1, 8, 5);
535  symbols[4] = polar(1, 8, 2);
536  symbols[5] = polar(1, 8, 7);
537  symbols[6] = polar(1, 8, 3);
538  symbols[7] = polar(1, 8, 6);
539  make_lut_from_symbols(mer);
540  printf("cstln_lut::cstln_lut: PSK8\n");
541  break;
542  case APSK16:
543  {
544  // Default gamma for non-DVB-S2 applications.
545  if (gamma1 == 0)
546  gamma1 = 2.57;
547  // EN 302 307, section 5.4.3
548  float r1 = sqrtf(4.0f / (1.0f + 3.0f * gamma1 * gamma1));
549  float r2 = gamma1 * r1;
550  amp_max = r2;
551  nrotations = 4;
552  nsymbols = 16;
553  symbols = new complex<signed char>[nsymbols];
554  symbols[0] = polar(r2, 12, 1.5);
555  symbols[1] = polar(r2, 12, 10.5);
556  symbols[2] = polar(r2, 12, 4.5);
557  symbols[3] = polar(r2, 12, 7.5);
558  symbols[4] = polar(r2, 12, 0.5);
559  symbols[5] = polar(r2, 12, 11.5);
560  symbols[6] = polar(r2, 12, 5.5);
561  symbols[7] = polar(r2, 12, 6.5);
562  symbols[8] = polar(r2, 12, 2.5);
563  symbols[9] = polar(r2, 12, 9.5);
564  symbols[10] = polar(r2, 12, 3.5);
565  symbols[11] = polar(r2, 12, 8.5);
566  symbols[12] = polar(r1, 4, 0.5);
567  symbols[13] = polar(r1, 4, 3.5);
568  symbols[14] = polar(r1, 4, 1.5);
569  symbols[15] = polar(r1, 4, 2.5);
570  make_lut_from_symbols(mer);
571  printf("cstln_lut::cstln_lut: APSK16: gamma1=%f r1=%f r2=%f\n", gamma1, r1, r2);
572  break;
573  }
574  case APSK32:
575  {
576  // Default gammas for non-DVB-S2 applications.
577  if (gamma1 == 0)
578  gamma1 = 2.53;
579  if (gamma2 == 0)
580  gamma2 = 4.30;
581  // EN 302 307, section 5.4.3
582  float r1 = sqrtf(
583  8.0f / (1.0f + 3.0f * gamma1 * gamma1 + 4 * gamma2 * gamma2));
584  float r2 = gamma1 * r1;
585  float r3 = gamma2 * r1;
586  amp_max = r3;
587  nrotations = 4;
588  nsymbols = 32;
589  symbols = new complex<signed char>[nsymbols];
590  symbols[0] = polar(r2, 12, 1.5);
591  symbols[1] = polar(r2, 12, 2.5);
592  symbols[2] = polar(r2, 12, 10.5);
593  symbols[3] = polar(r2, 12, 9.5);
594  symbols[4] = polar(r2, 12, 4.5);
595  symbols[5] = polar(r2, 12, 3.5);
596  symbols[6] = polar(r2, 12, 7.5);
597  symbols[7] = polar(r2, 12, 8.5);
598  symbols[8] = polar(r3, 16, 1);
599  symbols[9] = polar(r3, 16, 3);
600  symbols[10] = polar(r3, 16, 14);
601  symbols[11] = polar(r3, 16, 12);
602  symbols[12] = polar(r3, 16, 6);
603  symbols[13] = polar(r3, 16, 4);
604  symbols[14] = polar(r3, 16, 9);
605  symbols[15] = polar(r3, 16, 11);
606  symbols[16] = polar(r2, 12, 0.5);
607  symbols[17] = polar(r1, 4, 0.5);
608  symbols[18] = polar(r2, 12, 11.5);
609  symbols[19] = polar(r1, 4, 3.5);
610  symbols[20] = polar(r2, 12, 5.5);
611  symbols[21] = polar(r1, 4, 1.5);
612  symbols[22] = polar(r2, 12, 6.5);
613  symbols[23] = polar(r1, 4, 2.5);
614  symbols[24] = polar(r3, 16, 0);
615  symbols[25] = polar(r3, 16, 2);
616  symbols[26] = polar(r3, 16, 15);
617  symbols[27] = polar(r3, 16, 13);
618  symbols[28] = polar(r3, 16, 7);
619  symbols[29] = polar(r3, 16, 5);
620  symbols[30] = polar(r3, 16, 8);
621  symbols[31] = polar(r3, 16, 10);
622  make_lut_from_symbols(mer);
623  printf("cstln_lut::cstln_lut: APSK32: gamma1=%f gamma2=%f, r1=%f r2=%f r3=%f\n", gamma1, gamma2, r1, r2, r3);
624  break;
625  }
626  case APSK64E:
627  {
628  // Default gammas for non-DVB-S2 applications.
629  if (gamma1 == 0)
630  gamma1 = 2.4;
631  if (gamma2 == 0)
632  gamma2 = 4.3;
633  if (gamma3 == 0)
634  gamma3 = 7.0;
635  // EN 302 307-2, section 5.4.5, Table 13e
636  float r1 = sqrtf(
637  64.0f / (4.0f + 12.0f * gamma1 * gamma1 + 20.0f * gamma2 * gamma2 + 28.0f * gamma3 * gamma3));
638  float r2 = gamma1 * r1;
639  float r3 = gamma2 * r1;
640  float r4 = gamma3 * r1;
641  amp_max = r4;
642  nrotations = 4;
643  nsymbols = 64;
644  symbols = new complex<signed char>[nsymbols];
645  polar2(0, r4, 1.0 / 4, 7.0 / 4, 3.0 / 4, 5.0 / 4);
646  polar2(4, r4, 13.0 / 28, 43.0 / 28, 15.0 / 28, 41.0 / 28);
647  polar2(8, r4, 1.0 / 28, 55.0 / 28, 27.0 / 28, 29.0 / 28);
648  polar2(12, r1, 1.0 / 4, 7.0 / 4, 3.0 / 4, 5.0 / 4);
649  polar2(16, r4, 9.0 / 28, 47.0 / 28, 19.0 / 28, 37.0 / 28);
650  polar2(20, r4, 11.0 / 28, 45.0 / 28, 17.0 / 28, 39.0 / 28);
651  polar2(24, r3, 1.0 / 20, 39.0 / 20, 19.0 / 20, 21.0 / 20);
652  polar2(28, r2, 1.0 / 12, 23.0 / 12, 11.0 / 12, 13.0 / 12);
653  polar2(32, r4, 5.0 / 28, 51.0 / 28, 23.0 / 28, 33.0 / 28);
654  polar2(36, r3, 9.0 / 20, 31.0 / 20, 11.0 / 20, 29.0 / 20);
655  polar2(40, r4, 3.0 / 28, 53.0 / 28, 25.0 / 28, 31.0 / 28);
656  polar2(44, r2, 5.0 / 12, 19.0 / 12, 7.0 / 12, 17.0 / 12);
657  polar2(48, r3, 1.0 / 4, 7.0 / 4, 3.0 / 4, 5.0 / 4);
658  polar2(52, r3, 7.0 / 20, 33.0 / 20, 13.0 / 20, 27.0 / 20);
659  polar2(56, r3, 3.0 / 20, 37.0 / 20, 17.0 / 20, 23.0 / 20);
660  polar2(60, r2, 1.0 / 4, 7.0 / 4, 3.0 / 4, 5.0 / 4);
661  make_lut_from_symbols(mer);
662  printf("cstln_lut::cstln_lut: APSK64E: gamma1=%f gamma2=%f, gamm3=%f r1=%f r2=%f r3=%f r4=%f\n", gamma1, gamma2, gamma3, r1, r2, r3, r4);
663  break;
664  }
665  case QAM16:
666  amp_max = 0;
667  make_qam(16);
668  break;
669  case QAM64:
670  amp_max = 1;
671  make_qam(64);
672  break;
673  case QAM256:
674  amp_max = 1;
675  make_qam(256);
676  break;
677  default:
678  fail("Constellation not implemented");
679  }
680  }
681 
682  struct result
683  {
684  SOFTSYMB ss;
685  s_angle phase_error;
686  uint8_t symbol; // Nearest symbol, useful for C&T recovery
687  };
688 
689  inline result *lookup(float I, float Q)
690  {
691  // Handling of overflows beyond the lookup table:
692  // - For BPSK/QPSK/8PSK we only care about the phase,
693  // so the following is harmless and improves locking at low SNR.
694  // - For amplitude modulations this is not appropriate.
695  // However, if there is enough noise to cause overflow,
696  // demodulation would probably fail anyway.
697  //
698  // Comment-out for better throughput at high SNR.
699 #if 1
700  while (I < -128 || I > 127 || Q < -128 || Q > 127)
701  {
702  I *= 0.5;
703  Q *= 0.5;
704  }
705 #endif
706  return &lut[(u8)(s8)I][(u8)(s8)Q];
707  }
708 
709  inline result *lookup(int I, int Q)
710  {
711  // Ignore wrapping modulo 256
712  return &lut[(u8)I][(u8)Q];
713  }
714 
715  private:
716  complex<signed char> polar(float r, int n, float i)
717  {
718  float a = i * 2 * M_PI / n;
719  return complex<signed char>(r * cosf(a) * cstln_amp,
720  r * sinf(a) * cstln_amp);
721  }
722 
723  // Helper function for some constellation tables
724  void polar2(int i, float r, float a0, float a1, float a2, float a3)
725  {
726  float a[] = {a0, a1, a2, a3};
727 
728  for (int j = 0; j < 4; ++j)
729  {
730  float phi = a[j] * M_PI;
731  symbols[i + j] = complex<signed char>(r * cosf(phi) * cstln_amp,
732  r * sinf(phi) * cstln_amp);
733  }
734  }
735 
736  void make_qam(int n)
737  {
738  nrotations = 4;
739  nsymbols = n;
740  symbols = new complex<signed char>[nsymbols];
741  int m = sqrtl(n);
742  float scale;
743 
744  { // Average power in first quadrant with unit grid
745  int q = m / 2;
746  float avgpower = 2 * (q * 0.25 + (q - 1) * q / 2.0 + (q - 1) * q * (2 * q - 1) / 6.0) / q;
747  scale = 1.0 / sqrtf(avgpower);
748  }
749  // Arbitrary mapping
750 
751  int s = 0;
752 
753  for (int x = 0; x < m; ++x)
754  {
755  for (int y = 0; y < m; ++y)
756  {
757  float I = x - (float)(m - 1) / 2;
758  float Q = y - (float)(m - 1) / 2;
759  symbols[s].re = I * scale * cstln_amp;
760  symbols[s].im = Q * scale * cstln_amp;
761  ++s;
762  }
763  }
764 
765  make_lut_from_symbols(20); // TBD
766  }
767 
768  result lut[R][R];
769 
770  void make_lut_from_symbols(float mer)
771  {
772  // Note: Excessively low values of MER will break 16APSK and 32APSK.
773  float sigma = cstln_amp * pow(10.0, (-mer / 20));
774 
775  // Precomputed values.
776  // Shared scope so that we don't have to reset dists2[nsymbols..] to -1.
777  struct full_ss fss;
778 
779  for (int s = 0; s < 256; ++s)
780  fss.dists2[s] = -1;
781 
782  for (int I = -R / 2; I < R / 2; ++I)
783  {
784  for (int Q = -R / 2; Q < R / 2; ++Q)
785  {
786  // Nearest symbol
787  fss.nearest = 0;
788  fss.dists2[0] = 65535;
789  // Conditional probabilities:
790  // Sum likelyhoods from all candidate symbols.
791  //
792  // P(TX[b]==B | RX==IQ) =
793  // sum(S=0..nsymbols-1, P(TX[b]==B | RX==IQ && TXs==S))
794  //
795  // P(TX[b] == B | RX==IQ && TX==S) =
796  // P(TX[b]==B && RX==IQ && TX==S) / P(RX==IQ && TX==S)
797  float probs[8][2];
798  memset(probs, 0, sizeof(probs));
799 
800  for (int s = 0; s < nsymbols; ++s)
801  {
802  float d2 = ((I - symbols[s].re) * (I - symbols[s].re) + (Q - symbols[s].im) * (Q - symbols[s].im));
803 
804  if (d2 < fss.dists2[fss.nearest])
805  fss.nearest = s;
806 
807  fss.dists2[s] = d2;
808  float p = expf(-d2 / (2 * sigma * sigma)) / (sqrtf(2 * M_PI) * sigma);
809 
810  for (int bit = 0; bit < 8; ++bit)
811  {
812  probs[bit][(s >> bit) & 1] += p;
813  }
814  }
815 
816  // Normalize
817  for (int b = 0; b < 8; ++b)
818  {
819  float p = probs[b][1] / (probs[b][0] + probs[b][1]);
820  // Avoid trouble when sigma is unrealistically low.
821  if (!isnormal(p))
822  p = 0;
823  fss.p[b] = p;
824  }
825 
826  result *pr = &lut[I & (R - 1)][Q & (R - 1)];
827  to_softsymb(&fss, &pr->ss);
828  // Always record nearest symbol and phase error for C&T.
829  pr->symbol = fss.nearest;
830  float ph_symbol = atan2f(symbols[pr->symbol].im,
831  symbols[pr->symbol].re);
832  float ph_err = atan2f(Q, I) - ph_symbol;
833  pr->phase_error = (int32_t)(ph_err * 65536 / (2 * M_PI)); // Mod 65536
834  }
835  }
836  }
837 
838  public:
839  void dump(FILE *f)
840  {
841  int bps = log2(nsymbols);
842  fprintf(f, "P5\n%d %d\n255\n", R, R * (bps + 1));
843 
844  for (int bit = 0; bit < bps + 1; ++bit)
845  {
846  for (int Q = R / 2 - 1; Q >= -R / 2; --Q)
847  {
848  for (int I = -R / 2; I < R / 2; ++I)
849  {
850  result *pr = &lut[I & (R - 1)][Q & (R - 1)];
851  uint8_t v;
852  if (bit < bps)
853  v = softsymb_to_dump(pr->ss, bit);
854  else
855  v = 128 + pr->phase_error / 64;
856  // Highlight the constellation symbols.
857  for (int s = 0; s < nsymbols; ++s)
858  {
859  if (symbols[s].re == I && symbols[s].im == Q)
860  v ^= 128;
861  }
862 
863  fputc(v, f);
864  }
865  }
866  }
867  }
868 
869  // Convert soft metric to Hamming distance
870  void harden()
871  {
872  for (int i = 0; i < R; ++i)
873  {
874  for (int q = 0; q < R; ++q)
875  softsymb_harden(&lut[i][q].ss);
876  }
877  }
878 
882 };
883 // cstln_lut
884 
885 // SAMPLER INTERFACE FOR CSTLN_RECEIVER
886 
887 template <typename T>
889 {
891  {
892  }
893  virtual complex<T> interp(const complex<T> *pin, float mu, float phase) = 0;
894  virtual void update_freq(float freqw, int period = 1)
895  {
896  (void) freqw;
897  (void) period;
898  } // 65536 = 1 Hz
899  virtual int readahead() = 0;
900 };
901 
902 // NEAREST-SAMPLE SAMPLER FOR CSTLN_RECEIVER
903 // Suitable for bandpass-filtered, oversampled signals only
904 
905 template <typename T>
907 {
908  int readahead()
909  {
910  return 0;
911  }
912 
913  complex<T> interp(const complex<T> *pin, float mu, float phase)
914  {
915  (void) mu;
916  return pin[0] * trig.expi(-phase);
917  }
918 
919  private:
921 };
922 // nearest_sampler
923 
924 // LINEAR SAMPLER FOR CSTLN_RECEIVER
925 
926 template <typename T>
928 {
929  int readahead()
930  {
931  return 1;
932  }
933 
934  complex<T> interp(const complex<T> *pin, float mu, float phase)
935  {
936  // Derotate pin[0] and pin[1]
937  complex<T> s0 = pin[0] * trig.expi(-phase);
938  complex<T> s1 = pin[1] * trig.expi(-(phase + freqw));
939  // Interpolate linearly
940  return s0 * (1 - mu) + s1 * mu;
941  }
942 
943  void update_freq(float _freqw, int period = 1)
944  {
945  (void) period;
946  freqw = _freqw;
947  }
948 
949  private:
951  float freqw;
952 };
953 // linear_sampler
954 
955 // FIR SAMPLER FOR CSTLN_RECEIVER
956 
957 template <typename T, typename Tc>
959 {
960  fir_sampler(int _ncoeffs, Tc *_coeffs, int _subsampling = 1) : ncoeffs(_ncoeffs),
961  coeffs(_coeffs),
962  subsampling(_subsampling),
963  shifted_coeffs(new complex<T>[ncoeffs]),
964  update_freq_phase(0)
965  {
966  do_update_freq(0); // In case application never calls update_freq()
967  }
968 
969  int readahead()
970  {
971  return ncoeffs - 1;
972  }
973 
974  complex<T> interp(const complex<T> *pin, float mu, float phase)
975  {
976  // Apply FIR filter with subsampling
977  complex<T> acc(0, 0);
978  complex<T> *pc = shifted_coeffs + (int)((1 - mu) * subsampling);
979  complex<T> *pcend = shifted_coeffs + ncoeffs;
980 
981  if (subsampling == 1)
982  {
983  // Special case for heavily oversampled signals,
984  // where filtering is expensive.
985  // gcc-4.9.2 can vectorize this form with NEON on ARM.
986  while (pc < pcend)
987  acc += (*pc++) * (*pin++);
988  }
989  else
990  {
991  // Not vectorized because the coefficients are not
992  // guaranteed to be contiguous in memory.
993  for (; pc < pcend; pc += subsampling, ++pin)
994  acc += (*pc) * (*pin);
995  }
996 
997  // Derotate
998  return trig.expi(-phase) * acc;
999  }
1000 
1001  void update_freq(float freqw, int period)
1002  {
1003  // Throttling: Update one coeff per 16 processed samples,
1004  // to keep the overhead of freq tracking below about 10%.
1005  update_freq_phase -= period;
1006 
1007  if (update_freq_phase <= 0)
1008  {
1009  update_freq_phase = ncoeffs * 16;
1010  do_update_freq(freqw);
1011  }
1012  }
1013 
1014  private:
1015  void do_update_freq(float freqw)
1016  {
1017  float f = freqw / subsampling;
1018  for (int i = 0; i < ncoeffs; ++i)
1019  shifted_coeffs[i] = trig.expi(-f * (i - ncoeffs / 2)) * coeffs[i];
1020  }
1021 
1023  int ncoeffs;
1024  Tc *coeffs;
1028 };
1029 // fir_sampler
1030 
1031 // CONSTELLATION RECEIVER
1032 
1033 // Linear interpolation: good enough for 1.2 samples/symbol,
1034 // but higher oversampling is recommended.
1035 
1036 template <typename T, typename SOFTSYMB>
1038 {
1041  unsigned long meas_decimation; // Measurement rate
1042  float omega, min_omega, max_omega; // Samples per symbol
1043  float freqw, min_freqw, max_freqw; // Freq offs (65536 = 1 Hz)
1045  bool allow_drift; // Follow carrier beyond safe limits
1046  static const unsigned int chunk_size = 128;
1047  float kest;
1048 
1050  sampler_interface<T> *_sampler,
1051  pipebuf<complex<T>> &_in,
1052  pipebuf<SOFTSYMB> &_out,
1053  pipebuf<float> *_freq_out = NULL,
1054  pipebuf<float> *_ss_out = NULL,
1055  pipebuf<float> *_mer_out = NULL,
1056  pipebuf<cf32> *_cstln_out = NULL) : runnable(sch, "Constellation receiver"),
1057  sampler(_sampler),
1058  cstln(NULL),
1059  meas_decimation(1048576),
1060  pll_adjustment(1.0),
1061  allow_drift(false),
1062  kest(0.01),
1063  in(_in),
1064  out(_out, chunk_size),
1065  est_insp(cstln_amp * cstln_amp),
1066  agc_gain(1),
1067  mu(0),
1068  phase(0),
1069  est_sp(0),
1070  est_ep(0),
1071  meas_count(0)
1072  {
1073  set_omega(1);
1074  set_freq(0);
1075  freq_out = _freq_out ? new pipewriter<float>(*_freq_out) : NULL;
1076  ss_out = _ss_out ? new pipewriter<float>(*_ss_out) : NULL;
1077  mer_out = _mer_out ? new pipewriter<float>(*_mer_out) : NULL;
1078  cstln_out = _cstln_out ? new pipewriter<cf32>(*_cstln_out) : NULL;
1079  memset(hist, 0, sizeof(hist));
1080  }
1081 
1082  void set_omega(float _omega, float tol = 10e-6)
1083  {
1084  omega = _omega;
1085  min_omega = omega * (1 - tol);
1086  max_omega = omega * (1 + tol);
1087  update_freq_limits();
1088  }
1089 
1090  void set_freq(float freq)
1091  {
1092  freqw = freq * 65536;
1093  update_freq_limits();
1094  refresh_freq_tap();
1095  }
1096 
1097  void set_allow_drift(bool d)
1098  {
1099  allow_drift = d;
1100  }
1101 
1103  {
1104  // Prevent PLL from crossing +-SR/n/2 and locking at +-SR/n.
1105  int n = 4;
1106 
1107  if (cstln)
1108  {
1109  switch (cstln->nsymbols)
1110  {
1111  case 2:
1112  n = 2;
1113  break; // BPSK
1114  case 4:
1115  n = 4;
1116  break; // QPSK
1117  case 8:
1118  n = 8;
1119  break; // 8PSK
1120  case 16:
1121  n = 12;
1122  break; // 16APSK
1123  case 32:
1124  n = 16;
1125  break; // 32APSK
1126  default:
1127  n = 4;
1128  break;
1129  }
1130  }
1131 
1132  min_freqw = freqw - 65536 / max_omega / n / 2;
1133  max_freqw = freqw + 65536 / max_omega / n / 2;
1134  }
1135 
1136  void run()
1137  {
1138  if (!cstln)
1139  fail("constellation not set");
1140 
1141  // Magic constants that work with the qa recordings.
1142  float freq_alpha = 0.04;
1143  float freq_beta = 0.0012 / omega * pll_adjustment;
1144  float gain_mu = 0.02 / (cstln_amp * cstln_amp) * 2;
1145 
1146  int max_meas = chunk_size / meas_decimation + 1;
1147  // Large margin on output_size because mu adjustments
1148  // can lead to more than chunk_size/min_omega symbols.
1149  while (in.readable() >= chunk_size + sampler->readahead() && out.writable() >= chunk_size && (!freq_out || freq_out->writable() >= max_meas) && (!ss_out || ss_out->writable() >= max_meas) && (!mer_out || mer_out->writable() >= max_meas) && (!cstln_out || cstln_out->writable() >= max_meas))
1150  {
1151  sampler->update_freq(freqw, chunk_size);
1152 
1153  complex<T> *pin = in.rd(), *pin0 = pin, *pend = pin + chunk_size;
1154  SOFTSYMB *pout = out.wr(), *pout0 = pout;
1155 
1156  // These are scoped outside the loop for SS and MER estimation.
1157  complex<float> sg{0.0f, 0.0f}; // Symbol before AGC;
1158  complex<float> s; // For MER estimation and constellation viewer
1159  complex<signed char> *cstln_point = NULL;
1160 
1161  while (pin < pend)
1162  {
1163  // Here mu is the time of the next symbol counted from 0 at pin.
1164  if (mu < 1)
1165  {
1166  // Here 0<=mu<1 is the fractional time of the next symbol
1167  // between pin and pin+1.
1168  sg = sampler->interp(pin, mu, phase + mu * freqw);
1169  s = sg * agc_gain;
1170 
1171  // Constellation look-up
1172  typename cstln_lut<SOFTSYMB, 256>::result *cr =
1173  cstln->lookup(s.re, s.im);
1174  *pout = cr->ss;
1175  ++pout;
1176 
1177  // PLL
1178  phase += cr->phase_error * freq_alpha;
1179  freqw += cr->phase_error * freq_beta;
1180 
1181  // Modified Mueller and Müller
1182  // mu[k]=real((c[k]-c[k-2])*conj(p[k-1])-(p[k]-p[k-2])*conj(c[k-1]))
1183  // =dot(c[k]-c[k-2],p[k-1]) - dot(p[k]-p[k-2],c[k-1])
1184  // p = received signals
1185  // c = decisions (constellation points)
1186  hist[2] = hist[1];
1187  hist[1] = hist[0];
1188  hist[0].p.re = s.re;
1189  hist[0].p.im = s.im;
1190  cstln_point = &cstln->symbols[cr->symbol];
1191  hist[0].c.re = cstln_point->re;
1192  hist[0].c.im = cstln_point->im;
1193  float muerr = ((hist[0].p.re - hist[2].p.re) * hist[1].c.re + (hist[0].p.im - hist[2].p.im) * hist[1].c.im) - ((hist[0].c.re - hist[2].c.re) * hist[1].p.re + (hist[0].c.im - hist[2].c.im) * hist[1].p.im);
1194  float mucorr = muerr * gain_mu;
1195  const float max_mucorr = 0.1;
1196  // TBD Optimize out statically
1197  if (mucorr < -max_mucorr)
1198  mucorr = -max_mucorr;
1199  if (mucorr > max_mucorr)
1200  mucorr = max_mucorr;
1201  mu += mucorr;
1202  mu += omega; // Next symbol time;
1203  } // mu<1
1204 
1205  // Next sample
1206  ++pin;
1207  --mu;
1208  phase += freqw;
1209  } // chunk_size
1210 
1211  in.read(pin - pin0);
1212  out.written(pout - pout0);
1213 
1214  // Normalize phase so that it never exceeds 32 bits.
1215  // Max freqw is 2^31/65536/chunk_size = 256 Hz
1216  // (this may happen with leandvb --drift --decim).
1217  phase = fmodf(phase, 65536); // Rounding direction irrelevant
1218 
1219  if (cstln_point)
1220  {
1221  // Output the last interpolated PSK symbol, max once per chunk_size
1222  if (cstln_out)
1223  cstln_out->write(s);
1224 
1225  // AGC
1226  // For APSK we must do AGC on the symbols, not the whole signal.
1227  // TODO Use a better estimator at low SNR.
1228  float insp = sg.re * sg.re + sg.im * sg.im;
1229  est_insp = insp * kest + est_insp * (1 - kest);
1230 
1231  if (est_insp)
1232  agc_gain = cstln_amp / gen_sqrt(est_insp);
1233 
1234  // SS and MER
1235  complex<float> ev(s.re - cstln_point->re,
1236  s.im - cstln_point->im);
1237  float sig_power, ev_power;
1238 
1239  if (cstln->nsymbols == 2)
1240  {
1241  // Special case for BPSK: Ignore quadrature component of noise.
1242  // TBD Projection on I axis assumes BPSK at 45°
1243  float sig_real = (cstln_point->re + cstln_point->im) * 0.707;
1244  float ev_real = (ev.re + ev.im) * 0.707;
1245  sig_power = sig_real * sig_real;
1246  ev_power = ev_real * ev_real;
1247  }
1248  else
1249  {
1250  sig_power = (int)cstln_point->re * cstln_point->re + (int)cstln_point->im * cstln_point->im;
1251  ev_power = ev.re * ev.re + ev.im * ev.im;
1252  }
1253 
1254  est_sp = sig_power * kest + est_sp * (1 - kest);
1255  est_ep = ev_power * kest + est_ep * (1 - kest);
1256  }
1257 
1258  // This is best done periodically ouside the inner loop,
1259  // but will cause non-deterministic output.
1260 
1261  if (!allow_drift)
1262  {
1263  if (freqw < min_freqw || freqw > max_freqw)
1264  freqw = (max_freqw + min_freqw) / 2;
1265  }
1266 
1267  // Output measurements
1268 
1269  refresh_freq_tap();
1270 
1271  meas_count += pin - pin0;
1272 
1273  while (meas_count >= meas_decimation)
1274  {
1275  meas_count -= meas_decimation;
1276  if (freq_out)
1277  freq_out->write(freq_tap);
1278  if (ss_out)
1279  ss_out->write(sqrtf(est_insp));
1280  if (mer_out)
1281  mer_out->write(
1282  est_ep ? 10 * logf(est_sp / est_ep) / logf(10) : 0);
1283  }
1284 
1285  } // Work to do
1286  }
1287 
1288  float freq_tap;
1289 
1291  {
1292  freq_tap = freqw / 65536;
1293  }
1294 
1295  private:
1296  struct
1297  {
1298  complex<float> p; // Received symbol
1299  complex<float> c; // Matched constellation point
1300  } hist[3];
1301 
1304  float est_insp, agc_gain;
1305  float mu; // PSK time expressed in clock ticks
1306  float phase; // 65536=2pi
1307  // Signal estimation
1308  float est_sp; // Estimated RMS signal power
1309  float est_ep; // Estimated RMS error vector power
1310  unsigned long meas_count;
1311  pipewriter<float> *freq_out, *ss_out, *mer_out;
1313 };
1314 
1315 // FAST QPSK RECEIVER
1316 
1317 // Optimized for u8 input, no AGC, uses phase information only.
1318 // Outputs hard symbols.
1319 
1320 template <typename T>
1322 {
1323  typedef u8 hardsymbol;
1324  unsigned long meas_decimation; // Measurement rate
1325  float omega, min_omega, max_omega; // Samples per symbol
1326  signed long freqw, min_freqw, max_freqw; // Freq offs (angle per sample)
1328  bool allow_drift; // Follow carrier beyond safe limits
1329  static const unsigned int chunk_size = 128;
1330 
1332  pipebuf<complex<T>> &_in,
1333  pipebuf<hardsymbol> &_out,
1334  pipebuf<float> *_freq_out = NULL,
1335  pipebuf<complex<T>> *_cstln_out = NULL) : runnable(sch, "Fast QPSK receiver"),
1336  meas_decimation(1048576),
1337  pll_adjustment(1.0),
1338  allow_drift(false),
1339  in(_in),
1340  out(_out, chunk_size),
1341  mu(0),
1342  phase(0),
1343  meas_count(0)
1344  {
1345  set_omega(1);
1346  set_freq(0);
1347  freq_out = _freq_out ? new pipewriter<float>(*_freq_out) : NULL;
1348  cstln_out = _cstln_out ? new pipewriter<complex<T>>(*_cstln_out) : NULL;
1349  memset(hist, 0, sizeof(hist));
1350  init_lookup_tables();
1351  }
1352 
1353  void set_omega(float _omega, float tol = 10e-6)
1354  {
1355  omega = _omega;
1356  min_omega = omega * (1 - tol);
1357  max_omega = omega * (1 + tol);
1358  update_freq_limits();
1359  }
1360 
1361  void set_freq(float freq)
1362  {
1363  freqw = freq * 65536;
1364  update_freq_limits();
1365  }
1366 
1368  {
1369  // Prevent PLL from locking at +-symbolrate/4.
1370  // TODO The +-SR/8 limit is suitable for QPSK only.
1371  min_freqw = freqw - 65536 / max_omega / 8;
1372  max_freqw = freqw + 65536 / max_omega / 8;
1373  }
1374 
1375  static const int RLUT_BITS = 8;
1376  static const int RLUT_ANGLES = 1 << RLUT_BITS;
1377 
1378  void run()
1379  {
1380  // Magic constants that work with the qa recordings.
1381  signed long freq_alpha = 0.04 * 65536;
1382  signed long freq_beta = 0.0012 * 256 * 65536 / omega * pll_adjustment;
1383 
1384  if (!freq_beta)
1385  fail("Excessive oversampling");
1386 
1387  float gain_mu = 0.02 / (cstln_amp * cstln_amp) * 2;
1388 
1389  int max_meas = chunk_size / meas_decimation + 1;
1390  // Largin margin on output_size because mu adjustments
1391  // can lead to more than chunk_size/min_omega symbols.
1392 
1393  while (in.readable() >= chunk_size + 1 && // +1 for interpolation
1394  out.writable() >= chunk_size && (!freq_out || freq_out->writable() >= max_meas) && (!cstln_out || cstln_out->writable() >= max_meas))
1395  {
1396  complex<T> *pin = in.rd(), *pin0 = pin, *pend = pin + chunk_size;
1397  hardsymbol *pout = out.wr(), *pout0 = pout;
1398 
1399  cu8 s;
1400  u_angle symbol_arg = 0; // Exported for constellation viewer
1401 
1402  while (pin < pend)
1403  {
1404  // Here mu is the time of the next symbol counted from 0 at pin.
1405  if (mu < 1)
1406  {
1407  // Here 0<=mu<1 is the fractional time of the next symbol
1408  // between pin and pin+1.
1409 
1410  // Derotate and interpolate
1411 #if 0 /* Phase only (does not work)
1412  Careful with the float/signed/unsigned casts */
1413  u_angle a0 = fast_arg(pin[0]) - phase;
1414  u_angle a1 = fast_arg(pin[1]) - (phase+freqw);
1415  s_angle da = a1 - a0;
1416  symbol_arg = a0 + (s_angle)(da*mu);
1417  s = arg_to_symbol(symbol_arg);
1418 #elif 1 // Linear by lookup-table. 1.2M on bench3bishs
1419  polar *p0 = &lut_polar[pin[0].re][pin[0].im];
1420  u_angle a0 = (u_angle)(p0->a - phase) >> (16 - RLUT_BITS);
1421  cu8 *p0r = &lut_rect[a0][p0->r >> 1];
1422  polar *p1 = &lut_polar[pin[1].re][pin[1].im];
1423  u_angle a1 = (u_angle)(p1->a - (phase + freqw)) >> (16 - RLUT_BITS);
1424  cu8 *p1r = &lut_rect[a1][p1->r >> 1];
1425  s.re = (int)(p0r->re + (p1r->re - p0r->re) * mu);
1426  s.im = (int)(p0r->im + (p1r->im - p0r->im) * mu);
1427  symbol_arg = fast_arg(s);
1428 #else // Linear floating-point, for reference
1429  float a0 = -(int)phase * M_PI / 32768;
1430  float cosa0 = cosf(a0), sina0 = sinf(a0);
1432  p0r(((float)pin[0].re - 128) * cosa0 - ((float)pin[0].im - 128) * sina0,
1433  ((float)pin[0].re - 128) * sina0 + ((float)pin[0].im - 128) * cosa0);
1434  float a1 = -(int)(phase + freqw) * M_PI / 32768;
1435  float cosa1 = cosf(a1), sina1 = sinf(a1);
1437  p1r(((float)pin[1].re - 128) * cosa1 - ((float)pin[1].im - 128) * sina1,
1438  ((float)pin[1].re - 128) * sina1 + ((float)pin[1].im - 128) * cosa1);
1439  s.re = (int)(128 + p0r.re + (p1r.re - p0r.re) * mu);
1440  s.im = (int)(128 + p0r.im + (p1r.im - p0r.im) * mu);
1441  symbol_arg = fast_arg(s);
1442 #endif
1443 
1444  int quadrant = symbol_arg >> 14;
1445  static unsigned char quadrant_to_symbol[4] = {0, 2, 3, 1};
1446  *pout = quadrant_to_symbol[quadrant];
1447  ++pout;
1448 
1449  // PLL
1450  s_angle phase_error = (s_angle)(symbol_arg & 16383) - 8192;
1451  phase += (phase_error * freq_alpha + 32768) >> 16;
1452  freqw += (phase_error * freq_beta + 32768 * 256) >> 24;
1453 
1454  // Modified Mueller and Müller
1455  // mu[k]=real((c[k]-c[k-2])*conj(p[k-1])-(p[k]-p[k-2])*conj(c[k-1]))
1456  // =dot(c[k]-c[k-2],p[k-1]) - dot(p[k]-p[k-2],c[k-1])
1457  // p = received signals
1458  // c = decisions (constellation points)
1459  hist[2] = hist[1];
1460  hist[1] = hist[0];
1461 #define HIST_FLOAT 0
1462 #if HIST_FLOAT
1463  hist[0].p.re = (float)s.re - 128;
1464  hist[0].p.im = (float)s.im - 128;
1465 
1466  cu8 cp = arg_to_symbol((symbol_arg & 49152) + 8192);
1467  hist[0].c.re = (float)cp.re - 128;
1468  hist[0].c.im = (float)cp.im - 128;
1469 
1470  float muerr =
1471  ((hist[0].p.re - hist[2].p.re) * hist[1].c.re +
1472  (hist[0].p.im - hist[2].p.im) * hist[1].c.im) -
1473  ((hist[0].c.re - hist[2].c.re) * hist[1].p.re +
1474  (hist[0].c.im - hist[2].c.im) * hist[1].p.im);
1475 #else
1476  hist[0].p = s;
1477  hist[0].c = arg_to_symbol((symbol_arg & 49152) + 8192);
1478 
1479  int muerr =
1480  ((signed char)(hist[0].p.re - hist[2].p.re) * ((int)hist[1].c.re - 128) + (signed char)(hist[0].p.im - hist[2].p.im) * ((int)hist[1].c.im - 128)) - ((signed char)(hist[0].c.re - hist[2].c.re) * ((int)hist[1].p.re - 128) + (signed char)(hist[0].c.im - hist[2].c.im) * ((int)hist[1].p.im - 128));
1481 #endif
1482  float mucorr = muerr * gain_mu;
1483  const float max_mucorr = 0.1;
1484  // TBD Optimize out statically
1485  if (mucorr < -max_mucorr)
1486  mucorr = -max_mucorr;
1487  if (mucorr > max_mucorr)
1488  mucorr = max_mucorr;
1489  mu += mucorr;
1490  mu += omega; // Next symbol time;
1491  } // mu<1
1492 
1493  // Next sample
1494  ++pin;
1495  --mu;
1496  phase += freqw;
1497  } // chunk_size
1498 
1499  in.read(pin - pin0);
1500  out.written(pout - pout0);
1501 
1502  if (symbol_arg && cstln_out)
1503  // Output the last interpolated PSK symbol, max once per chunk_size
1504  cstln_out->write(s);
1505 
1506  // This is best done periodically ouside the inner loop,
1507  // but will cause non-deterministic output.
1508 
1509  if (!allow_drift)
1510  {
1511  if (freqw < min_freqw || freqw > max_freqw)
1512  freqw = (max_freqw + min_freqw) / 2;
1513  }
1514 
1515  // Output measurements
1516 
1517  meas_count += pin - pin0;
1518 
1519  while (meas_count >= meas_decimation)
1520  {
1521  meas_count -= meas_decimation;
1522 
1523  if (freq_out)
1524  freq_out->write((float)freqw / 65536);
1525  }
1526 
1527  } // Work to do
1528  }
1529 
1530  private:
1531  struct polar
1532  {
1533  u_angle a;
1534  unsigned char r;
1535  } lut_polar[256][256];
1536 
1537  u_angle fast_arg(const cu8 &c)
1538  {
1539  // TBD read cu8 as u16 index, same endianness as in init()
1540  return lut_polar[c.re][c.im].a;
1541  }
1542 
1543  cu8 lut_rect[RLUT_ANGLES][256];
1544  cu8 lut_sincos[65536];
1545 
1546  cu8 arg_to_symbol(u_angle a)
1547  {
1548  return lut_sincos[a];
1549  }
1550 
1552  {
1553  for (int i = 0; i < 256; ++i)
1554  {
1555  for (int q = 0; q < 256; ++q)
1556  {
1557  // Don't cast float to unsigned directly
1558  lut_polar[i][q].a = (s_angle)(atan2f(q - 128, i - 128) * 65536 / (2 * M_PI));
1559  lut_polar[i][q].r = (int)hypotf(i - 128, q - 128);
1560  }
1561  }
1562 
1563  for (unsigned long a = 0; a < 65536; ++a)
1564  {
1565  float f = 2 * M_PI * a / 65536;
1566  lut_sincos[a].re = 128 + cstln_amp * cosf(f);
1567  lut_sincos[a].im = 128 + cstln_amp * sinf(f);
1568  }
1569 
1570  for (int a = 0; a < RLUT_ANGLES; ++a)
1571  {
1572  for (int r = 0; r < 256; ++r)
1573  {
1574  lut_rect[a][r].re = (int)(128 + r * cos(2 * M_PI * a / RLUT_ANGLES));
1575  lut_rect[a][r].im = (int)(128 + r * sin(2 * M_PI * a / RLUT_ANGLES));
1576  }
1577  }
1578  }
1579 
1580  struct
1581  {
1582 #if HIST_FLOAT
1583  complex<float> p; // Received symbol
1584  complex<float> c; // Matched constellation point
1585 #else
1586  cu8 p; // Received symbol
1587  cu8 c; // Matched constellation point
1588 #endif
1589  } hist[3];
1590 
1593  float mu; // PSK time expressed in clock ticks. TBD fixed point.
1594  u_angle phase;
1595  unsigned long meas_count;
1598 };
1599 // fast_qpsk_receiver
1600 
1601 // CONSTELLATION TRANSMITTER
1602 
1603 // Maps symbols to I/Q points.
1604 
1605 template <typename Tout, int Zout>
1607 {
1609 
1611  pipebuf<u8> &_in,
1612  pipebuf<complex<Tout>> &_out) : runnable(sch, "cstln_transmitter"),
1613  in(_in),
1614  out(_out),
1615  cstln(0)
1616  {
1617  }
1618 
1619  void run()
1620  {
1621  if (!cstln)
1622  fail("constellation not set");
1623 
1624  int count = min(in.readable(), out.writable());
1625  u8 *pin = in.rd(), *pend = pin + count;
1626  complex<Tout> *pout = out.wr();
1627 
1628  for (; pin < pend; ++pin, ++pout)
1629  {
1630  complex<signed char> *cp = &cstln->symbols[*pin];
1631  pout->re = Zout + cp->re;
1632  pout->im = Zout + cp->im;
1633  }
1634 
1635  in.read(count);
1636  out.written(count);
1637  }
1638 
1639  private:
1642 };
1643 // cstln_transmitter
1644 
1645 // FREQUENCY SHIFTER
1646 
1647 // Resolution is sample_freq/65536.
1648 
1649 template <typename T>
1651 {
1653  pipebuf<complex<T>> &_in,
1654  pipebuf<complex<T>> &_out,
1655  float freq) : runnable(sch, "rotator"),
1656  in(_in),
1657  out(_out),
1658  index(0)
1659  {
1660  int ifreq = freq * 65536;
1661  if (sch->debug)
1662  fprintf(stderr, "Rotate: req=%f real=%f\n", freq, ifreq / 65536.0);
1663 
1664  for (int i = 0; i < 65536; ++i)
1665  {
1666  lut_cos[i] = cosf(2 * M_PI * i * ifreq / 65536);
1667  lut_sin[i] = sinf(2 * M_PI * i * ifreq / 65536);
1668  }
1669  }
1670 
1671  void run()
1672  {
1673  unsigned long count = min(in.readable(), out.writable());
1674  complex<T> *pin = in.rd(), *pend = pin + count;
1675  complex<T> *pout = out.wr();
1676 
1677  for (; pin < pend; ++pin, ++pout, ++index)
1678  {
1679  float c = lut_cos[index];
1680  float s = lut_sin[index];
1681  pout->re = pin->re * c - pin->im * s;
1682  pout->im = pin->re * s + pin->im * c;
1683  }
1684 
1685  in.read(count);
1686  out.written(count);
1687  }
1688 
1689  private:
1692  float lut_cos[65536];
1693  float lut_sin[65536];
1694  unsigned short index; // Current phase
1695 };
1696 // rotator
1697 
1698 // SPECTRUM-BASED CNR ESTIMATOR
1699 
1700 // Assumes that the spectrum is as follows:
1701 //
1702 // ---|--noise---|-roll-off-|---carrier+noise----|-roll-off-|---noise--|---
1703 // | (bw/2) | (bw) | (bw/2) | (bw) | (bw/2) |
1704 //
1705 // Maximum roll-off 0.5
1706 
1707 template <typename T>
1709 {
1711  pipebuf<complex<T>> &_in,
1712  pipebuf<float> &_out,
1713  float _bandwidth, int nfft = 4096) : runnable(sch, "cnr_fft"),
1714  bandwidth(_bandwidth),
1715  freq_tap(NULL),
1716  tap_multiplier(1),
1717  decimation(1048576),
1718  kavg(0.1),
1719  in(_in),
1720  out(_out),
1721  fft(nfft),
1722  avgpower(NULL),
1723  phase(0)
1724  {
1725  if (bandwidth > 0.25)
1726  fail("CNR estimator requires Fsampling > 4x Fsignal");
1727  }
1728 
1729  float bandwidth;
1730  float *freq_tap, tap_multiplier;
1732  float kavg;
1733 
1734  void run()
1735  {
1736  while (in.readable() >= fft.n && out.writable() >= 1)
1737  {
1738  phase += fft.n;
1739 
1740  if (phase >= decimation)
1741  {
1742  phase -= decimation;
1743  do_cnr();
1744  }
1745 
1746  in.read(fft.n);
1747  }
1748  }
1749 
1750  private:
1751  void do_cnr()
1752  {
1753  float center_freq = freq_tap ? *freq_tap * tap_multiplier : 0;
1754  int icf = floor(center_freq * fft.n + 0.5);
1755  complex<T> *data = new complex<T>[fft.n];
1756  memcpy(data, in.rd(), fft.n * sizeof(data[0]));
1757  fft.inplace(data, true);
1758  T *power = new T[fft.n];
1759 
1760  for (int i = 0; i < fft.n; ++i)
1761  power[i] = data[i].re * data[i].re + data[i].im * data[i].im;
1762 
1763  if (!avgpower)
1764  {
1765  // Initialize with first spectrum
1766  avgpower = new T[fft.n];
1767  memcpy(avgpower, power, fft.n * sizeof(avgpower[0]));
1768  }
1769 
1770  // Accumulate and low-pass filter
1771  for (int i = 0; i < fft.n; ++i)
1772  avgpower[i] = avgpower[i] * (1 - kavg) + power[i] * kavg;
1773 
1774  int bwslots = (bandwidth / 4) * fft.n;
1775 
1776  if (!bwslots)
1777  {
1778  delete[] power;
1779  delete[] data;
1780  return;
1781  }
1782 
1783  // Measure carrier+noise in center band
1784  float c2plusn2 = avgslots(icf - bwslots, icf + bwslots);
1785  // Measure noise left and right of roll-off zones
1786  float n2 = (avgslots(icf - bwslots * 4, icf - bwslots * 3) + avgslots(icf + bwslots * 3, icf + bwslots * 4)) / 2;
1787  float c2 = c2plusn2 - n2;
1788  float cnr = (c2 > 0 && n2 > 0) ? 10 * logf(c2 / n2) / logf(10) : -50;
1789  out.write(cnr);
1790  delete[] power;
1791  delete[] data;
1792  }
1793 
1794  float avgslots(int i0, int i1)
1795  { // i0 <= i1
1796  T s = 0;
1797 
1798  for (int i = i0; i <= i1; ++i)
1799  s += avgpower[i & (fft.n - 1)];
1800 
1801  return s / (i1 - i0 + 1);
1802  }
1803 
1808  int phase;
1809 };
1810 // cnr_fft
1811 
1812 template <typename T, int NFFT>
1814 {
1816  pipebuf<complex<T>> &_in,
1817  pipebuf<float[NFFT]> &_out) : runnable(sch, "spectrum"),
1818  decimation(1048576),
1819  kavg(0.1),
1820  decim(1), in(_in),
1821  out(_out),
1822  fft(NFFT),
1823  avgpower(NULL),
1824  phase(0)
1825  {
1826  }
1827 
1829  float kavg;
1830  int decim;
1831 
1832  void run()
1833  {
1834  while (in.readable() >= fft.n * decim && out.writable() >= 1)
1835  {
1836  phase += fft.n * decim;
1837 
1838  if (phase >= decimation)
1839  {
1840  phase -= decimation;
1841  do_spectrum();
1842  }
1843 
1844  in.read(fft.n * decim);
1845  }
1846  }
1847 
1848  private:
1850  {
1851  complex<T> data[fft.n];
1852 
1853  if (decim == 1)
1854  {
1855  memcpy(data, in.rd(), fft.n * sizeof(data[0]));
1856  }
1857  else
1858  {
1859  complex<T> *pin = in.rd();
1860 
1861  for (int i = 0; i < fft.n; ++i, pin += decim)
1862  data[i] = *pin;
1863  }
1864 
1865  fft.inplace(data, true);
1866  float power[NFFT];
1867 
1868  for (int i = 0; i < fft.n; ++i)
1869  power[i] = (float)data[i].re * data[i].re + (float)data[i].im * data[i].im;
1870 
1871  if (!avgpower)
1872  {
1873  // Initialize with first spectrum
1874  avgpower = new float[fft.n];
1875  memcpy(avgpower, power, fft.n * sizeof(avgpower[0]));
1876  }
1877 
1878  // Accumulate and low-pass filter
1879  for (int i = 0; i < fft.n; ++i)
1880  avgpower[i] = avgpower[i] * (1 - kavg) + power[i] * kavg;
1881 
1882  // Reuse power[]
1883  for (int i = 0; i < fft.n / 2; ++i)
1884  {
1885  power[i] = 10 * log10f(avgpower[NFFT / 2 + i]);
1886  power[NFFT / 2 + i] = 10 * log10f(avgpower[i]);
1887  }
1888 
1889  memcpy(out.wr(), power, sizeof(power[0]) * NFFT);
1890  out.written(1);
1891  }
1892 
1897  int phase;
1898 };
1899 // spectrum
1900 
1901 } // namespace leansdr
1902 
1903 #endif // LEANSDR_SDR_H
short int16_t
Definition: rtptypes_win.h:43
float kavg
Definition: sdr.h:1732
int decimation
Definition: sdr.h:48
bool llr_harden(llr_t v)
Definition: sdr.h:446
unsigned long window_size
Definition: sdr.h:226
simple_agc(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< complex< T >> &_out)
Definition: sdr.h:335
cstln_receiver(scheduler *sch, sampler_interface< T > *_sampler, pipebuf< complex< T >> &_in, pipebuf< SOFTSYMB > &_out, pipebuf< float > *_freq_out=NULL, pipebuf< float > *_ss_out=NULL, pipebuf< float > *_mer_out=NULL, pipebuf< cf32 > *_cstln_out=NULL)
Definition: sdr.h:1049
spectrum(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< float[NFFT]> &_out)
Definition: sdr.h:1815
cstln_lut< SOFTSYMB, 256 > * cstln
Definition: sdr.h:1040
slot * __slots
Definition: sdr.h:213
uint16_t u_angle
Definition: sdr.h:384
signed long min_freqw
Definition: sdr.h:1326
pipewriter< cf32 > * cstln_out
Definition: sdr.h:1312
pipereader< complex< T > > in
Definition: sdr.h:320
ss_estimator(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< T > &_out)
Definition: sdr.h:229
Fixed< IntType, IntBits > cos(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2271
void set_omega(float _omega, float tol=10e-6)
Definition: sdr.h:1353
cfft_engine< T > fft
Definition: sdr.h:1806
int m_typeCode
Definition: sdr.h:879
unsigned long decimation
Definition: sdr.h:269
int update_freq_phase
Definition: sdr.h:1027
void to_softsymb(const full_ss *fss, hard_ss *ss)
Definition: sdr.cpp:55
complex< signed char > polar(float r, int n, float i)
Definition: sdr.h:716
complex< f32 > cf32
Definition: sdr.h:34
void set_freq(float freq)
Definition: sdr.h:1361
void run()
Definition: sdr.h:1671
fast_qpsk_receiver(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< hardsymbol > &_out, pipebuf< float > *_freq_out=NULL, pipebuf< complex< T >> *_cstln_out=NULL)
Definition: sdr.h:1331
void detect()
Definition: sdr.h:98
uint8_t nearest
Definition: sdr.h:409
void run()
Definition: sdr.h:350
fir_sampler(int _ncoeffs, Tc *_coeffs, int _subsampling=1)
Definition: sdr.h:960
float f32
Definition: sdr.h:28
void make_qam(int n)
Definition: sdr.h:736
pipereader< complex< T > > in
Definition: sdr.h:201
unsigned char u8
Definition: framework.h:453
void run()
Definition: sdr.h:1832
unsigned long meas_decimation
Definition: sdr.h:1041
cf32 * shifted_coeffs
Definition: sdr.h:1026
cstln_lut< hard_ss, 256 > * cstln
Definition: sdr.h:1608
void process()
Definition: sdr.h:169
uint16_t discr2
Definition: sdr.h:435
bool m_setByModcod
Definition: sdr.h:881
#define M_PI
Definition: rdsdemod.cpp:27
unsigned long decimation
Definition: sdr.h:227
complex< float > p
Definition: sdr.h:1298
complex< float > * ej
Definition: sdr.h:209
pipewriter< SOFTSYMB > out
Definition: sdr.h:1303
Fixed< IntType, IntBits > floor(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2301
pipewriter< T > out
Definition: sdr.h:261
pipewriter< complex< Tout > > out
Definition: sdr.h:1641
float kavg
Definition: sdr.h:1829
cfft_engine< T > fft
Definition: sdr.h:1895
pipewriter< float > * mer_out
Definition: sdr.h:1596
pipewriter< cu8 > * cstln_out
Definition: sdr.h:1597
void do_spectrum()
Definition: sdr.h:1849
pipewriter< float[NFFT]> out
Definition: sdr.h:1894
cnr_fft(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< float > &_out, float _bandwidth, int nfft=4096)
Definition: sdr.h:1710
int16_t s_angle
Definition: sdr.h:385
unsigned long window_size
Definition: sdr.h:268
pipereader< complex< T > > in
Definition: sdr.h:260
float amp_max
Definition: sdr.h:477
virtual complex< T > interp(const complex< T > *pin, float mu, float phase)=0
void make_lut_from_symbols(float mer)
Definition: sdr.h:770
complex< float > estim
Definition: sdr.h:208
void set_allow_drift(bool d)
Definition: sdr.h:1097
unsigned char uint8_t
Definition: rtptypes_win.h:42
void softsymb_harden(llr_ss *ss)
Definition: sdr.cpp:20
complex< u8 > cu8
Definition: sdr.h:30
pipewriter< float > * ss_out
Definition: sdr.h:1311
unsigned long meas_decimation
Definition: sdr.h:1324
int16_t cost
Definition: sdr.h:395
T gen_abs(T x)
uint8_t softsymb_to_dump(const llr_ss &ss, int bit)
Definition: sdr.cpp:38
unsigned short uint16_t
Definition: rtptypes_win.h:44
float out_rms
Definition: sdr.h:330
void fail(const char *s)
Definition: framework.cpp:11
pipereader< u8 > in
Definition: sdr.h:1640
rotator(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< complex< T >> &_out, float freq)
Definition: sdr.h:1652
pipereader< complex< T > > in
Definition: sdr.h:1690
uint16_t dists2[256]
Definition: sdr.h:410
unsigned long meas_count
Definition: sdr.h:1595
complex< T > interp(const complex< T > *pin, float mu, float phase)
Definition: sdr.h:934
result * lookup(float I, float Q)
Definition: sdr.h:689
void run()
Definition: sdr.h:80
Fixed< IntType, IntBits > sin(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2265
int32_t i
Definition: decimators.h:244
pipewriter< T > out_ss
Definition: sdr.h:321
char int8_t
Definition: rtptypes_win.h:41
cstln_lut(cstln_base::predef type, float mer=10, float gamma1=0, float gamma2=0, float gamma3=0)
Definition: sdr.h:487
pipereader< complex< T > > in
Definition: sdr.h:347
complex< T > interp(const complex< T > *pin, float mu, float phase)
Definition: sdr.h:913
complex< float > c
Definition: sdr.h:1299
int decimation
Definition: sdr.h:1828
uint8_t hard_ss
Definition: sdr.h:423
pipewriter< complex< T > > out
Definition: sdr.h:348
float p[8]
Definition: sdr.h:411
pipewriter< float > out
Definition: sdr.h:1805
int int32_t
Definition: rtptypes_win.h:45
u_angle fast_arg(const cu8 &c)
Definition: sdr.h:1537
virtual void update_freq(float freqw, int period=1)
Definition: sdr.h:894
pipereader< complex< T > > in
Definition: sdr.h:1893
void written(unsigned long n)
Definition: framework.h:308
void set_omega(float _omega, float tol=10e-6)
Definition: sdr.h:1082
pipereader< complex< T > > in
Definition: sdr.h:1302
void polar2(int i, float r, float a0, float a1, float a2, float a3)
Definition: sdr.h:724
auto_notch(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< complex< T >> &_out, int _nslots, T _agc_rms_setpoint)
Definition: sdr.h:51
signed char s8
Definition: framework.h:456
float estimated
Definition: sdr.h:332
virtual ~sampler_interface()
Definition: sdr.h:890
uint8_t nearest
Definition: sdr.h:436
complex< s16 > cs16
Definition: sdr.h:33
pipewriter< complex< T > > out
Definition: sdr.h:202
float avgslots(int i0, int i1)
Definition: sdr.h:1794
int m_rateCode
Definition: sdr.h:880
sampler_interface< T > * sampler
Definition: sdr.h:1039
unsigned long phase
Definition: sdr.h:262
pipewriter< complex< T > > out
Definition: sdr.h:1691
void update_freq_limits()
Definition: sdr.h:1102
complex< s8 > cs8
Definition: sdr.h:31
T gen_sqrt(T x)
unsigned short index
Definition: sdr.h:1694
complex< float > * expj
Definition: sdr.h:209
void do_cnr()
Definition: sdr.h:1751
scheduler * sch
Definition: framework.h:199
result * lookup(int I, int Q)
Definition: sdr.h:709
cfft_engine< float > fft
Definition: sdr.h:200
std::complex< Fixed< IntType, IntBits > > polar(const Fixed< IntType, IntBits > &rho, const Fixed< IntType, IntBits > &theta)
Definition: fixed.h:2409
complex< u16 > cu16
Definition: sdr.h:32
void run()
Definition: sdr.h:1734
void update_freq(float freqw, int period)
Definition: sdr.h:1001
uint8_t i0
Definition: decimators.h:244
pipewriter< hardsymbol > out
Definition: sdr.h:1592
complex< int8_t > * symbols
Definition: sdr.h:478
float bandwidth
Definition: sdr.h:1729
const int n
Definition: dsp.h:64
void refresh_freq_tap()
Definition: sdr.h:1290
unsigned long meas_count
Definition: sdr.h:1310
const float cstln_amp
Definition: sdr.h:404
complex< T > interp(const complex< T > *pin, float mu, float phase)
Definition: sdr.h:974
cu8 arg_to_symbol(u_angle a)
Definition: sdr.h:1546
void harden()
Definition: sdr.h:870
int readahead()
Definition: sdr.h:969
void set_freq(float freq)
Definition: sdr.h:1090
virtual int readahead()=0
uint8_t symbol
Definition: sdr.h:396
unsigned long phase
Definition: sdr.h:322
ss_amp_estimator(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< T > &_out_ss, pipebuf< T > &_out_ampmin, pipebuf< T > &_out_ampmax)
Definition: sdr.h:271
void dump(FILE *f)
Definition: sdr.h:839
int8_t llr_t
Definition: sdr.h:444
pipereader< cu8 > in
Definition: sdr.h:1591
pipereader< complex< T > > in
Definition: sdr.h:1804
void do_update_freq(float freqw)
Definition: sdr.h:1015
void update_freq(float _freqw, int period=1)
Definition: sdr.h:943
uint8_t i1
Definition: decimators.h:245
T * avgpower
Definition: sdr.h:1807
cstln_transmitter(scheduler *sch, pipebuf< u8 > &_in, pipebuf< complex< Tout >> &_out)
Definition: sdr.h:1610
void read(unsigned long n)
Definition: framework.h:367
float tap_multiplier
Definition: sdr.h:1730
int decimation
Definition: sdr.h:1731
void inplace(complex< T > *data, bool reverse=false)
Definition: dsp.h:89
T min(const T &x, const T &y)
Definition: framework.h:440