|
| 1 | +#pragma once |
| 2 | +// Active Constellation Extension (ACE) PAPR reduction. |
| 3 | +// |
| 4 | +// WHY ACE (not clipping, not tone reservation) — measured on the rig 2026-07-13: |
| 5 | +// - Crude amplitude clipping recovers ~4.5 dB of average power (real, |
| 6 | +// hardware-only benefit under peak-normalization) but injects EVM that |
| 7 | +// CRATERS clean channels (−58%: it damages a signal the RX decodes fine). |
| 8 | +// - Tone reservation needs ~4-6 spare in-band carriers; our spectrum is |
| 9 | +// packed (59 carriers in ~61 in-band bins, only DC+1 edge spare), so it |
| 10 | +// would require a frame-geometry redesign. |
| 11 | +// ACE fits a packed spectrum: it extends OUTER constellation points OUTWARD |
| 12 | +// (away from their decision boundaries) to cancel time-domain peaks. The RX |
| 13 | +// sees points FARTHER from where it would misread them, so BER is |
| 14 | +// equal-or-better — never worse. No reserved carriers, no frame change, no |
| 15 | +// RX change (the demapper's decision regions are unbounded on the outward |
| 16 | +// side). The two invariants below are what make it safe; the unit test |
| 17 | +// (test_papr_ace) proves both hold. |
| 18 | +// |
| 19 | +// INVARIANT 1 (safety): every data carrier moves ONLY in its outward-safe |
| 20 | +// region — radially outward for PSK, axis-outward on an outer level for QAM. |
| 21 | +// No point moves toward any decision boundary. Pilots and empty bins are |
| 22 | +// never touched (channel estimation and the reserved bins are preserved). |
| 23 | +// INVARIANT 2 (benefit): post-PAPR <= pre-PAPR (peaks only shrink). |
| 24 | +// |
| 25 | +// Algorithm: clip-and-project (Krongold & Jones 2003). Clip the time signal, |
| 26 | +// take the clip residual's spectrum, add back ONLY the outward-safe component |
| 27 | +// on data carriers, iterate a few times. The projection is what confines the |
| 28 | +// correction to the safe region. |
| 29 | + |
| 30 | +#include "ultra/types.hpp" |
| 31 | +#include "ultra/dsp.hpp" |
| 32 | +#include "ofdm/demodulator_constants.hpp" |
| 33 | + |
| 34 | +#include <algorithm> |
| 35 | +#include <cmath> |
| 36 | +#include <vector> |
| 37 | + |
| 38 | +namespace ultra { |
| 39 | +namespace ofdm { |
| 40 | +namespace papr_ace { |
| 41 | + |
| 42 | +struct AceResult { |
| 43 | + bool applied = false; |
| 44 | + float pre_papr_db = 0.0f; |
| 45 | + float post_papr_db = 0.0f; |
| 46 | + int iterations = 0; |
| 47 | + size_t data_carriers_extended = 0; |
| 48 | +}; |
| 49 | + |
| 50 | +// Outer-level magnitude (in scaled constellation units) per QAM axis: the |
| 51 | +// point is "outer" on an axis if its coordinate sits beyond the last inner |
| 52 | +// decision boundary, so extending it further out crosses no boundary. |
| 53 | +// 16QAM I/Q levels {-3,-1,1,3}·s → outer at ±3·s, inner/outer boundary at 2·s. |
| 54 | +// 32QAM (cross): I {-3..3}·s, Q {-7..7}·s → outer boundary 2·s (I) / 6·s (Q). |
| 55 | +struct QamAxisSpec { |
| 56 | + float scale = 0.0f; // constellation scale (1/sqrt(Es)) |
| 57 | + float i_outer_boundary = 0.0f; // |I| beyond this ⇒ outer on I |
| 58 | + float q_outer_boundary = 0.0f; // |Q| beyond this ⇒ outer on Q |
| 59 | + bool valid = false; |
| 60 | +}; |
| 61 | + |
| 62 | +// QAM16_SCALE lives in modulator.cpp (not the shared demod constants); mirror |
| 63 | +// the value here (1/sqrt(10)) so ACE and the TX constellation agree exactly. |
| 64 | +inline constexpr float kAceQam16Scale = 0.3162277660168379f; // 1/sqrt(10) |
| 65 | + |
| 66 | +inline QamAxisSpec qamAxisSpec(Modulation mod) { |
| 67 | + using namespace demod_constants; |
| 68 | + QamAxisSpec s; |
| 69 | + if (mod == Modulation::QAM16) { |
| 70 | + s.scale = kAceQam16Scale; |
| 71 | + s.i_outer_boundary = 2.0f * kAceQam16Scale; // between level 1 and 3 |
| 72 | + s.q_outer_boundary = 2.0f * kAceQam16Scale; |
| 73 | + s.valid = true; |
| 74 | + } else if (mod == Modulation::QAM32) { |
| 75 | + s.scale = QAM32_SCALE; |
| 76 | + s.i_outer_boundary = 2.0f * QAM32_SCALE; // I levels ±1,±3 → outer ±3 |
| 77 | + s.q_outer_boundary = 6.0f * QAM32_SCALE; // Q levels ±1..±7 → outer ±7 |
| 78 | + s.valid = true; |
| 79 | + } |
| 80 | + return s; |
| 81 | +} |
| 82 | + |
| 83 | +inline bool isPskMod(Modulation mod) { |
| 84 | + // All constellations whose decisions are by phase sector (points on a |
| 85 | + // circle) → safe extension is radially outward. QAM8 is the 8PSK |
| 86 | + // constellation (types.hpp:69), so it belongs here, not with rectangular QAM. |
| 87 | + return mod == Modulation::DBPSK || mod == Modulation::BPSK || |
| 88 | + mod == Modulation::DQPSK || mod == Modulation::QPSK || |
| 89 | + mod == Modulation::D8PSK || mod == Modulation::QAM8; |
| 90 | +} |
| 91 | + |
| 92 | +// Project a desired correction `delta` for a data symbol `s` onto its |
| 93 | +// outward-safe region. Returns the allowed part (added to s stays decodable). |
| 94 | +inline Complex aceProject(Complex s, Complex delta, Modulation mod, |
| 95 | + const QamAxisSpec& qam) { |
| 96 | + if (isPskMod(mod)) { |
| 97 | + // PSK decisions are by phase sector: safe = radially OUTWARD only |
| 98 | + // (increase magnitude, never rotate). Project delta onto the unit |
| 99 | + // vector s/|s| and keep only the outward (positive-radial) part. |
| 100 | + const float mag = std::abs(s); |
| 101 | + if (mag < 1e-9f) return Complex(0.0f, 0.0f); |
| 102 | + const Complex u = s / mag; |
| 103 | + const float radial = delta.real() * u.real() + delta.imag() * u.imag(); |
| 104 | + if (radial <= 0.0f) return Complex(0.0f, 0.0f); |
| 105 | + return u * radial; |
| 106 | + } |
| 107 | + if (qam.valid) { |
| 108 | + // QAM decisions are per-axis at integer half-levels: safe = extend an |
| 109 | + // OUTER-level coordinate further out on that axis only. Inner |
| 110 | + // coordinates are pinned (moving them crosses a boundary). |
| 111 | + float dI = 0.0f, dQ = 0.0f; |
| 112 | + if (s.real() > qam.i_outer_boundary && delta.real() > 0.0f) dI = delta.real(); |
| 113 | + else if (s.real() < -qam.i_outer_boundary && delta.real() < 0.0f) dI = delta.real(); |
| 114 | + if (s.imag() > qam.q_outer_boundary && delta.imag() > 0.0f) dQ = delta.imag(); |
| 115 | + else if (s.imag() < -qam.q_outer_boundary && delta.imag() < 0.0f) dQ = delta.imag(); |
| 116 | + return Complex(dI, dQ); |
| 117 | + } |
| 118 | + return Complex(0.0f, 0.0f); // unknown mod → no extension (safe no-op) |
| 119 | +} |
| 120 | + |
| 121 | +inline float paprDb(const std::vector<Complex>& x) { |
| 122 | + double sum_sq = 0.0, peak_sq = 0.0; |
| 123 | + for (const auto& v : x) { |
| 124 | + const double p = static_cast<double>(v.real()) * v.real() + |
| 125 | + static_cast<double>(v.imag()) * v.imag(); |
| 126 | + sum_sq += p; |
| 127 | + if (p > peak_sq) peak_sq = p; |
| 128 | + } |
| 129 | + if (sum_sq <= 0.0) return 0.0f; |
| 130 | + const double mean = sum_sq / static_cast<double>(x.size()); |
| 131 | + return static_cast<float>(10.0 * std::log10(peak_sq / mean)); |
| 132 | +} |
| 133 | + |
| 134 | +// In-place ACE on the freq-domain symbol. `clip_ratio` is the peak/RMS target |
| 135 | +// in linear (e.g. 10^(target_db/20)); `mu` is the per-iteration step (0..1); |
| 136 | +// `max_iters` typically 3-4. Data carriers only; pilots/empty bins untouched. |
| 137 | +inline AceResult applyAce(std::vector<Complex>& freq_domain, |
| 138 | + const std::vector<int>& data_indices, |
| 139 | + Modulation mod, FFT& fft, |
| 140 | + float clip_ratio, int max_iters, float mu) { |
| 141 | + AceResult r; |
| 142 | + const size_t N = freq_domain.size(); |
| 143 | + if (N == 0 || data_indices.empty() || clip_ratio <= 0.0f) return r; |
| 144 | + const QamAxisSpec qam = qamAxisSpec(mod); |
| 145 | + if (!isPskMod(mod) && !qam.valid) return r; // unsupported mod |
| 146 | + |
| 147 | + std::vector<Complex> x(N), excess(N), E(N); |
| 148 | + fft.inverse(freq_domain.data(), x.data()); |
| 149 | + r.pre_papr_db = paprDb(x); |
| 150 | + |
| 151 | + // FIXED clip level from the INITIAL RMS: ACE extends points outward, which |
| 152 | + // raises RMS — recomputing A each iteration would chase the rising RMS and |
| 153 | + // stop clipping (self-defeating; measured 0.1 dB). Holding A at the target |
| 154 | + // peak level lets successive iterations keep pushing the peak down toward it. |
| 155 | + double init_sq = 0.0; |
| 156 | + for (const auto& v : x) init_sq += static_cast<double>(v.real()) * v.real() + |
| 157 | + static_cast<double>(v.imag()) * v.imag(); |
| 158 | + const float init_rms = std::sqrt(static_cast<float>(init_sq / static_cast<double>(N))); |
| 159 | + const float A = clip_ratio * init_rms; |
| 160 | + if (A <= 0.0f) return r; |
| 161 | + |
| 162 | + for (int it = 0; it < max_iters; ++it) { |
| 163 | + bool any = false; |
| 164 | + for (size_t i = 0; i < N; ++i) { |
| 165 | + const float m = std::abs(x[i]); |
| 166 | + if (m > A) { |
| 167 | + excess[i] = x[i] * ((m - A) / m); // the part above the clip level |
| 168 | + any = true; |
| 169 | + } else { |
| 170 | + excess[i] = Complex(0.0f, 0.0f); |
| 171 | + } |
| 172 | + } |
| 173 | + if (!any) break; |
| 174 | + |
| 175 | + // X-domain rep of the time residual: forward(excess). This FFT |
| 176 | + // convention has forward∘inverse = identity (the modulator uses |
| 177 | + // inverse for X→time), so no 1/N — E[bin] is already in freq_domain's |
| 178 | + // units, ready to add to a data carrier. |
| 179 | + fft.forward(excess.data(), E.data()); |
| 180 | + size_t extended = 0; |
| 181 | + for (int bin : data_indices) { |
| 182 | + const Complex s = freq_domain[bin]; |
| 183 | + // Subtract the residual on this carrier (reduces the peak), but only |
| 184 | + // the outward-safe component. |
| 185 | + const Complex delta_desired = E[bin] * (-mu); |
| 186 | + const Complex allowed = aceProject(s, delta_desired, mod, qam); |
| 187 | + if (allowed.real() != 0.0f || allowed.imag() != 0.0f) { |
| 188 | + freq_domain[bin] = s + allowed; |
| 189 | + ++extended; |
| 190 | + } |
| 191 | + } |
| 192 | + r.data_carriers_extended = std::max(r.data_carriers_extended, extended); |
| 193 | + fft.inverse(freq_domain.data(), x.data()); |
| 194 | + r.iterations = it + 1; |
| 195 | + if (extended == 0) break; // nothing safe left to do |
| 196 | + } |
| 197 | + |
| 198 | + r.post_papr_db = paprDb(x); |
| 199 | + r.applied = r.iterations > 0; |
| 200 | + return r; |
| 201 | +} |
| 202 | + |
| 203 | +} // namespace papr_ace |
| 204 | +} // namespace ofdm |
| 205 | +} // namespace ultra |
0 commit comments