55 T _agc_rms_setpoint) :
runnable(sch,
"auto_notch"),
56 decimation(1024 * 4096),
68 for (
int s = 0; s <
nslots; ++s)
86 if (
phase >= decimation)
102 float m0 = 0, m2 = 0;
104 for (
int i = 0;
i <
fft.
n; ++
i)
108 m2 += (float)pin[
i].re * pin[
i].re + (
float)pin[
i].
im * pin[
i].
im;
119 fprintf(stderr,
"(pow %f max %f)", rms, m0);
125 float *amp =
new float[
fft.
n];
127 for (
int i = 0;
i <
fft.
n; ++
i)
128 amp[
i] = hypotf(data[
i].re, data[
i].im);
134 for (
int i = 0;
i <
fft.
n; ++
i)
135 if (amp[
i] > amp[iamax])
141 fprintf(stderr,
"%s: slot %d new peak %d -> %d\n",
name, (
int)(s -
__slots), s->i, iamax);
148 for (
int i = 0;
i <
fft.
n; ++
i)
151 s->expj[
i].re = cosf(a);
152 s->expj[
i].im = sinf(a);
161 if (iamax + 1 <
fft.
n)
176 for (; pin < pend; ++pin, ++pout)
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);
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);
223 template <
typename T>
240 while (
in.readable() >= window_size &&
out.writable() >= 1)
242 phase += window_size;
244 if (
phase >= decimation)
250 for (; p < pend; ++p)
251 s += (
float)p->
re * p->
re + (float)p->
im * p->
im;
253 out.write(sqrtf(s / window_size));
255 in.read(window_size);
265 template <
typename T>
280 out_ampmin(_out_ampmin),
281 out_ampmax(_out_ampmax),
288 while (
in.readable() >= window_size && out_ss.writable() >= 1 && out_ampmin.writable() >= 1 && out_ampmax.writable() >= 1)
290 phase += window_size;
292 if (
phase >= decimation)
297 float amin = 1e38, amax = 0;
299 for (; p < pend; ++p)
301 float mag2 = (float)p->
re * p->
re + (
float)p->
im * p->
im;
303 float mag = sqrtf(mag2);
310 out_ss.write(sqrtf(s2 / window_size));
311 out_ampmin.write(amin);
312 out_ampmax.write(amax);
315 in.read(window_size);
327 template <
typename T>
333 static const int chunk_size = 128;
342 out(_out, chunk_size)
357 for (; pin < pend; ++pin)
358 amp2 += pin->
re * pin->
re + pin->
im * pin->
im;
365 estimated = estimated * (1 - bw) + amp2 * bw;
366 float gain = estimated ? out_rms / sqrtf(estimated) : 0;
369 float bwcomp = 1 - bw;
371 for (; pin < pend; ++pin, ++pout)
433 static const int MAX_SYMBOLS = 4;
476 static const char *names[];
484 template <
typename SOFTSYMB,
int R>
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");
509 make_lut_from_symbols(mer);
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");
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");
548 float r1 = sqrtf(4.0f / (1.0f + 3.0f * gamma1 * gamma1));
549 float r2 = gamma1 * r1;
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);
583 8.0f / (1.0f + 3.0f * gamma1 * gamma1 + 4 * gamma2 * gamma2));
584 float r2 = gamma1 * r1;
585 float r3 = gamma2 * r1;
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);
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;
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);
678 fail(
"Constellation not implemented");
700 while (I < -128 || I > 127 || Q < -128 || Q > 127)
712 return &lut[(
u8)I][(
u8)Q];
718 float a = i * 2 *
M_PI / n;
724 void polar2(
int i,
float r,
float a0,
float a1,
float a2,
float a3)
726 float a[] = {a0, a1, a2, a3};
728 for (
int j = 0; j < 4; ++j)
730 float phi = a[j] *
M_PI;
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);
753 for (
int x = 0; x < m; ++x)
755 for (
int y = 0; y < m; ++y)
757 float I = x - (float)(m - 1) / 2;
758 float Q = y - (float)(m - 1) / 2;
765 make_lut_from_symbols(20);
773 float sigma = cstln_amp * pow(10.0, (-mer / 20));
779 for (
int s = 0; s < 256; ++s)
782 for (
int I = -R / 2; I < R / 2; ++I)
784 for (
int Q = -R / 2; Q < R / 2; ++Q)
798 memset(probs, 0,
sizeof(probs));
800 for (
int s = 0; s < nsymbols; ++s)
802 float d2 = ((I - symbols[s].re) * (I - symbols[s].re) + (Q - symbols[s].im) * (Q - symbols[s].im));
808 float p = expf(-d2 / (2 * sigma * sigma)) / (sqrtf(2 *
M_PI) * sigma);
810 for (
int bit = 0; bit < 8; ++bit)
812 probs[bit][(s >> bit) & 1] += p;
817 for (
int b = 0; b < 8; ++b)
819 float p = probs[b][1] / (probs[b][0] + probs[b][1]);
826 result *pr = &lut[I & (R - 1)][Q & (R - 1)];
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));
841 int bps = log2(nsymbols);
842 fprintf(f,
"P5\n%d %d\n255\n", R, R * (bps + 1));
844 for (
int bit = 0; bit < bps + 1; ++bit)
846 for (
int Q = R / 2 - 1; Q >= -R / 2; --Q)
848 for (
int I = -R / 2; I < R / 2; ++I)
850 result *pr = &lut[I & (R - 1)][Q & (R - 1)];
855 v = 128 + pr->phase_error / 64;
857 for (
int s = 0; s < nsymbols; ++s)
859 if (symbols[s].re == I && symbols[s].im == Q)
872 for (
int i = 0;
i < R; ++
i)
874 for (
int q = 0; q < R; ++q)
887 template <
typename T>
899 virtual int readahead() = 0;
905 template <
typename T>
916 return pin[0] * trig.expi(-phase);
926 template <
typename T>
938 complex<T> s1 = pin[1] * trig.expi(-(phase + freqw));
940 return s0 * (1 - mu) + s1 * mu;
957 template <
typename T,
typename Tc>
960 fir_sampler(
int _ncoeffs, Tc *_coeffs,
int _subsampling = 1) : ncoeffs(_ncoeffs),
962 subsampling(_subsampling),
963 shifted_coeffs(new
complex<T>[ncoeffs]),
978 complex<T> *pc = shifted_coeffs + (int)((1 - mu) * subsampling);
981 if (subsampling == 1)
987 acc += (*pc++) * (*pin++);
993 for (; pc < pcend; pc += subsampling, ++pin)
994 acc += (*pc) * (*pin);
998 return trig.expi(-phase) * acc;
1005 update_freq_phase -= period;
1007 if (update_freq_phase <= 0)
1009 update_freq_phase = ncoeffs * 16;
1010 do_update_freq(freqw);
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];
1036 template <
typename T,
typename SOFTSYMB>
1046 static const unsigned int chunk_size = 128;
1059 meas_decimation(1048576),
1060 pll_adjustment(1.0),
1064 out(_out, chunk_size),
1065 est_insp(cstln_amp * cstln_amp),
1079 memset(hist, 0,
sizeof(hist));
1085 min_omega = omega * (1 - tol);
1086 max_omega = omega * (1 + tol);
1087 update_freq_limits();
1092 freqw = freq * 65536;
1093 update_freq_limits();
1132 min_freqw = freqw - 65536 / max_omega / n / 2;
1133 max_freqw = freqw + 65536 / max_omega / n / 2;
1139 fail(
"constellation not set");
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;
1146 int max_meas = chunk_size / meas_decimation + 1;
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))
1153 complex<T> *pin =
in.rd(), *pin0 = pin, *pend = pin + chunk_size;
1154 SOFTSYMB *pout =
out.wr(), *pout0 = pout;
1168 sg = sampler->
interp(pin, mu,
phase + mu * freqw);
1178 phase += cr->phase_error * freq_alpha;
1179 freqw += cr->phase_error * freq_beta;
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;
1197 if (mucorr < -max_mucorr)
1198 mucorr = -max_mucorr;
1199 if (mucorr > max_mucorr)
1200 mucorr = max_mucorr;
1211 in.read(pin - pin0);
1212 out.written(pout - pout0);
1223 cstln_out->write(s);
1228 float insp = sg.re * sg.re + sg.im * sg.im;
1229 est_insp = insp * kest + est_insp * (1 - kest);
1232 agc_gain = cstln_amp /
gen_sqrt(est_insp);
1236 s.
im - cstln_point->
im);
1237 float sig_power, ev_power;
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;
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;
1254 est_sp = sig_power * kest + est_sp * (1 - kest);
1255 est_ep = ev_power * kest + est_ep * (1 - kest);
1263 if (freqw < min_freqw || freqw > max_freqw)
1264 freqw = (max_freqw + min_freqw) / 2;
1271 meas_count += pin - pin0;
1273 while (meas_count >= meas_decimation)
1275 meas_count -= meas_decimation;
1277 freq_out->write(freq_tap);
1279 ss_out->write(sqrtf(est_insp));
1282 est_ep ? 10 * logf(est_sp / est_ep) / logf(10) : 0);
1292 freq_tap = freqw / 65536;
1320 template <
typename T>
1329 static const unsigned int chunk_size = 128;
1336 meas_decimation(1048576),
1337 pll_adjustment(1.0),
1340 out(_out, chunk_size),
1349 memset(hist, 0,
sizeof(hist));
1350 init_lookup_tables();
1356 min_omega = omega * (1 - tol);
1357 max_omega = omega * (1 + tol);
1358 update_freq_limits();
1363 freqw = freq * 65536;
1364 update_freq_limits();
1371 min_freqw = freqw - 65536 / max_omega / 8;
1372 max_freqw = freqw + 65536 / max_omega / 8;
1375 static const int RLUT_BITS = 8;
1376 static const int RLUT_ANGLES = 1 << RLUT_BITS;
1381 signed long freq_alpha = 0.04 * 65536;
1382 signed long freq_beta = 0.0012 * 256 * 65536 / omega * pll_adjustment;
1385 fail(
"Excessive oversampling");
1387 float gain_mu = 0.02 / (cstln_amp *
cstln_amp) * 2;
1389 int max_meas = chunk_size / meas_decimation + 1;
1393 while (
in.readable() >= chunk_size + 1 &&
1394 out.writable() >= chunk_size && (!freq_out || freq_out->writable() >= max_meas) && (!cstln_out || cstln_out->writable() >= max_meas))
1396 complex<T> *pin =
in.rd(), *pin0 = pin, *pend = pin + chunk_size;
1397 hardsymbol *pout =
out.wr(), *pout0 = pout;
1400 u_angle symbol_arg = 0;
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];
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 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);
1444 int quadrant = symbol_arg >> 14;
1445 static unsigned char quadrant_to_symbol[4] = {0, 2, 3, 1};
1446 *pout = quadrant_to_symbol[quadrant];
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;
1461 #define HIST_FLOAT 0 1463 hist[0].p.re = (float)s.
re - 128;
1464 hist[0].p.im = (
float)s.
im - 128;
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;
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);
1477 hist[0].c = arg_to_symbol((symbol_arg & 49152) + 8192);
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));
1482 float mucorr = muerr * gain_mu;
1483 const float max_mucorr = 0.1;
1485 if (mucorr < -max_mucorr)
1486 mucorr = -max_mucorr;
1487 if (mucorr > max_mucorr)
1488 mucorr = max_mucorr;
1499 in.read(pin - pin0);
1500 out.written(pout - pout0);
1502 if (symbol_arg && cstln_out)
1504 cstln_out->write(s);
1511 if (freqw < min_freqw || freqw > max_freqw)
1512 freqw = (max_freqw + min_freqw) / 2;
1517 meas_count += pin - pin0;
1519 while (meas_count >= meas_decimation)
1521 meas_count -= meas_decimation;
1524 freq_out->write((
float)freqw / 65536);
1535 } lut_polar[256][256];
1540 return lut_polar[c.
re][c.
im].a;
1543 cu8 lut_rect[RLUT_ANGLES][256];
1544 cu8 lut_sincos[65536];
1548 return lut_sincos[a];
1553 for (
int i = 0;
i < 256; ++
i)
1555 for (
int q = 0; q < 256; ++q)
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);
1563 for (
unsigned long a = 0; a < 65536; ++a)
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);
1570 for (
int a = 0; a < RLUT_ANGLES; ++a)
1572 for (
int r = 0; r < 256; ++r)
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));
1605 template <
typename Tout,
int Zout>
1622 fail(
"constellation not set");
1624 int count =
min(
in.readable(),
out.writable());
1625 u8 *pin =
in.rd(), *pend = pin + count;
1628 for (; pin < pend; ++pin, ++pout)
1631 pout->
re = Zout + cp->
re;
1632 pout->
im = Zout + cp->
im;
1649 template <
typename T>
1655 float freq) :
runnable(sch,
"rotator"),
1660 int ifreq = freq * 65536;
1662 fprintf(stderr,
"Rotate: req=%f real=%f\n", freq, ifreq / 65536.0);
1664 for (
int i = 0;
i < 65536; ++
i)
1666 lut_cos[
i] = cosf(2 *
M_PI *
i * ifreq / 65536);
1667 lut_sin[
i] = sinf(2 *
M_PI *
i * ifreq / 65536);
1673 unsigned long count =
min(
in.readable(),
out.writable());
1677 for (; pin < pend; ++pin, ++pout, ++index)
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;
1692 float lut_cos[65536];
1693 float lut_sin[65536];
1707 template <
typename T>
1713 float _bandwidth,
int nfft = 4096) :
runnable(sch,
"cnr_fft"),
1714 bandwidth(_bandwidth),
1725 if (bandwidth > 0.25)
1726 fail(
"CNR estimator requires Fsampling > 4x Fsignal");
1736 while (
in.readable() >=
fft.
n &&
out.writable() >= 1)
1740 if (
phase >= decimation)
1753 float center_freq = freq_tap ? *freq_tap * tap_multiplier : 0;
1754 int icf =
floor(center_freq *
fft.
n + 0.5);
1756 memcpy(data,
in.rd(),
fft.
n *
sizeof(data[0]));
1758 T *power =
new T[
fft.
n];
1760 for (
int i = 0;
i <
fft.
n; ++
i)
1761 power[
i] = data[
i].re * data[
i].re + data[
i].im * data[
i].im;
1766 avgpower =
new T[
fft.
n];
1767 memcpy(avgpower, power,
fft.
n *
sizeof(avgpower[0]));
1771 for (
int i = 0;
i <
fft.
n; ++
i)
1772 avgpower[
i] = avgpower[
i] * (1 - kavg) + power[
i] * kavg;
1774 int bwslots = (bandwidth / 4) *
fft.
n;
1784 float c2plusn2 = avgslots(icf - bwslots, icf + bwslots);
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;
1798 for (
int i = i0;
i <=
i1; ++
i)
1799 s += avgpower[
i & (
fft.
n - 1)];
1801 return s / (i1 - i0 + 1);
1812 template <
typename T,
int NFFT>
1834 while (
in.readable() >=
fft.
n * decim &&
out.writable() >= 1)
1838 if (
phase >= decimation)
1855 memcpy(data,
in.rd(),
fft.
n *
sizeof(data[0]));
1861 for (
int i = 0;
i <
fft.
n; ++
i, pin += decim)
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;
1874 avgpower =
new float[
fft.
n];
1875 memcpy(avgpower, power,
fft.
n *
sizeof(avgpower[0]));
1879 for (
int i = 0;
i <
fft.
n; ++
i)
1880 avgpower[
i] = avgpower[
i] * (1 - kavg) + power[
i] * kavg;
1883 for (
int i = 0;
i <
fft.
n / 2; ++
i)
1885 power[
i] = 10 * log10f(avgpower[NFFT / 2 +
i]);
1886 power[NFFT / 2 +
i] = 10 * log10f(avgpower[
i]);
1889 memcpy(
out.wr(), power,
sizeof(power[0]) * NFFT);
1903 #endif // LEANSDR_SDR_H
unsigned long window_size
simple_agc(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< complex< T >> &_out)
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)
spectrum(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< float[NFFT]> &_out)
cstln_lut< SOFTSYMB, 256 > * cstln
pipewriter< cf32 > * cstln_out
pipereader< complex< T > > in
ss_estimator(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< T > &_out)
Fixed< IntType, IntBits > cos(Fixed< IntType, IntBits > const &x)
void set_omega(float _omega, float tol=10e-6)
void to_softsymb(const full_ss *fss, hard_ss *ss)
complex< signed char > polar(float r, int n, float i)
void set_freq(float freq)
fast_qpsk_receiver(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< hardsymbol > &_out, pipebuf< float > *_freq_out=NULL, pipebuf< complex< T >> *_cstln_out=NULL)
fir_sampler(int _ncoeffs, Tc *_coeffs, int _subsampling=1)
pipereader< complex< T > > in
unsigned long meas_decimation
cstln_lut< hard_ss, 256 > * cstln
pipewriter< SOFTSYMB > out
Fixed< IntType, IntBits > floor(Fixed< IntType, IntBits > const &x)
pipewriter< complex< Tout > > out
pipewriter< float > * mer_out
pipewriter< cu8 > * cstln_out
pipewriter< float[NFFT]> out
cnr_fft(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< float > &_out, float _bandwidth, int nfft=4096)
unsigned long window_size
pipereader< complex< T > > in
virtual complex< T > interp(const complex< T > *pin, float mu, float phase)=0
void make_lut_from_symbols(float mer)
void set_allow_drift(bool d)
void softsymb_harden(llr_ss *ss)
pipewriter< float > * ss_out
unsigned long meas_decimation
uint8_t softsymb_to_dump(const llr_ss &ss, int bit)
rotator(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< complex< T >> &_out, float freq)
pipereader< complex< T > > in
complex< T > interp(const complex< T > *pin, float mu, float phase)
result * lookup(float I, float Q)
Fixed< IntType, IntBits > sin(Fixed< IntType, IntBits > const &x)
cstln_lut(cstln_base::predef type, float mer=10, float gamma1=0, float gamma2=0, float gamma3=0)
pipereader< complex< T > > in
complex< T > interp(const complex< T > *pin, float mu, float phase)
pipewriter< complex< T > > out
u_angle fast_arg(const cu8 &c)
virtual void update_freq(float freqw, int period=1)
pipereader< complex< T > > in
void written(unsigned long n)
void set_omega(float _omega, float tol=10e-6)
pipereader< complex< T > > in
void polar2(int i, float r, float a0, float a1, float a2, float a3)
auto_notch(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< complex< T >> &_out, int _nslots, T _agc_rms_setpoint)
virtual ~sampler_interface()
pipewriter< complex< T > > out
float avgslots(int i0, int i1)
sampler_interface< T > * sampler
void update_freq_limits()
pipewriter< complex< T > > out
void update_freq_limits()
void init_lookup_tables()
result * lookup(int I, int Q)
std::complex< Fixed< IntType, IntBits > > polar(const Fixed< IntType, IntBits > &rho, const Fixed< IntType, IntBits > &theta)
void update_freq(float freqw, int period)
pipewriter< hardsymbol > out
complex< int8_t > * symbols
complex< T > interp(const complex< T > *pin, float mu, float phase)
cu8 arg_to_symbol(u_angle a)
void set_freq(float freq)
virtual int readahead()=0
ss_amp_estimator(scheduler *sch, pipebuf< complex< T >> &_in, pipebuf< T > &_out_ss, pipebuf< T > &_out_ampmin, pipebuf< T > &_out_ampmax)
pipereader< complex< T > > in
void do_update_freq(float freqw)
void update_freq(float _freqw, int period=1)
cstln_transmitter(scheduler *sch, pipebuf< u8 > &_in, pipebuf< complex< Tout >> &_out)
void read(unsigned long n)
void inplace(complex< T > *data, bool reverse=false)
T min(const T &x, const T &y)