diff --git a/.claude/skills/comnumpy-tutorial/SKILL.md b/.claude/skills/comnumpy-tutorial/SKILL.md index e9b6121..152f6bc 100644 --- a/.claude/skills/comnumpy-tutorial/SKILL.md +++ b/.claude/skills/comnumpy-tutorial/SKILL.md @@ -59,9 +59,30 @@ then learns the wrong lesson -- they came to see how the library is used. of **one** chain are not a Python loop: grow the stimulus a leading batch axis -- `(n_trials, N)` -- and every block broadcasts over it, stochastic blocks draw independently per trial, and the pooled metric - is exactly the mean of the per-trial rates (D51; the exception is a - channel whose realization is frozen at construction -- a new draw per - trial is a new instance, and that loop stays visible). When the study + is exactly the mean of the per-trial rates (D51). + **Averaging over a random channel (fading, MIMO)** is the same idea + with the draws stacked into the configuration: draw a **fixed** set + of channels once -- `rayleigh_channel(size=K)`, one draw = one batch + row, the SAME K draws at every sweep point, never a draw count that + grows with the SNR (the reader compares curves, not draw counts) -- + hand the stack to the channel block AND to its detector at + construction, and sweep the noise variance with `monte_carlo`: no + simulation loop at all. A detector that weights by the noise + (`MMSE`, `OSIC`) receives the swept variance too, by zipping two + dotted paths: `monte_carlo(chain, ("noise.sigma2", + "detector.sigma2"), list(zip(values, values)), ...)`. + **On the page, interleave**: one technique is one unit -- its idea + and equations, then its `Sequential` and its `monte_carlo` call + quoted together as one `literalinclude` block. The reader meets each + method whole: theory, chain, sweep. + **Never a loop over the techniques.** A `for name, chain in (...)` + around two or three `monte_carlo` calls is over-engineering: write + the calls out, one explicit statement per chain, filling + `curves[name]` by name. A loop earns its place at ten iterations, + not at two. The same goes for storage: two curves are a dict literal + of two `np.zeros`, not an init loop; and a runtime sweep is ONE loop + over the swept variable with each chain timed explicitly in the + body, never a nested loop over the techniques. When the study is *around* a chain -- several detectors on one frame, several receivers on one propagation -- write the loop out, on **one storage convention**: declare the methods and the metrics first; pre-allocate @@ -115,34 +136,23 @@ then learns the wrong lesson -- they came to see how the library is used. ### Naming the chain -Three rules, and they are about where the reader looks and what the -signature promises. - -1. **A function only when the chain is built more than once** -- several - values of a parameter, several variants, a Monte-Carlo. A chain built - once and used once is written inline: wrapping four blocks in a - function called once, with no argument, is ceremony, and it puts an - indirection between the reader and their first chain. -2. **Always `get_()`**, one prefix, the noun saying what comes - back. Which noun follows from how the page is organized: - - | the page… | names | - |---|---| - | has one chain | `get_chain(...)` | - | cuts one chain in two | `get_channel()` and `get_receiver()` | - | compares two structurally different chains | `get_uncoded_chain()`, `get_coded_chain(soft)` | - | builds a sub-assembly, not a whole chain | the part: `get_transmitter()` | - - Two chains that differ by a block or two are **one** function with a - parameter, not two functions. Two chains that differ in structure -- - four blocks against six -- are two functions, because merging them - would hide behind a flag exactly what the page is comparing. -3. **What varies is a parameter, never a module global.** Constants of - the page (the fibre, the constellation, the roll-off) may stay - global; anything the page *varies* goes through the signature. A - function whose signature announces three arguments and silently reads - eleven more is a closure wearing a `def`, and the reader cannot tell - from the call site what changed. +**Default: no function. One `Sequential` per technique, written out in +full, before the sweep.** A page that compares SISO, Alamouti and MRC +writes three `Sequential([...])` blocks -- even if they share most of +their blocks -- because the reader must see each chain whole; a factory +with a `kind` flag hides behind the flag exactly what the page is +comparing, and a factory called once with no argument is ceremony. The +batch contract (D51) is what makes this possible: the channel draws are +stacked into the chain at construction (see the Monte-Carlo bullet), so +nothing needs rebuilding inside a loop -- there is no loop. + +A `def` remains legitimate in exactly one case: the chain is genuinely +rebuilt with different *structure* along the page (a DBP receiver with +a varying step count, a coded chain against an uncoded one built in +several places). Then it is `get_()`, the noun saying what comes +back (`get_chain`, `get_receiver`, `get_transmitter`), and what varies +is a parameter of the signature, never a module global read silently +from inside. ### The rest of the script @@ -171,6 +181,10 @@ signature promises. one and 9.1 ms in the other. - State the estimator's floor when a curve hits one (an SER of 1e-4 from 8192 symbols is the estimator, not the link). +- **No analysis codas.** No `np.interp` readings of a closed form, no + diversity-order polyfits, no measured/theory ratio tables, no mermaid + dumps at the end of a script: they drown the reader. What a figure + shows is said in one sentence of prose next to the figure. ## Register diff --git a/CHANGELOG.md b/CHANGELOG.md index 4cdfa83..03bd35d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,78 @@ one release; there is no compatibility layer. | `core.metrics.calculate_acpr` | `compute_acpr` — it was the only `calculate_*` in the library, against 17 `compute_*` | | `core.metrics.compute_effective_SNR`, `ofdm.metrics.compute_PAPR` | `compute_effective_snr`, `compute_papr` — the two capitalized outliers among functions otherwise all lowercase (`compute_ser`, `compute_ber`, `compute_evm`, `compute_ccdf`, `compute_mi`) | +### Changed — every tutorial script held to the writing doctrine + +The conventions sweep closed over the whole tutorial set. The coding +page: three chains written out in place of a factory pair, the LDPC +iteration count varied by `set_params(decoder__n_iter=)` on one chain +instead of rebuilding it in a loop, two `print_data` tables in place of +hand-rolled aligned prints. The AWGN page: the chain inline, the +constellation-order study reconfiguring it with `set_params`, the +diagram section gone (its feature demo lives in the getting-started +page). The OFDM page: two explicit `monte_carlo` calls and a runtime +loop over the swept lengths only. The MIMO cost study and the PAPR +references: unrolled. The shaping and GN-model scripts lose their +mermaid dumps and stray `plt.show`; the GN-model factory stays -- its +source genuinely changes structure (a Gaussian stimulus against a +mapped one) -- and so do the two DBP scripts' factories, the one +sanctioned chain split, for the Monte-Carlo cost of re-propagating the +fibre. The interpolation the shaping page keeps reads its own computed +capacity curves, not a closed form. + +### Changed — the Alamouti and MIMO tutorials rewritten on the batch + +Both pages carried simulation scaffolding that drowned the reader: a +factory building five chains behind a `kind` flag, per-SNR draw counts, +zipped channel sweeps, interpolation codas reading numbers off closed +forms, mermaid dumps. They are rewritten on one doctrine, now in the +tutorial skill. The page opens on **one** chain, one channel draw, no +batch -- the problem shown plainly. The simulation then draws a +**fixed** stack of channels once (`rayleigh_channel(size=K)`, the same +K draws at every SNR point: the reader compares curves, not draw +counts), builds **one `Sequential` per technique** -- written out in +full, no factory -- with the stack on the channel block and on the +detector, and runs **one `monte_carlo` per technique** sweeping the +noise variance, zipped into the detector when it weights by it (MMSE, +OSIC): no simulation loop at all. On the MIMO page each technique is +one unit -- equations, chain, sweep -- interleaved. To make it +possible, `SpaceTimeDecoder` accepts a stacked H (channel k decodes +frame k, locked in `tests/mimo/test_stacked_channel.py`) and both +Alamouti blocks join the BROADCAST register of the batch contract. + +### Changed — the 91 docstrings audited against their code + +Five parallel reviewers read every `Processor` docstring against its +`forward()` -- equations, Axes lines, parameter lists, doctests -- and +62 findings were fixed. Four were **code** bugs the prose exposed: +`Resampler` resampled axis 0 (scipy's default) while promising the last +axis; `SelectiveMIMOChannel` validated its input against `N_r` instead +of `N_t`, rejecting every non-square channel's correct frames; +`BlindCFOCompensator.history` stored 16 times the offset (renamed +`history_`, D23, and fixed to store the omega_0 iterates it documents); +and `Laser` drew its "random" initial phase from a Gaussian of std 2 pi +instead of uniformly on [-pi, pi). + +The rest is the docstrings catching up with the code. Wrong equations +straightened: `BWFilter`'s mask (off by 2 vs its own parameters), +`AmplitudeDemapper`'s decision (on |Re y|, per quadrature lane, not +|y|), the Rice path of `TappedDelayLineChannel` (the gamma_0 scaling +the code applies), the ambiguous alpha_dB conversion in the optical +helpers. Stale claims deleted: seven phantom pre-D46 parameters on +`FiberLink` and `DBP` (now one `fiber: FiberSpec` entry), the unusable +`norm` field of `FrequencyDomainEqualizer` (removed -- it never reached +the FFT), 1-D-only Axes lines on blocks that broadcast (CD, Kerr, +EDFA, BlindPhaseTracker), "multiplied" for a division, `"time"` as a +method that never existed. The MIMO detectors' stored decisions are +`S_` (D23), the CMA equalizer is written W -- H is the channel +everywhere else -- and `BlindPhaseSearchCompensator` gained the +Attributes section its `phase_` deserved. The batch registers were +corrected on the way: `DelayRemover` and `BlindPhaseTracker` are +verified BROADCAST, `Normalizer` moved to EXEMPT -- its gain is one +scalar over the whole array by documented design, and its docstring now +warns that a batch is pooled. `PtsPaprReductor` validates its (T, F) +layout in `prepare()` instead of dying on a bare IndexError. + ### Added — the batch promise is a ratchet over the catalogue A convention over 91 blocks is only worth the sweep that checks it. diff --git a/CONVENTIONS.md b/CONVENTIONS.md index fb1ef1f..b40c502 100644 --- a/CONVENTIONS.md +++ b/CONVENTIONS.md @@ -41,7 +41,7 @@ Every block belongs to exactly one category, stated in its docstring ("Axes:" line): 1. **Element-wise** — operates pointwise, shape-agnostic - (`AWGN`, `SymbolMapper`, `HardClipper`). + (`AWGN`, `SymbolMapper`, `Clipper`). 2. **Axis -1** — hardcodes `axis=-1` and broadcasts over everything else (filters, `FFTProcessor`, `IFFTProcessor`, `CyclicPrefixer`). 3. **Declared axis** — requires a specific axis and validates it in diff --git a/docs/tutorials/alamouti.rst b/docs/tutorials/alamouti.rst index db028b9..bdc6cc2 100644 --- a/docs/tutorials/alamouti.rst +++ b/docs/tutorials/alamouti.rst @@ -53,14 +53,14 @@ We start with the imports and the parameters, QPSK over a Rayleigh channel: .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py :language: python - :lines: 1-31 + :lines: 1-29 Let us check the law on draws from :func:`~comnumpy.mimo.utils.rayleigh_channel`: .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py :language: python - :lines: 37-55 + :lines: 31-53 .. image:: img/one_shot_alamouti_fig1.png :width: 100% @@ -175,45 +175,29 @@ diversity costs no bandwidth**, and no feedback. Implementation ^^^^^^^^^^^^^^ -The three schemes differ by what is put on the antennas and what is done with -what comes back, so one function builds them all: +One channel draw, one chain, transmitter to decision. The scaling by +:math:`1/\sqrt{N_t}` matters for the comparison that follows: two antennas +each transmitting :math:`|s|^2` would spend twice the power of a single +antenna, a 3 dB advantage that has nothing to do with coding. Splitting the +power keeps the comparison about diversity alone. .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py :language: python - :lines: 60-91 - -The scaling by :math:`1/\sqrt{N_t}` matters for the comparison that follows. -Two antennas each transmitting :math:`|s|^2` would spend twice the power of a -single antenna, which is a 3 dB advantage that has nothing to do with coding. -Splitting the power keeps the comparison about diversity alone. - -.. mermaid:: mermaid/alamouti.mmd - -The diagram is not drawn by hand. It is what the chain says about itself -- -``chain.to_mermaid()`` (decision D33c) -- exported by the script, so the -block names are the ones the code uses and a dashed outline marks a tapped -block. - -One-shot simulation -""""""""""""""""""" - -.. literalinclude:: ../../examples/mimo/one_shot_alamouti.py - :language: python - :lines: 94-107 + :lines: 55-84 .. code:: one-shot SER: 0.0600 # block id output shape dtype time ms - 0 SymbolGenerator tx (1000,) int64 0.04 + 0 SymbolGenerator tx (1000,) int64 0.03 1 SymbolMapper symbol_mapper (1000,) complex128 0.00 2 Amplifier signal_amplifier (1000,) complex128 0.01 - 3 SpaceTimeEncoder space_time_encoder (2, 1000) complex128 0.10 + 3 SpaceTimeEncoder space_time_encoder (2, 1000) complex128 0.06 4 FlatMIMOChannel channel (1, 1000) complex128 0.01 - 5 AWGN noise (1, 1000) complex128 0.06 - 6 SpaceTimeDecoder detector (1000,) complex128 0.12 - 7 Amplifier signal_amplifier_2 (1000,) complex128 0.01 - 8 SymbolDemapper symbol_demapper (1000,) int64 0.08 + 5 AWGN noise (1, 1000) complex128 0.05 + 6 SpaceTimeDecoder detector (1000,) complex128 0.08 + 7 Amplifier signal_amplifier_2 (1000,) complex128 0.00 + 8 SymbolDemapper symbol_demapper (1000,) int64 0.06 The shape column is the code at work: 1000 symbols enter, the encoder spreads them over ``(2, 1000)`` -- two antennas, one thousand channel uses, so rate 1 @@ -233,63 +217,67 @@ the equations above. Monte Carlo Evaluation ^^^^^^^^^^^^^^^^^^^^^^ -Averaging over fading means running the chain once per channel realization, -which is a sweep whose parameter is the channel. :func:`~comnumpy.monte_carlo.sweep` -takes several dotted parameter names at once and zips them, so one sweep point -sets the channel the signal goes through **and** the channel the detector -inverts: +Averaging over fading takes many channel draws, and the draws are a +**batch**: the *same* 5000 realizations serve every SNR point, one draw per +row. Each scheme is one chain, built once with its stack of channels -- the +channel block propagates draw :math:`k` on frame :math:`k`, and the detector +holds the same stack: .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py :language: python - :lines: 112-169 - -The three curves are three readings of one closed form, -:func:`~comnumpy.core.metrics.compute_ser_rayleigh_psk`: :math:`L` branches -evaluated at the per-branch SNR, with a transmit scheme dividing that SNR by -:math:`N_t` because it splits its power over the antennas. ``plot_error_rate`` -draws measurements as hollow markers and their references as lines of the same -colour, so a pair reads as one statement: + :lines: 86-126 -.. image:: img/one_shot_alamouti_fig3.png - :width: 100% - :align: center +From there the sweep needs no simulation loop at all: +:func:`~comnumpy.monte_carlo.monte_carlo` moves the noise variance, and +everything else -- the channels, the chains -- is frozen: -The two statements this tutorial is about -- the slope and the 3 dB -- are -exact, so they are read off the closed form rather than fitted to the points; -the simulation is what confronts them: +.. literalinclude:: ../../examples/mimo/one_shot_alamouti.py + :language: python + :lines: 128-143 .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py :language: python - :lines: 171-210 + :lines: 145-164 .. code:: - 1 Tx, 1 Rx (no diversity) measured / closed form 0.99 0.97 0.94 0.90 0.83 0.86 - Alamouti, 2 Tx, 1 Rx measured / closed form 0.98 0.92 0.88 0.77 0.84 0.65 - MRC, 1 Tx, 2 Rx measured / closed form 0.98 0.91 0.82 0.72 0.73 0.28 - 1 Tx, 1 Rx (no diversity) diversity order 1.00 - Alamouti, 2 Tx, 1 Rx diversity order 2.00 - MRC, 1 Tx, 2 Rx diversity order 2.00 - SNR for SER = 0.001: MRC 15.6 dB, Alamouti 18.6 dB, gap 3.01 dB (10log10(N_t) = 3.01 dB) - -**The slope.** The single-antenna scheme has diversity order 1, the two -others 2: the error rate falls twice as fast, which is the whole reason a -space-time code exists. - -**The 3 dB.** Reaching a symbol error rate of :math:`10^{-3}` costs 15.6 dB -with two receive antennas and 18.6 dB with Alamouti, a gap of **3.01 dB** -against :math:`10\log_{10} N_t = 3.01` dB. That is the price of transmitting -*blind*: the receiver knows the channel and weights its branches by -:math:`h_i^{*}`, the transmitter cannot and splits its power evenly. Alamouti -buys the full diversity order anyway; it pays only in array gain, not in -slope. - -**The ratios.** The measurement tracks the closed form to a few percent at low -SNR and drifts below it as the SNR grows, down to 0.28 for the steepest curve: -6000 channel draws no longer sample the deep fades that dominate the average -there. Nothing is wrong with either side -- -``validation/mimo_diversity_ber.py`` runs the same three schemes with up to -80 000 draws and lands within 4.4 % of the same curves. + SER + snr_dB 1 Tx, 1 Rx (no diversity) Alamouti, 2 Tx, 1 Rx MRC, 1 Tx, 2 Rx + ------------------------------------------------------------------------ + 4 0.22261 1.778e-01 8.162e-02 + 8 0.11186 6.430e-02 2.155e-02 + 12 0.04976 1.675e-02 4.200e-03 + 16 0.02082 3.373e-03 5.700e-04 + 20 0.00843 5.525e-04 8.000e-05 + 24 0.00312 6.250e-05 2.500e-06 + +.. image:: img/one_shot_alamouti_fig3.png + :width: 100% + :align: center + +The three curves are three readings of one closed form, +:func:`~comnumpy.core.metrics.compute_ser_rayleigh_psk`: :math:`L` branches +evaluated at the per-branch SNR, with a transmit scheme dividing that SNR by +:math:`N_t` because it splits its power over the antennas. +``plot_error_rate`` draws measurements as hollow markers and their references +as lines of the same colour, so a pair reads as one statement -- and the two +statements this tutorial is about are on the figure: + +**The slope.** The single-antenna curve loses one decade of error rate per +10 dB; the two others fall twice as fast. That is diversity order 2, the +whole reason a space-time code exists. + +**The 3 dB.** The Alamouti and MRC curves are parallel, separated +horizontally by :math:`10\log_{10} N_t = 3` dB. That is the price of +transmitting *blind*: the receiver knows the channel and weights its +branches by :math:`h_i^{*}`; the transmitter cannot, and splits its power +evenly. Alamouti buys the full diversity order anyway -- it pays only in +array gain, never in slope. + +The last MRC point sits on the estimator's floor: 2.5e-6 out of 400 000 +symbols is one error, not a rate. ``validation/mimo_diversity_ber.py`` runs +the same three schemes with up to 80 000 draws and lands within 4.4 % of the +closed forms. Beyond two antennas diff --git a/docs/tutorials/awgn.rst b/docs/tutorials/awgn.rst index 6d416c2..0ffebb5 100644 --- a/docs/tutorials/awgn.rst +++ b/docs/tutorials/awgn.rst @@ -81,7 +81,7 @@ this page builds one per constellation order. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 25-42 + :lines: 25-33 The processors are: @@ -99,19 +99,6 @@ The processors are: - ``SymbolDemapper`` Maps received noisy constellation points back to integers. -The chain, as the chain itself describes it: - -.. mermaid:: mermaid/awgn_chain.mmd - -The diagram above is not drawn by hand. It is what the chain says about -itself -- ``chain.to_mermaid()`` (decision D33c) -- exported by the -script, so the block names are the ones the code uses and a dashed -outline marks a tapped block: - -.. literalinclude:: ../../examples/simple/monte_carlo_awgn.py - :language: python - :lines: 129-135 - Monte Carlo Simulation """""""""""""""""""""" @@ -122,7 +109,7 @@ against the transmitted symbols. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 44-53 + :lines: 35-44 Three chain services appear there, and they are the ones every study is made of. @@ -141,7 +128,7 @@ at every point: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 55-66 + :lines: 46-58 .. code:: @@ -209,7 +196,7 @@ two expressions above transcribe line for line: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 69-82 + :lines: 60-73 From the constellation """""""""""""""""""""" @@ -222,7 +209,7 @@ by :math:`k` above happens inside. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 84-89 + :lines: 75-81 .. code:: @@ -244,7 +231,7 @@ share a colour, so a pair reads as one statement. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 91-97 + :lines: 83-89 .. image:: img/monte_carlo_awgn.png :width: 100% @@ -272,7 +259,7 @@ not guess (decision D41), so it happens at the call site: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 99-120 + :lines: 91-114 The metric changes with the axis. Against :math:`E_b/N_0` the natural quantity is the bit error rate, and ``compute_ber`` needs the symbol width, which @@ -291,7 +278,7 @@ disagreeing with theory. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 120-127 + :lines: 116-123 .. image:: img/monte_carlo_awgn_orders.png :width: 100% diff --git a/docs/tutorials/coding.rst b/docs/tutorials/coding.rst index 23424fb..83931dc 100644 --- a/docs/tutorials/coding.rst +++ b/docs/tutorials/coding.rst @@ -61,7 +61,7 @@ simple -- and the codes below all have rate 1/2. .. literalinclude:: ../../examples/simple/channel_coding.py :language: python - :lines: 22-31 + :lines: 21-29 That helper deserves a line of its own, because it is where most coded simulations go wrong. The channel sees *symbols*, so it is parameterized by @@ -106,23 +106,29 @@ a path in a **trellis**. That trellis is the object the decoder searches. .. literalinclude:: ../../examples/simple/channel_coding.py :language: python - :lines: 33-55 + :lines: 32-46 .. code:: generators ('0o133', '0o171') K = 7 states = 64 rate = 0.5 4 bits in -> [1 1 0 1 0 0 0 1 1 0 1 0 0 0 1 0 0 1 1 1] (with the tail) free distance d_free = 10 - d 10 11 12 13 14 15 - a_d 11 0 38 0 193 0 - beta_d 36 0 211 0 1404 0 + distance spectrum + d a_d beta_d + --------------- + 10 11 36 + 11 0 0 + 12 38 211 + 13 0 0 + 14 193 1404 + 15 0 0 Four input bits produced twenty output bits rather than eight: the encoder is **terminated**, so :math:`K - 1 = 6` zero bits are appended to flush the register and bring it back to the all-zero state. On a short block that overhead is visible; on a realistic one it is not. -The last three lines are the code's **distance spectrum**, and they are worth +The table is the code's **distance spectrum**, and they are worth more than they look. :math:`d_{\mathrm{free}} = 10` is the smallest Hamming distance between two distinct code sequences; :math:`a_d` counts the error events at distance :math:`d`, and :math:`\beta_d` counts the *information* @@ -168,17 +174,21 @@ wrong side. The information thrown away by deciding too early is worth about .. literalinclude:: ../../examples/simple/channel_coding.py :language: python - :lines: 57-89 - -The chain, as the chain itself describes it: - -.. mermaid:: mermaid/channel_coding.mmd + :lines: 49-90 .. code:: - uncoded 7.86e-02 5.68e-02 3.79e-02 2.28e-02 1.31e-02 6.45e-03 2.80e-03 7.75e-04 (0.0 s) - hard-decision Viterbi 3.71e-01 2.43e-01 1.14e-01 3.15e-02 7.52e-03 9.75e-04 0.00e+00 0.00e+00 (8.1 s) - soft-decision Viterbi 1.40e-01 4.32e-02 3.63e-03 3.50e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 (8.9 s) + BER + Eb/N0 [dB] uncoded hard-decision Viterbi soft-decision Viterbi + ------------------------------------------------------------------- + 0 7.865e-02 3.708e-01 1.401e-01 + 1 5.680e-02 2.434e-01 4.320e-02 + 2 3.787e-02 1.143e-01 3.625e-03 + 3 2.278e-02 3.150e-02 3.500e-04 + 4 1.305e-02 7.525e-03 0.000e+00 + 5 6.450e-03 9.750e-04 0.000e+00 + 6 2.800e-03 0.000e+00 0.000e+00 + 7 7.750e-04 0.000e+00 0.000e+00 Three readings, and the first one is a warning. @@ -195,7 +205,7 @@ needs 5 dB and soft decoding a little under 3 -- the classic figure for this code. Nothing was added to the receiver except *not throwing information away*: the same trellis, the same recursion, a different branch metric. -**And the last two columns are not results.** Zero errors in 40 000 bits does +**And the zero entries are not results.** Zero errors in 40 000 bits does not mean the BER is zero; it means it is somewhere below roughly :math:`1/40000 = 2.5 \times 10^{-5}` and this simulation cannot see it. The figure drops those points rather than drawing them at zero on a logarithmic @@ -224,11 +234,21 @@ closed form: no simulation, no random draw, and no floor at .. literalinclude:: ../../examples/simple/channel_coding.py :language: python - :lines: 91-101 + :lines: 92-106 .. code:: - union bound 1.82e+01 9.32e-01 2.85e-02 6.81e-04 1.87e-05 4.43e-07 5.61e-09 2.70e-11 + BER + Eb/N0 [dB] soft-decision Viterbi union bound + ---------------------------------------------- + 0 1.401e-01 1.817e+01 + 1 4.320e-02 9.323e-01 + 2 3.625e-03 2.846e-02 + 3 3.500e-04 6.814e-04 + 4 0.000e+00 1.869e-05 + 5 0.000e+00 4.427e-07 + 6 0.000e+00 5.609e-09 + 7 0.000e+00 2.702e-11 .. image:: img/channel_coding_fig1.png :width: 100% @@ -279,17 +299,26 @@ by construction, there is no hard-decision variant of it. .. literalinclude:: ../../examples/simple/channel_coding.py :language: python - :lines: 103-128 + :lines: 108-139 .. code:: LDPC: H is 1020 x 2040, k = 1022 information bits, rate = 0.501, column weight 3, row weight 6 - LDPC 5 iterations 1.70e-01 1.09e-01 3.31e-02 2.30e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 - LDPC 25 iterations 1.70e-01 1.14e-01 8.78e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 + BER + Eb/N0 [dB] LDPC (2040, 1022), 5 iterations LDPC (2040, 1022), 25 iterations + ----------------------------------------------------------------------------- + 0 0.17021 0.16955 + 1 0.10908 0.11424 + 2 0.03312 0.00878 + 3 0.00230 0.00000 + 4 0.00000 0.00000 + 5 0.00000 0.00000 + 6 0.00000 0.00000 + 7 0.00000 0.00000 .. literalinclude:: ../../examples/simple/channel_coding.py :language: python - :lines: 130-143 + :lines: 141-148 .. image:: img/channel_coding_fig2.png :width: 100% diff --git a/docs/tutorials/gn_model.rst b/docs/tutorials/gn_model.rst index 6c7d73c..8071586 100644 --- a/docs/tutorials/gn_model.rst +++ b/docs/tutorials/gn_model.rst @@ -252,8 +252,6 @@ Everything above is a prediction. Now we propagate samples and find out. :language: python :lines: 117-152 -.. mermaid:: mermaid/gn_model.mmd - ``launch_amplitude(power_W, polarizations=2)`` is the factor of two of the warning above, named rather than written: :math:`\sqrt{P/2}` of field per polarization for a channel power :math:`P`. @@ -402,11 +400,6 @@ enhanced GN model exists to recover that gap; this library does not implement it, and the table above is the honest statement of what that costs -- about a decibel for 16QAM, less as the format grows. -.. literalinclude:: ../../examples/optical/gn_model.py - :language: python - :lines: 205-209 - - Going further ^^^^^^^^^^^^^ diff --git a/docs/tutorials/img/monte_carlo_mimo_fig3.png b/docs/tutorials/img/monte_carlo_mimo_fig3.png index 643dd82..ac5adfc 100644 Binary files a/docs/tutorials/img/monte_carlo_mimo_fig3.png and b/docs/tutorials/img/monte_carlo_mimo_fig3.png differ diff --git a/docs/tutorials/img/monte_carlo_mimo_fig4.png b/docs/tutorials/img/monte_carlo_mimo_fig4.png index aa506e0..e189b73 100644 Binary files a/docs/tutorials/img/monte_carlo_mimo_fig4.png and b/docs/tutorials/img/monte_carlo_mimo_fig4.png differ diff --git a/docs/tutorials/img/one_shot_alamouti_fig3.png b/docs/tutorials/img/one_shot_alamouti_fig3.png index 86facc6..adbf8cc 100644 Binary files a/docs/tutorials/img/one_shot_alamouti_fig3.png and b/docs/tutorials/img/one_shot_alamouti_fig3.png differ diff --git a/docs/tutorials/img/one_shot_ofdm_fig5.png b/docs/tutorials/img/one_shot_ofdm_fig5.png index 6309ce4..7b99b1a 100644 Binary files a/docs/tutorials/img/one_shot_ofdm_fig5.png and b/docs/tutorials/img/one_shot_ofdm_fig5.png differ diff --git a/docs/tutorials/mermaid/alamouti.mmd b/docs/tutorials/mermaid/alamouti.mmd deleted file mode 100644 index 2626bed..0000000 --- a/docs/tutorials/mermaid/alamouti.mmd +++ /dev/null @@ -1,20 +0,0 @@ -flowchart LR - tx["SymbolGenerator"] - symbol_mapper["SymbolMapper"] - signal_amplifier["Amplifier"] - space_time_encoder["SpaceTimeEncoder"] - channel["FlatMIMOChannel"] - noise["AWGN"] - detector["SpaceTimeDecoder"] - signal_amplifier_2["Amplifier"] - symbol_demapper["SymbolDemapper"] - tx --> symbol_mapper - symbol_mapper --> signal_amplifier - signal_amplifier --> space_time_encoder - space_time_encoder --> channel - channel --> noise - noise --> detector - detector --> signal_amplifier_2 - signal_amplifier_2 --> symbol_demapper - classDef tapped stroke-dasharray: 4 2 - class tx,noise,detector tapped \ No newline at end of file diff --git a/docs/tutorials/mermaid/awgn_chain.mmd b/docs/tutorials/mermaid/awgn_chain.mmd deleted file mode 100644 index df05c43..0000000 --- a/docs/tutorials/mermaid/awgn_chain.mmd +++ /dev/null @@ -1,10 +0,0 @@ -flowchart LR - tx["SymbolGenerator"] - symbol_mapper["SymbolMapper"] - awgn_channel["AWGN"] - symbol_demapper["SymbolDemapper"] - tx --> symbol_mapper - symbol_mapper --> awgn_channel - awgn_channel --> symbol_demapper - classDef tapped stroke-dasharray: 4 2 - class tx tapped \ No newline at end of file diff --git a/docs/tutorials/mermaid/gn_model.mmd b/docs/tutorials/mermaid/gn_model.mmd deleted file mode 100644 index e5f15b2..0000000 --- a/docs/tutorials/mermaid/gn_model.mmd +++ /dev/null @@ -1,20 +0,0 @@ -flowchart LR - source["SymbolGenerator"] - tx["SymbolMapper"] - launch["Amplifier"] - upsampler["Upsampler"] - srrcfilter["SRRCFilter"] - fibre["FiberLink"] - dbp["DBP"] - srrcfilter_2["SRRCFilter"] - downsampler["Downsampler"] - source --> tx - tx --> launch - launch --> upsampler - upsampler --> srrcfilter - srrcfilter --> fibre - fibre --> dbp - dbp --> srrcfilter_2 - srrcfilter_2 --> downsampler - classDef tapped stroke-dasharray: 4 2 - class tx tapped \ No newline at end of file diff --git a/docs/tutorials/mermaid/mimo_zf.mmd b/docs/tutorials/mermaid/mimo_zf.mmd deleted file mode 100644 index 1d5078d..0000000 --- a/docs/tutorials/mermaid/mimo_zf.mmd +++ /dev/null @@ -1,12 +0,0 @@ -flowchart LR - tx["SymbolGenerator"] - symbol_mapper["SymbolMapper"] - channel["FlatMIMOChannel"] - noise["AWGN"] - detector["LinearDetector"] - tx --> symbol_mapper - symbol_mapper --> channel - channel --> noise - noise --> detector - classDef tapped stroke-dasharray: 4 2 - class tx,noise tapped \ No newline at end of file diff --git a/docs/tutorials/mermaid/mrc.mmd b/docs/tutorials/mermaid/mrc.mmd deleted file mode 100644 index 007d021..0000000 --- a/docs/tutorials/mermaid/mrc.mmd +++ /dev/null @@ -1,12 +0,0 @@ -flowchart LR - tx["SymbolGenerator"] - symbol_mapper["SymbolMapper"] - channel["FlatMIMOChannel"] - noise["AWGN"] - detector["LinearDetector"] - tx --> symbol_mapper - symbol_mapper --> channel - channel --> noise - noise --> detector - classDef tapped stroke-dasharray: 4 2 - class tx tapped \ No newline at end of file diff --git a/docs/tutorials/mermaid/shaping_pas.mmd b/docs/tutorials/mermaid/shaping_pas.mmd deleted file mode 100644 index 4d39eb9..0000000 --- a/docs/tutorials/mermaid/shaping_pas.mmd +++ /dev/null @@ -1,6 +0,0 @@ -flowchart LR - bits["SymbolGenerator"] - dm["DistributionMatcher"] - signs["AmplitudeMapper"] - bits --> dm - dm --> signs \ No newline at end of file diff --git a/docs/tutorials/mimo.rst b/docs/tutorials/mimo.rst index d9c13ed..6507eef 100644 --- a/docs/tutorials/mimo.rst +++ b/docs/tutorials/mimo.rst @@ -7,9 +7,10 @@ compare on it. .. note:: - **Before you start.** :doc:`awgn` introduced ``monte_carlo``, which is used - here to average over channel realizations rather than over noise. The - chain is the same object; only the number of antennas changes. + **Before you start.** :doc:`awgn` introduced ``monte_carlo``. Here the + averaging over fading is a *batch*: a fixed stack of channel draws is + built into each chain, and the sweep only moves the noise. The chain + is the same object; only the number of antennas changes. **What you'll learn:** @@ -45,42 +46,24 @@ problem, and it is what the rest of this tutorial is about. Implementation """""""""""""" -We start with the imports: +We start with the imports and the parameters: .. literalinclude:: ../../examples/mimo/one_shot_mimo.py :language: python - :lines: 1-21 + :lines: 1-32 -The channel is drawn once, with a seed: the first half of the tutorial is -about **one** realization, so it must be the same one on every run. The link -is a single ``Sequential`` closed by a detector, so comparing five detectors -means comparing five chains that differ by their last block alone: +The channel is drawn once, with a seed: this first half of the tutorial is +about **one** realization, so it must be the same one on every run. One +chain, transmitter to decision, closed by the most direct detector -- zero +forcing: .. literalinclude:: ../../examples/mimo/one_shot_mimo.py :language: python - :lines: 23-63 - -.. mermaid:: mermaid/mimo_zf.mmd - -The diagram is not drawn by hand. It is what the chain says about itself -- -``chain.to_mermaid()`` (decision D33c) -- exported by the script, so the -block names are the ones the code uses and a dashed outline marks a tapped -block: - -.. literalinclude:: ../../examples/mimo/one_shot_mimo.py - :language: python - :lines: 185-188 - -Each chain is given the same seed before running, so the five numbers below -differ by the detector alone -- same symbols, same noise, same channel: + :lines: 34-47 .. code:: - * detector ZF : ser=0.0025 - * detector MMSE : ser=0.0025 - * detector OSIC : ser=0.0010 - * detector ML : ser=0.0005 - * detector SD : ser=0.0005 + ZF, one channel draw: SER = 0.0025 Received Constellations """"""""""""""""""""""" @@ -89,21 +72,18 @@ Let us look at what each receive antenna sees, read from the ``"noise"`` tap: .. literalinclude:: ../../examples/mimo/one_shot_mimo.py :language: python - :lines: 65-72 + :lines: 49-56 .. image:: img/monte_carlo_mimo_fig1.png :width: 100% :align: center -No constellation is visible on any antenna: each one receives a noisy -superposition of the two transmitted streams. This is the problem the -detector has to solve. - -Applying the pseudo-inverse of the channel matrix separates the streams: +No constellation is visible on any antenna: each one observes a mixture of +the two streams. The ZF estimator inverts that mixture: .. literalinclude:: ../../examples/mimo/one_shot_mimo.py :language: python - :lines: 74-81 + :lines: 58-65 .. image:: img/monte_carlo_mimo_fig2.png :width: 100% @@ -119,7 +99,19 @@ Detection Strategies ^^^^^^^^^^^^^^^^^^^^ The five detectors are five answers to one question: what to do with the -interference that the other streams put on top of the one being read. +interference that the other streams put on top of the one being read. One +channel realization proves nothing about them -- over fading, the error rate +is an *average*, dominated by the rare draws where the matrix is nearly +singular -- so each answer is judged over the **same** 200 channel draws, +stacked into a batch: one draw is one row (D51), the channel block +propagates draw :math:`k` on frame :math:`k`, and each detector holds the +same stack. The draws and the sweep are set up once: + +.. literalinclude:: ../../examples/mimo/one_shot_mimo.py + :language: python + :lines: 67-81 + +Each technique is then three things -- an idea, one chain, one sweep. Zero forcing (ZF) """"""""""""""""" @@ -141,6 +133,10 @@ enhancement**. It also costs diversity: each stream spends :math:`N_t-1` of its :math:`N_r` degrees of freedom cancelling the others, leaving a diversity order of :math:`N_r - N_t + 1` instead of :math:`N_r`. +.. literalinclude:: ../../examples/mimo/one_shot_mimo.py + :language: python + :lines: 83-94 + Minimum mean square error (MMSE) """""""""""""""""""""""""""""""" @@ -159,7 +155,13 @@ conditioned. As :math:`\sigma^2 \to 0` it *is* zero forcing; as :math:`\sigma^2 \to \infty` it becomes the matched filter :math:`\mathbf{H}^H`. In between it accepts a little residual interference in exchange for much less amplified noise. It does not buy diversity, which is -why the ZF and MMSE curves run parallel. +why the ZF and MMSE curves run parallel. Because the detector *weights by* +:math:`\sigma^2`, the sweep zips the noise variance into it, one pair of +dotted parameters per point: + +.. literalinclude:: ../../examples/mimo/one_shot_mimo.py + :language: python + :lines: 96-107 Ordered successive interference cancellation (OSIC) """"""""""""""""""""""""""""""""""""""""""""""""""" @@ -173,6 +175,10 @@ is why the ordering matters and why ``osic_type="sinr"`` sorts by post-detection SNR (the V-BLAST rule). The price is error propagation: a wrong decision is subtracted as if it were right. +.. literalinclude:: ../../examples/mimo/one_shot_mimo.py + :language: python + :lines: 109-121 + Maximum likelihood (ML) """"""""""""""""""""""" @@ -187,7 +193,12 @@ refuses that split and scores the candidate vectors jointly: Nothing is inverted, so nothing is amplified, and every receive antenna contributes to every stream: the diversity order is :math:`N_r`. The cost is that the minimum is taken over :math:`|\mathcal{M}|^{N_t}` vectors -- 16 here, -but 65 536 for 16-QAM on four streams. +but 65 536 for 16-QAM on four streams. ML compares distances and never needs +:math:`\sigma^2`, so its sweep moves the noise alone: + +.. literalinclude:: ../../examples/mimo/one_shot_mimo.py + :language: python + :lines: 123-133 Sphere decoding (SD) """""""""""""""""""" @@ -211,66 +222,57 @@ can be discarded**: it cannot contain the minimum. Only the lattice points inside a shrinking ball are ever visited, and nothing is approximated, so the decision is exactly the ML one. +.. literalinclude:: ../../examples/mimo/one_shot_mimo.py + :language: python + :lines: 135-145 + .. note :: Diversity orders are claims about a limit, and the SNR range simulated - below stops at 18 dB, which is not the asymptotic regime. They are + here stops at 18 dB, which is not the asymptotic regime. They are checked where a limit can be reached: ``validation/mimo_zf_ml_ber.py`` confirms diversity 1 for zero forcing and 2 for maximum likelihood on a 2x2 channel, against the closed forms of :mod:`comnumpy.core.metrics` (decision D7). -Monte Carlo Evaluation -^^^^^^^^^^^^^^^^^^^^^^ +Results +^^^^^^^ -A single channel realization proves nothing: over fading, the error rate is -an *average*, dominated by the rare draws where the matrix is nearly -singular. Averaging means one independent draw per frame, and the draws are -a **stack**: ``rayleigh_channel(size=n_channels)`` returns ``(200, 3, 2)``, -the channel block propagates frame :math:`k` through matrix :math:`k`, and -the detector -- handed the *same* stack -- decides frame :math:`k` against -matrix :math:`k`. The chain runs once per SNR point instead of once per -draw, and the pooled SER over equal-size frames is exactly the mean of the -per-draw rates: +The five curves come from the same dictionary, shown twice: .. literalinclude:: ../../examples/mimo/one_shot_mimo.py :language: python - :lines: 83-132 + :lines: 147-157 .. code:: SER - snr_dB ZF MMSE OSIC ML SD - --------------------------------------------------------- - 0 0.33086 0.28510 2.822e-01 2.791e-01 2.791e-01 - 3 0.20976 0.17386 1.595e-01 1.545e-01 1.545e-01 - 6 0.10669 0.08601 6.488e-02 5.980e-02 5.980e-02 - 9 0.04389 0.03476 1.890e-02 1.584e-02 1.584e-02 - 12 0.01510 0.01174 4.163e-03 2.488e-03 2.488e-03 - 15 0.00507 0.00394 8.000e-04 3.875e-04 3.875e-04 - 18 0.00201 0.00144 1.625e-04 1.625e-04 1.625e-04 + snr_dB ZF MMSE OSIC ML SD + ------------------------------------------------------- + 0 0.33427 0.28820 0.28735 2.851e-01 2.868e-01 + 3 0.21634 0.17955 0.16855 1.613e-01 1.620e-01 + 6 0.11209 0.09271 0.07311 6.785e-02 6.871e-02 + 9 0.04911 0.04034 0.02546 2.057e-02 2.095e-02 + 12 0.01876 0.01570 0.00779 4.700e-03 4.950e-03 + 15 0.00719 0.00583 0.00243 8.125e-04 6.875e-04 + 18 0.00286 0.00196 0.00100 1.000e-04 2.500e-05 .. image:: img/monte_carlo_mimo_fig3.png :width: 100% :align: center The ordering is the textbook one and it holds at every SNR: ML is the best, -OSIC follows closely, MMSE beats ZF everywhere, and the gap widens with the -SNR -- at 18 dB, ML is an order of magnitude below ZF. The reason is the -slope: with :math:`N_r = 3` receive antennas and :math:`N_t = 2` streams, ML -enjoys the full receive diversity while a linear detector spends part of it -cancelling the other stream. - -The sphere decoder's row is **the ML row, digit for digit**, which is the -claim it has to make: a pruning that changed a single decision would not be -maximum likelihood any more. The two curves therefore lie on top of each -other in the figure, and only their markers tell them apart. +OSIC follows, MMSE beats ZF everywhere, and the gap widens with the SNR -- +at 18 dB, ML is an order of magnitude below ZF. The reason is the slope: +with :math:`N_r = 3` receive antennas and :math:`N_t = 2` streams, ML enjoys +the full receive diversity while a linear detector spends part of it +cancelling the other stream. The ML and SD rows differ only by their noise +draws (each sweep reseeds its own chain): where they share a seed they are +digit for digit the same decision, which is the sphere decoder's claim. -Note which detectors are told the noise variance and which are not: MMSE and -OSIC weight by :math:`\sigma^2`, ZF ignores the noise by construction, and ML -only compares distances -- so neither of the last two has that parameter at -all, and ``set_params`` says so if you try. +The last rows sit near the estimator's floor: 2.5e-05 out of 40 000 symbols +per point is one error, not a rate. Computational Cost @@ -284,16 +286,16 @@ where that product no longer fits. We therefore time both detectors on .. literalinclude:: ../../examples/mimo/one_shot_mimo.py :language: python - :lines: 134-183 + :lines: 159-223 .. code:: - ML 0 dB 142.0 ms 65536 nodes SER 0.7900 - ML 9 dB 123.5 ms 65536 nodes SER 0.5044 - ML 18 dB 140.8 ms 65536 nodes SER 0.0112 - SD 0 dB 197.5 ms 91.8 nodes SER 0.7900 - SD 9 dB 36.3 ms 19.2 nodes SER 0.5044 - SD 18 dB 11.8 ms 5.8 nodes SER 0.0112 + ML 0 dB 176.7 ms 65536 nodes SER 0.7900 + ML 9 dB 113.7 ms 65536 nodes SER 0.5044 + ML 18 dB 119.0 ms 65536 nodes SER 0.0112 + SD 0 dB 163.4 ms 91.8 nodes SER 0.7900 + SD 9 dB 33.3 ms 19.2 nodes SER 0.5044 + SD 18 dB 10.8 ms 5.8 nodes SER 0.0112 .. image:: img/monte_carlo_mimo_fig4.png :width: 100% diff --git a/docs/tutorials/ofdm.rst b/docs/tutorials/ofdm.rst index 410a5e8..7cfe675 100644 --- a/docs/tutorials/ofdm.rst +++ b/docs/tutorials/ofdm.rst @@ -146,7 +146,7 @@ Results .. code:: - single carrier: SER 0.0117, 3669 ms + single carrier: SER 0.0117, 3713 ms .. image:: img/one_shot_ofdm_fig2.png :width: 100% @@ -202,7 +202,7 @@ Results .. code:: - OFDM : SER 0.0406, 1.00 ms (3670 times faster) + OFDM : SER 0.0406, 1.06 ms (3494 times faster) .. image:: img/one_shot_ofdm_fig3.png :width: 100% @@ -221,7 +221,7 @@ range of SNR values: .. literalinclude:: ../../examples/ofdm/one_shot_ofdm.py :language: python - :lines: 101-128 + :lines: 101-130 .. code:: @@ -268,17 +268,17 @@ block length grows: .. literalinclude:: ../../examples/ofdm/one_shot_ofdm.py :language: python - :lines: 129-165 + :lines: 131-166 .. code:: receiver runtime [ms], and their ratio block length N single carrier OFDM ratio --------------------------------------------- - 128 7.34 0.568 12.9 - 256 34.91 0.567 61.6 - 512 188.09 0.793 237.1 - 1024 1377.50 0.862 1598.3 + 128 8.16 0.805 10.1 + 256 33.58 0.733 45.8 + 512 219.69 0.890 246.8 + 1024 1491.35 0.945 1577.4 .. image:: img/one_shot_ofdm_fig5.png :width: 100% diff --git a/docs/tutorials/ofdm_papr.rst b/docs/tutorials/ofdm_papr.rst index f2dd948..cc12d84 100644 --- a/docs/tutorials/ofdm_papr.rst +++ b/docs/tutorials/ofdm_papr.rst @@ -45,14 +45,14 @@ are simply not seen. .. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py :language: python - :lines: 1-44 + :lines: 1-41 Let us generate four OFDM symbols and look at the instantaneous power, normalized by its own mean: .. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py :language: python - :lines: 47-63 + :lines: 43-57 .. image:: img/monte_carlo_ofdm_papr_fig1.png :width: 100% @@ -81,7 +81,7 @@ non-zero probability. .. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py :language: python - :lines: 69-95 + :lines: 59-89 .. image:: img/monte_carlo_ofdm_papr_fig2.png :width: 100% @@ -109,7 +109,7 @@ with those values, which saves reducing them by hand: .. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py :language: python - :lines: 98-104 + :lines: 91-98 .. code:: @@ -150,20 +150,10 @@ the oversampling rather than the effective count: the fit belongs to the library, with its domain of validity, rather than to each script that draws the curve. -The fitted constant can be avoided. Counting how often a Gaussian process -crosses a level gives an effective count that grows with the threshold -instead, - -.. math:: - - \mathrm{CCDF}(\gamma) \simeq 1 - \exp\left(-N_{sc} - \sqrt{\frac{\pi}{3}}\, \sqrt{\gamma}\, e^{-\gamma}\right) - -which the same function returns with ``method="level_crossing"``. It has no -fitted constant, but it describes the peak of the **continuous-time** -waveform, which a sampled one can only underestimate -- so it reads high on -sampled data, and it is a large-threshold approximation. The two are -compared against the measurement at the end of this page. +A fit-free alternative based on level-crossing counts exists -- the same +function returns it with ``method="level_crossing"`` -- but it describes the +peak of the **continuous-time** waveform, which sampled data can only +underestimate, so it reads high here. Implementation """""""""""""" @@ -172,7 +162,7 @@ We estimate the CCDF over 20 000 OFDM symbols, for 256 and 1024 subcarriers: .. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py :language: python - :lines: 109-145 + :lines: 100-152 Results """"""" @@ -184,34 +174,13 @@ Results .. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py :language: python - :lines: 146-150 + :lines: 154-158 .. code:: N_sc = 256: PAPR exceeded once in a thousand symbols above 11.30 dB N_sc = 1024: PAPR exceeded once in a thousand symbols above 11.72 dB -Which of the two models to believe is a question the measurement answers: - -.. literalinclude:: ../../examples/ofdm/monte_carlo_ofdm_papr.py - :language: python - :lines: 152-166 - -.. code:: - - N_sc level threshold effective level crossing - 256 1e-02 10.43 dB 1.2e-02 1.4e-02 - 256 1e-03 11.28 dB 1.1e-03 1.4e-03 - 1024 1e-02 11.00 dB 9.6e-03 1.2e-02 - 1024 1e-03 11.69 dB 1.1e-03 1.6e-03 - -The effective-count model is within 20 % of the measurement, the -level-crossing one 20 to 60 % above it -- consistently above, which is the -direction it must err in, since it is the continuous-time waveform that is -being modelled and the simulation only sees four samples per Nyquist -interval. Take the first to size an amplifier from sampled data, and the -second when the question is the analogue waveform itself. - The measured points sit on the closed form over the whole range. As expected, the probability of a large PAPR grows with the number of subcarriers -- but slowly: four times as many subcarriers cost only 0.4 dB, because the number diff --git a/docs/tutorials/shaping.rst b/docs/tutorials/shaping.rst index ba8e07a..335755d 100644 --- a/docs/tutorials/shaping.rst +++ b/docs/tutorials/shaping.rst @@ -46,7 +46,7 @@ variance -- and the Gaussian is the maximum-entropy law at fixed variance. .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 1-32 + :lines: 1-31 Why a Uniform QAM Is Not Optimal @@ -59,7 +59,7 @@ then the mutual information :math:`I(X;Y)` of that discrete input, which .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 34-62 + :lines: 33-61 .. image:: img/probabilistic_shaping_fig1.png :width: 100% @@ -117,7 +117,7 @@ transmit power**, which is the only comparison that means anything. .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 64-119 + :lines: 63-118 .. image:: img/probabilistic_shaping_fig2.png :width: 100% @@ -156,7 +156,7 @@ that lives there: .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 121-133 + :lines: 120-132 .. code:: @@ -189,7 +189,7 @@ the rate is :math:`k/n`. As :math:`n` grows, :math:`k/n \to H(X)`; at finite .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 135-148 + :lines: 134-147 .. code:: @@ -218,9 +218,7 @@ compose instead of fighting. .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 150-162 - -.. mermaid:: mermaid/shaping_pas.mmd + :lines: 149-159 Uniform bits go in, shaped symbols come out, so the check is to run it and count what comes out. The signs are equiprobable, so the law the signed @@ -235,7 +233,7 @@ which is the Maxwell-Boltzmann law on the full constellation at the same .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 164-191 + :lines: 161-188 .. image:: img/probabilistic_shaping_fig3.png :width: 100% @@ -263,7 +261,7 @@ the generalized mutual information, .. literalinclude:: ../../examples/simple/probabilistic_shaping.py :language: python - :lines: 193-205 + :lines: 190-202 .. code:: diff --git a/examples/mimo/one_shot_alamouti.py b/examples/mimo/one_shot_alamouti.py index 43af67f..12d6a10 100644 --- a/examples/mimo/one_shot_alamouti.py +++ b/examples/mimo/one_shot_alamouti.py @@ -1,12 +1,12 @@ """Rayleigh fading, receive diversity, and the Alamouti space-time code. -Run from this directory: it writes the tutorial's figures and diagrams -into ../../docs/tutorials/. +Run from this directory: it writes the tutorial's figures into +../../docs/tutorials/. """ import matplotlib.pyplot as plt import numpy as np -from comnumpy import monte_carlo, style +from comnumpy import monte_carlo, print_data, style from comnumpy.core import Sequential from comnumpy.core.generators import SymbolGenerator from comnumpy.core.mappers import SymbolDemapper, SymbolMapper @@ -22,13 +22,11 @@ style.use() img_dir = "../../docs/tutorials/img/" -mermaid_dir = "../../docs/tutorials/mermaid/" constellation = Constellation("PSK", 4) M = constellation.order code = get_code("alamouti") power = 1 / np.sqrt(code.n_tx) # split the power over the antennas -sigma2 = 0.2 # --- the fading channel ---------------------------------------------- # One Rayleigh coefficient is a complex Gaussian, so its squared modulus @@ -54,50 +52,29 @@ print(f"P[|h|^2 < {threshold:g}] = {np.mean(gain < threshold):.4f} " f"(1 - exp(-{threshold:g}) = {1 - np.exp(-threshold):.4f})") -# --- the three links -------------------------------------------------- - - -def get_chain(kind, H, sigma2=sigma2): - """Build the chain of one scheme, on a given channel realization. - - The three schemes differ by two blocks -- what is put on the - antennas, and what is done with what comes back -- so they are one - function and not three scripts. - """ - source = [SymbolGenerator(constellation.order, name="tx"), - SymbolMapper(constellation)] - channel = [FlatMIMOChannel(H, name="channel"), - AWGN(sigma2=sigma2, name="noise")] - match kind: - case "alamouti": - return Sequential([ - *source, Amplifier(power), SpaceTimeEncoder(code), *channel, - SpaceTimeDecoder(code, H=H, name="detector"), - Amplifier(1 / power), SymbolDemapper(constellation), - ], taps=["tx", "noise", "detector"], name="Alamouti 2x1") - case "linear": - # zero forcing on an (N_r, 1) channel *is* maximum ratio - # combining: the pseudo-inverse of a column is h^H / ||h||^2 - return Sequential([ - *source, *channel, - LinearDetector(constellation, H=H, name="detector"), - ], taps=["tx"], name=f"{H.shape[0]} Rx, {H.shape[1]} Tx") - case _: - raise ValueError(f"unknown link {kind!r}") - - -alamouti = get_chain("alamouti", rayleigh_channel(1, 2, seed=42)) -siso = get_chain("linear", rayleigh_channel(1, 1, seed=1)) -mrc = get_chain("linear", rayleigh_channel(2, 1, seed=2)) - # --- one shot --------------------------------------------------------- -alamouti.seed(7) # every stochastic block, reproducibly -s_rx = alamouti(500 * code.n_symbols) -print(f"\none-shot SER: {compute_ser(alamouti.tap('tx'), s_rx):.4f}") -alamouti.summary(500 * code.n_symbols) - -received = alamouti.tap("noise")[0] -combined = alamouti.tap("detector") / power +# One channel draw, one chain, transmitter to decision: what Alamouti +# combining does to a single receive antenna. +H_one = rayleigh_channel(1, 2, seed=42) +one_shot = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + Amplifier(power), + SpaceTimeEncoder(code), + FlatMIMOChannel(H_one, name="channel"), + AWGN(sigma2=0.2, name="noise"), + SpaceTimeDecoder(code, H=H_one, name="detector"), + Amplifier(1 / power), + SymbolDemapper(constellation), +], taps=["tx", "noise", "detector"], name="Alamouti 2x1") + +one_shot.seed(7) # every stochastic block, reproducibly +s_rx = one_shot(500 * code.n_symbols) +print(f"\none-shot SER: {compute_ser(one_shot.tap('tx'), s_rx):.4f}") +one_shot.summary(500 * code.n_symbols) + +received = one_shot.tap("noise")[0] +combined = one_shot.tap("detector") / power fig2, (ax_left, ax_right) = plt.subplots(nrows=1, ncols=2, figsize=(9, 4.2)) plot_iq(received, marker=".", ax=ax_left) ax_left.set_title("tap('noise'): the single receive antenna") @@ -107,51 +84,68 @@ def get_chain(kind, H, sigma2=sigma2): plt.savefig(f"{img_dir}/one_shot_alamouti_fig2.png") # --- averaging over fading -------------------------------------------- - - -def average_ser(chain, n_rx, n_tx, snr_dB, stimulus, n_channels, seed=0): - """Average one chain over independent quasi-static fading draws. - - The chain is built once. ``monte_carlo`` takes several dotted parameters - at once and zips them, so one sweep point sets the channel the - signal goes through *and* the channel the detector inverts -- which - is exactly what a fading realization is. - """ - rng = np.random.default_rng(seed) - draws = [] - for _ in range(n_channels): - H = rayleigh_channel(n_rx, n_tx, rng=rng) - draws.append((H, H)) - chain.set_params(noise__sigma2=10 ** (-snr_dB / 10)) - results = monte_carlo(chain, ("channel.H", "detector.H"), draws, - {"ser": compute_ser}, stimulus=stimulus, - reference="tx", seed=seed) - return float(np.mean(results["ser"])) - - -# The accuracy of an average over fading is set by the number of channel -# draws, not by the symbol count, so the draw count grows with the SNR. -snr_dB_list = np.arange(4, 25, 4) -draws = [1500, 2000, 2500, 3500, 5000, 6000] +# The SAME K channel draws serve every SNR point: one draw is one row +# of the batch (D51), each chain is built once with its stack -- the +# channel block propagates draw k on frame k, the detector holds the +# same stack -- and only the noise moves along the sweep. +K = 5000 n_symbols = 80 -links = { - "1 Tx, 1 Rx (no diversity)": (siso, 1, 1, (1, n_symbols)), - "Alamouti, 2 Tx, 1 Rx": (alamouti, 1, 2, n_symbols), - "MRC, 1 Tx, 2 Rx": (mrc, 2, 1, (1, n_symbols)), -} -# --- metrics, pre-allocated ------------------------------------------ -curves = {} -for name in links: - curves[name] = np.zeros(len(snr_dB_list)) - -# --- simulation loop ------------------------------------------------- -for name, (chain, n_rx, n_tx, stimulus) in links.items(): - for index, (snr_dB, count) in enumerate(zip(snr_dB_list, draws, - strict=True)): - curves[name][index] = average_ser(chain, n_rx, n_tx, snr_dB, - stimulus, count) +rng = np.random.default_rng(3) +H_siso = rayleigh_channel(1, 1, rng=rng, size=K) +H_alamouti = rayleigh_channel(1, 2, rng=rng, size=K) +H_mrc = rayleigh_channel(2, 1, rng=rng, size=K) + +siso = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_siso, name="channel"), + AWGN(sigma2=1.0, name="noise"), + # zero forcing on an (N_r, 1) channel *is* maximum ratio combining: + # the pseudo-inverse of a column is h^H / ||h||^2 + LinearDetector(constellation, H=H_siso, name="detector"), +], taps=["tx"], name="1 Tx, 1 Rx") + +alamouti = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + Amplifier(power), + SpaceTimeEncoder(code), + FlatMIMOChannel(H_alamouti, name="channel"), + AWGN(sigma2=1.0, name="noise"), + SpaceTimeDecoder(code, H=H_alamouti, name="detector"), + Amplifier(1 / power), + SymbolDemapper(constellation), +], taps=["tx"], name="Alamouti 2x1") + +mrc = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_mrc, name="channel"), + AWGN(sigma2=1.0, name="noise"), + LinearDetector(constellation, H=H_mrc, name="detector"), +], taps=["tx"], name="MRC 1x2") + +# --- the sweep: no simulation loop ----------------------------------- +# monte_carlo moves the noise variance; everything else is frozen in +# the chains above. +snr_dB_list = np.arange(4, 25, 4) +sigma2_list = 10.0 ** (-snr_dB_list / 10) -# --- results: figure against the closed forms ------------------------ +curves = {} +curves["1 Tx, 1 Rx (no diversity)"] = monte_carlo( + siso, "noise.sigma2", sigma2_list, {"ser": compute_ser}, + (K, 1, n_symbols), reference="tx", seed=1)["ser"] +curves["Alamouti, 2 Tx, 1 Rx"] = monte_carlo( + alamouti, "noise.sigma2", sigma2_list, {"ser": compute_ser}, + (K, n_symbols), reference="tx", seed=2)["ser"] +curves["MRC, 1 Tx, 2 Rx"] = monte_carlo( + mrc, "noise.sigma2", sigma2_list, {"ser": compute_ser}, + (K, 1, n_symbols), reference="tx", seed=3)["ser"] + +# --- results: table and figure against the closed forms -------------- +print() +print_data({"x": snr_dB_list, "curves": curves}, + xlabel="snr_dB", ylabel="SER") fine = np.linspace(snr_dB_list[0], snr_dB_list[-1], 100) per_bit = 10 ** (fine / 10) / np.log2(M) @@ -163,48 +157,8 @@ def average_ser(chain, n_rx, n_tx, snr_dB, stimulus, n_channels, seed=0): } ax = plot_error_rate(snr_dB_list, curves, theory=theory, x_theory=fine, ylabel="symbol error rate", - title=f"{M}-PSK over Rayleigh, equal transmit power") + title=f"{M}-PSK over Rayleigh, {K} channel draws, " + f"equal transmit power") ax.set_ylim(1e-5, 1) plt.tight_layout() plt.savefig(f"{img_dir}/one_shot_alamouti_fig3.png") - -# --- the two claims, read off the closed form ------------------------- -def closed_form(name, snr_dB): - """The Rayleigh SER this curve should follow, at equal transmit power.""" - power_split = code.n_tx if "Alamouti" in name else 1 - diversity = 1 if "no diversity" in name else 2 - return compute_ser_rayleigh_psk( - M, 10 ** (snr_dB / 10) / np.log2(M) / power_split, - diversity=diversity) - - -for name, values in curves.items(): - ratio = np.array(values) / closed_form(name, snr_dB_list) - line = f"{name:28s} measured / closed form " - for value in ratio: - line += f"{value:.2f} " - print(line) - -high = np.array([30.0, 40.0]) -for name in curves: - reference = closed_form(name, high) - slope = np.polyfit(high / 10, np.log10(reference), 1)[0] - print(f"{name:28s} diversity order {-slope:.2f}") - -grid_dB = np.linspace(0, 30, 601) -alamouti_snr = 10 ** (grid_dB / 10) / np.log2(M) -target = 1e-3 -of = {"MRC": compute_ser_rayleigh_psk(M, alamouti_snr, diversity=2), - "Alamouti": compute_ser_rayleigh_psk(M, alamouti_snr / code.n_tx, - diversity=2)} -needed = {} -for name, values in of.items(): - needed[name] = np.interp(np.log10(target), np.log10(values)[::-1], - grid_dB[::-1]) -print(f"SNR for SER = {target:g}: MRC {needed['MRC']:.1f} dB, Alamouti " - f"{needed['Alamouti']:.1f} dB, gap {needed['Alamouti'] - needed['MRC']:.2f} dB " - f"(10log10(N_t) = {10 * np.log10(code.n_tx):.2f} dB)") - -for diagram_name, diagram_chain in [("alamouti", alamouti), ("mrc", mrc)]: - with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: - stream.write(diagram_chain.to_mermaid()) diff --git a/examples/mimo/one_shot_mimo.py b/examples/mimo/one_shot_mimo.py index 2bd1b81..9167e70 100644 --- a/examples/mimo/one_shot_mimo.py +++ b/examples/mimo/one_shot_mimo.py @@ -1,9 +1,14 @@ +"""MIMO detection: one spatial mixture, five ways to undo it. + +Run from this directory: it writes the tutorial's figures into +../../docs/tutorials/. +""" import time import numpy as np import matplotlib.pyplot as plt -from comnumpy import print_data, style +from comnumpy import monte_carlo, print_data, style from comnumpy.core import Sequential from comnumpy.core.generators import SymbolGenerator from comnumpy.core.mappers import SymbolMapper @@ -23,46 +28,25 @@ N = 1000 N_r, N_t = 3, 2 constellation = Constellation("PSK", 4) +M = constellation.order sigma2 = 0.1 -H = rayleigh_channel(N_r, N_t, seed=0) - -def get_chain(detector): - """The same link, closed by one detector or another. - - Only the last block changes, which is the point: a detector is a - chain block like any other, and the comparison further down is four - chains that differ by that block alone. - """ - return Sequential([ - SymbolGenerator(constellation.order, name="tx"), - SymbolMapper(constellation), - FlatMIMOChannel(H, name="channel"), - AWGN(sigma2=sigma2, name="noise"), - detector, - ], taps=["tx", "noise"], name=f"{N_r}x{N_t} MIMO, {detector.name}") - - -detectors = { - "ZF": LinearDetector(constellation, H=H, method="zf", name="detector"), - "MMSE": LinearDetector(constellation, H=H, sigma2=sigma2, method="mmse", - name="detector"), - "OSIC": OrderedSuccessiveInterferenceCancellationDetector( - constellation, osic_type="sinr", H=H, sigma2=sigma2, - name="detector"), - "ML": MaximumLikelihoodDetector(constellation, H=H, name="detector"), - "SD": SphereDecoder(constellation, H=H, name="detector"), -} -chains = {} -for name, detector in detectors.items(): - chains[name] = get_chain(detector) - -for name, chain in chains.items(): - chain.seed(0) - detected = chain((N_t, N)) - print(f"* detector {name:5s}: ser={compute_ser(chain.tap('tx'), detected):.4f}") - -Y = chains["ZF"].tap("noise") +# --- the problem: one channel draw, one chain, zero forcing ----------- +H = rayleigh_channel(N_r, N_t, seed=0) +chain = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H, name="channel"), + AWGN(sigma2=sigma2, name="noise"), + LinearDetector(constellation, H=H, method="zf", name="detector"), +], taps=["tx", "noise"], name=f"{N_r}x{N_t} MIMO, ZF") + +chain.seed(0) +detected = chain((N_t, N)) +print(f"ZF, one channel draw: SER = " + f"{compute_ser(chain.tap('tx'), detected):.4f}") + +Y = chain.tap("noise") fig1, axes1 = plt.subplots(nrows=1, ncols=N_r, figsize=(4 * N_r, 4)) for index in range(N_r): plot_iq(Y[index, :], ax=axes1[index]) @@ -71,7 +55,7 @@ def get_chain(detector): axes1[index].set_ylim([-2, 2]) plt.savefig(f"{img_dir}/monte_carlo_mimo_fig1.png") -Z = detectors["ZF"].linear_estimator(Y) +Z = chain["detector"].linear_estimator(Y) fig2, axes2 = plt.subplots(nrows=1, ncols=N_t, figsize=(4 * N_t, 4)) for index in range(N_t): plot_iq(Z[index, :], reference=constellation, ax=axes2[index]) @@ -80,94 +64,150 @@ def get_chain(detector): axes2[index].set_ylim([-2, 2]) plt.savefig(f"{img_dir}/monte_carlo_mimo_fig2.png") -snr_dB_list = np.arange(0, 20, 3) -n_channels = 200 +# --- the detectors, averaged over fading ------------------------------ +# The SAME K channel draws serve every SNR point: one draw is one row +# of the batch (D51). Each technique below is one chain built with the +# stack -- on the channel block and on the detector -- followed by its +# own Monte-Carlo sweep: monte_carlo moves the noise variance, and the +# detectors that weight by sigma2 receive it too, through the zip of +# two parameter paths. No simulation loop. +K = 200 n_symbols = 200 - - -NOISE_AWARE = {"MMSE", "OSIC"} # the two detectors that weight by sigma2 - - -def average_ser(name, chain, snr_dB, seed=0): - """Average one chain over independent Rayleigh draws at one SNR. - - The draws are a *stack* (D51): channel k propagates frame k, the - detector holds the same stack and decides frame k against channel - k, and the chain runs once instead of once per draw. Pooled over - equal-size frames, the SER is exactly the mean of the per-draw - rates. - """ - stacked = rayleigh_channel(N_r, N_t, seed=seed, size=n_channels) - noise_variance = N_t * 10 ** (-snr_dB / 10) - params = {"noise.sigma2": noise_variance, - "channel.H": stacked, "detector.H": stacked} - if name in NOISE_AWARE: - params["detector.sigma2"] = noise_variance - chain.set_params(**params) - chain.seed(seed) - detected = chain((n_channels, N_t, n_symbols)) - return float(compute_ser(chain.tap("tx"), detected)) - - -# --- metrics, pre-allocated ------------------------------------------ +snr_dB_list = np.arange(0, 20, 3) +sigma2_list = N_t * 10.0 ** (-snr_dB_list / 10) +both = list(zip(sigma2_list, sigma2_list, strict=True)) +stimulus = (K, N_t, n_symbols) +H_stack = rayleigh_channel(N_r, N_t, seed=1, size=K) curves = {} -for name in chains: - curves[name] = np.zeros(len(snr_dB_list)) -# --- simulation loop ------------------------------------------------- -for name, chain in chains.items(): - for index, snr_dB in enumerate(snr_dB_list): - curves[name][index] = average_ser(name, chain, snr_dB) +# --- zero forcing ---------------------------------------------------- +zf = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_stack, name="channel"), + AWGN(sigma2=1.0, name="noise"), + LinearDetector(constellation, H=H_stack, method="zf", + name="detector"), +], taps=["tx"], name="ZF") +curves["ZF"] = monte_carlo( + zf, "noise.sigma2", sigma2_list, {"ser": compute_ser}, + stimulus, reference="tx", seed=1)["ser"] + +# --- MMSE ------------------------------------------------------------ +mmse = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_stack, name="channel"), + AWGN(sigma2=1.0, name="noise"), + LinearDetector(constellation, H=H_stack, method="mmse", sigma2=1.0, + name="detector"), +], taps=["tx"], name="MMSE") +curves["MMSE"] = monte_carlo( + mmse, ("noise.sigma2", "detector.sigma2"), both, {"ser": compute_ser}, + stimulus, reference="tx", seed=2)["ser"] + +# --- OSIC ------------------------------------------------------------ +osic = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_stack, name="channel"), + AWGN(sigma2=1.0, name="noise"), + OrderedSuccessiveInterferenceCancellationDetector( + constellation, osic_type="sinr", H=H_stack, sigma2=1.0, + name="detector"), +], taps=["tx"], name="OSIC") +curves["OSIC"] = monte_carlo( + osic, ("noise.sigma2", "detector.sigma2"), both, {"ser": compute_ser}, + stimulus, reference="tx", seed=3)["ser"] + +# --- maximum likelihood ---------------------------------------------- +ml = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_stack, name="channel"), + AWGN(sigma2=1.0, name="noise"), + MaximumLikelihoodDetector(constellation, H=H_stack, name="detector"), +], taps=["tx"], name="ML") +curves["ML"] = monte_carlo( + ml, "noise.sigma2", sigma2_list, {"ser": compute_ser}, + stimulus, reference="tx", seed=4)["ser"] + +# --- sphere decoder -------------------------------------------------- +sd = Sequential([ + SymbolGenerator(M, name="tx"), + SymbolMapper(constellation), + FlatMIMOChannel(H_stack, name="channel"), + AWGN(sigma2=1.0, name="noise"), + SphereDecoder(constellation, H=H_stack, name="detector"), +], taps=["tx"], name="SD") +curves["SD"] = monte_carlo( + sd, "noise.sigma2", sigma2_list, {"ser": compute_ser}, + stimulus, reference="tx", seed=5)["ser"] # --- results: table and figure --------------------------------------- +print() print_data({"x": snr_dB_list, "curves": curves}, xlabel="snr_dB", ylabel="SER") ax = plot_error_rate(snr_dB_list, curves, ylabel="SER", - title=f"{N_r}x{N_t} MIMO, {constellation.order}-" - f"{constellation.family}, " - f"{n_channels} channel draws per point") + title=f"{N_r}x{N_t} MIMO, {M}-PSK, " + f"{K} channel draws per point") ax.set_ylim(1e-4, 1) plt.tight_layout() plt.savefig(f"{img_dir}/monte_carlo_mimo_fig3.png") # --- what the pruning is worth, in seconds --------------------------- # 16-QAM on four streams: 65 536 candidates per vector, which is where -# the two detectors stop costing the same. +# the exhaustive search and the sphere decoder stop costing the same. print("\n16-QAM, 4x4: same decision, what it costs") big_constellation = Constellation("QAM", 16) big_H = rayleigh_channel(4, 4, seed=2) -big_decoder = SphereDecoder(big_constellation, H=big_H, name="detector") -big_detectors = { - "ML": MaximumLikelihoodDetector(big_constellation, H=big_H, name="detector"), - "SD": big_decoder, -} -# --- metrics, pre-allocated ------------------------------------------ -elapsed = {} -for name in big_detectors: - elapsed[name] = np.zeros(len(snr_dB_list)) -# --- simulation loop ------------------------------------------------- -for name, detector in big_detectors.items(): - big_chain = Sequential([ - SymbolGenerator(big_constellation.order, name="tx"), - SymbolMapper(big_constellation), - FlatMIMOChannel(big_H, name="channel"), - AWGN(sigma2=1.0, name="noise"), - detector, - ], taps=["tx"], name=f"4x4 MIMO, {name}") - for index, snr_dB in enumerate(snr_dB_list): - big_chain.seed(4) - big_chain.set_params(noise__sigma2=4 * 10 ** (-snr_dB / 10)) - start = time.perf_counter() - detected = big_chain((4, 400)) - elapsed[name][index] = (time.perf_counter() - start) * 1e3 - errors = compute_ser(big_chain.tap("tx"), detected) - nodes = (f"{big_decoder.nodes_:7.1f} nodes" if name == "SD" - else f"{16 ** 4:7d} nodes") - print(f" {name:2s} {snr_dB:2d} dB {elapsed[name][index]:8.1f} ms " - f"{nodes} SER {errors:.4f}") +ml_16 = Sequential([ + SymbolGenerator(big_constellation.order, name="tx"), + SymbolMapper(big_constellation), + FlatMIMOChannel(big_H, name="channel"), + AWGN(sigma2=1.0, name="noise"), + MaximumLikelihoodDetector(big_constellation, H=big_H, + name="detector"), +], taps=["tx"], name="4x4 MIMO, ML") + +sd_16 = Sequential([ + SymbolGenerator(big_constellation.order, name="tx"), + SymbolMapper(big_constellation), + FlatMIMOChannel(big_H, name="channel"), + AWGN(sigma2=1.0, name="noise"), + SphereDecoder(big_constellation, H=big_H, name="detector"), +], taps=["tx"], name="4x4 MIMO, SD") + +elapsed = {"ML": np.zeros(len(snr_dB_list)), + "SD": np.zeros(len(snr_dB_list))} +# --- simulation loop ------------------------------------------------- +# A wall time is measured per point, so the sweep stays a loop -- one +# loop, over the SNR, with both chains timed explicitly in the body. +for index, snr_dB in enumerate(snr_dB_list): + sigma2_point = 4 * 10 ** (-snr_dB / 10) + + ml_16.seed(4) + ml_16.set_params(noise__sigma2=sigma2_point) + start = time.perf_counter() + detected = ml_16((4, 400)) + elapsed["ML"][index] = (time.perf_counter() - start) * 1e3 + errors = compute_ser(ml_16.tap("tx"), detected) + print(f" ML {snr_dB:2d} dB {elapsed['ML'][index]:8.1f} ms " + f"{16 ** 4:7d} nodes SER {errors:.4f}") + + sd_16.seed(4) + sd_16.set_params(noise__sigma2=sigma2_point) + start = time.perf_counter() + detected = sd_16((4, 400)) + elapsed["SD"][index] = (time.perf_counter() - start) * 1e3 + errors = compute_ser(sd_16.tap("tx"), detected) + print(f" SD {snr_dB:2d} dB {elapsed['SD'][index]:8.1f} ms " + f"{sd_16['detector'].nodes_:7.1f} nodes SER {errors:.4f}") + +# --- results: figure ------------------------------------------------- # A runtime is not an error rate, so it does not go through # plot_error_rate and there is no style kind for it either: what this # figure needs is a log ordinate and a grid, which is what it says. @@ -181,10 +221,3 @@ def average_ser(name, chain, snr_dB, seed=0): ax.legend() plt.tight_layout() plt.savefig(f"{img_dir}/monte_carlo_mimo_fig4.png") - -mermaid_dir = "../../docs/tutorials/mermaid/" -for diagram_name, diagram_chain in [("mimo_zf", chains["ZF"])]: - with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: - stream.write(diagram_chain.to_mermaid()) - -plt.show() diff --git a/examples/ofdm/monte_carlo_ofdm_papr.py b/examples/ofdm/monte_carlo_ofdm_papr.py index 7f3d36e..17b0845 100644 --- a/examples/ofdm/monte_carlo_ofdm_papr.py +++ b/examples/ofdm/monte_carlo_ofdm_papr.py @@ -26,28 +26,22 @@ constellation = Constellation("PSK", 4) -def get_transmitter(n_sub, oversampling=os, name="ofdm"): - """The OFDM transmitter, up to the IFFT. - - Oversampling is zero padding: the data occupy ``n_sub`` of the - ``oversampling * n_sub`` bins, so the IFFT interpolates between the - Nyquist samples and the peaks stop hiding between them. - """ - carrier_type = np.zeros(oversampling * n_sub) - carrier_type[:n_sub] = 1 - return Sequential([ - SymbolGenerator(constellation.order, name="tx"), - SymbolMapper(constellation), - Serial2Parallel(n_sub, name="s2p"), - CarrierAllocator(carrier_type=carrier_type, name="carrier_allocator"), - IFFTProcessor(), - ], name=name) - - # --- what the transmitter puts out ----------------------------------- -chain = get_transmitter(N_sc) -chain.seed(0) -blocks = chain(4 * N_sc) # four OFDM symbols +# One transmitter, up to the IFFT. Oversampling is zero padding: the +# data occupy N_sc of the os * N_sc bins, so the IFFT interpolates +# between the Nyquist samples and the peaks stop hiding between them. +carrier_1024 = np.zeros(os * N_sc) +carrier_1024[:N_sc] = 1 +ofdm_1024 = Sequential([ + SymbolGenerator(constellation.order, name="tx"), + SymbolMapper(constellation), + Serial2Parallel(N_sc, name="s2p"), + CarrierAllocator(carrier_type=carrier_1024, name="carrier_allocator"), + IFFTProcessor(), +], name="ofdm 1024") + +ofdm_1024.seed(0) +blocks = ofdm_1024(4 * N_sc) # four OFDM symbols signal = np.ravel(blocks) # C-order flatten of the (T, F) blocks power = np.abs(signal) ** 2 / np.mean(np.abs(signal) ** 2) @@ -66,7 +60,7 @@ def get_transmitter(n_sub, oversampling=os, name="ofdm"): # Each sample is a sum of N_sc independent terms, so the central limit # theorem applies: the samples are complex Gaussian, and their power is # exponential. Nothing about OFDM is needed beyond that. -long_run = np.ravel(chain(200 * N_sc)) +long_run = np.ravel(ofdm_1024(200 * N_sc)) normalized = long_run / np.std(long_run) gain = np.abs(long_run) ** 2 / np.mean(np.abs(long_run) ** 2) @@ -104,21 +98,35 @@ def get_transmitter(n_sub, oversampling=os, name="ofdm"): f"{compute_papr(blocks, unit='dB', axis=-1, reduction='mean'):.2f} dB") # --- the distribution ------------------------------------------------- -# The CCDF is estimated over many symbols, in batches, because 20 000 -# OFDM symbols at 4096 samples do not have to exist at the same time. +# Two subcarrier counts, so two chains -- the 1024 one already exists. +# The batch needs no loop: Serial2Parallel stacks the OFDM symbols on +# the leading axis, so 20 000 symbols are one call and compute_papr +# reads one value per row. +n_symbols = 20000 threshold_dB = np.arange(4, 14, 0.1) -n_batches, batch = 20, 1000 -measured, reference = {}, {} -for n_sub in (256, 1024): - transmitter = get_transmitter(n_sub) - transmitter.seed(1) - values = np.concatenate([ - compute_papr(transmitter(batch * n_sub), unit="dB", axis=-1) - for _ in range(n_batches)]) - sorted_dB, ccdf = compute_ccdf(values) - measured[f"$N_{{sc}}$ = {n_sub}"] = (sorted_dB, ccdf) - reference[f"$N_{{sc}}$ = {n_sub}"] = compute_papr_ccdf_theo( - threshold_dB, n_sub, oversampling=os, unit="dB") + +carrier_256 = np.zeros(os * 256) +carrier_256[:256] = 1 +ofdm_256 = Sequential([ + SymbolGenerator(constellation.order, name="tx"), + SymbolMapper(constellation), + Serial2Parallel(256, name="s2p"), + CarrierAllocator(carrier_type=carrier_256, name="carrier_allocator"), + IFFTProcessor(), +], name="ofdm 256") + +ofdm_256.seed(1) +papr_256 = compute_papr(ofdm_256(n_symbols * 256), unit="dB", axis=-1) +ofdm_1024.seed(1) +papr_1024 = compute_papr(ofdm_1024(n_symbols * N_sc), unit="dB", axis=-1) + +measured = {"$N_{sc}$ = 256": compute_ccdf(papr_256), + "$N_{sc}$ = 1024": compute_ccdf(papr_1024)} +reference = { + "$N_{sc}$ = 256": compute_papr_ccdf_theo(threshold_dB, 256, + oversampling=os, unit="dB"), + "$N_{sc}$ = 1024": compute_papr_ccdf_theo(threshold_dB, 1024, + oversampling=os, unit="dB")} # Markers are placed on a logarithmic grid of the ordinate: spacing them # evenly in index would crowd the top of the curve and leave the tail -- @@ -149,18 +157,3 @@ def get_transmitter(n_sub, oversampling=os, name="ofdm"): print(f"N_sc = {n_sub:4d}: PAPR exceeded once in a thousand symbols above " f"{solved:.2f} dB") -# --- the two models, against the measurement ------------------------- -# "effective" replaces the sample count by a fitted multiple of N_sc; -# "level_crossing" replaces the fitted constant by a term that grows -# with the threshold, at the price of describing the *continuous* -# waveform, whose peak a sampled one can only underestimate. -print("\nCCDF models against the measurement") -print("N_sc level threshold effective level crossing") -for name, (sorted_dB, ccdf) in measured.items(): - n_sub = int(name.split("= ")[1].rstrip("$")) - for level in (1e-2, 1e-3): - index = int(np.clip(ccdf.size * (1 - level), 0, ccdf.size - 1)) - at = sorted_dB[index] - print(f"{n_sub:5d} {ccdf[index]:6.0e} {at:8.2f} dB " - f"{compute_papr_ccdf_theo(at, n_sub, oversampling=os, unit='dB'):11.1e}" - f"{compute_papr_ccdf_theo(at, n_sub, unit='dB', method='level_crossing'):16.1e}") diff --git a/examples/ofdm/one_shot_ofdm.py b/examples/ofdm/one_shot_ofdm.py index 81c0c58..348cd40 100644 --- a/examples/ofdm/one_shot_ofdm.py +++ b/examples/ofdm/one_shot_ofdm.py @@ -109,12 +109,14 @@ snr_list = np.arange(6, 22, 2) n_trials = 2 -# --- simulation loop ------------------------------------------------- +# --- the sweep: no simulation loop ----------------------------------- measured = {} -for name, chain in (("single carrier", sc_chain), ("OFDM", ofdm_chain)): - measured[name] = monte_carlo(chain, "data_rx.snr_dB", snr_list, - {"ser": compute_ser}, (n_trials, N), - reference="data_tx", seed=1)["ser"] +measured["single carrier"] = monte_carlo( + sc_chain, "data_rx.snr_dB", snr_list, {"ser": compute_ser}, + (n_trials, N), reference="data_tx", seed=1)["ser"] +measured["OFDM"] = monte_carlo( + ofdm_chain, "data_rx.snr_dB", snr_list, {"ser": compute_ser}, + (n_trials, N), reference="data_tx", seed=1)["ser"] # --- results: table and figure --------------------------------------- @@ -128,20 +130,19 @@ # --- what it costs --------------------------------------------------- # A chain records the wall time of its last pass in `elapsed_`, so the -# run is also the measurement. +# run is also the measurement. One loop, over the swept lengths; the +# two chains are timed explicitly inside it. lengths = np.array([128, 256, 512, 1024]) +runtime = {"single carrier": np.zeros(len(lengths)), + "OFDM": np.zeros(len(lengths))} -# --- metrics, pre-allocated ------------------------------------------ -runtime = {} -for name in ("single carrier", "OFDM"): - runtime[name] = np.zeros(len(lengths)) - -# --- simulation loop ------------------------------------------------- for index, length in enumerate(lengths): - for name, chain in (("single carrier", sc_chain), ("OFDM", ofdm_chain)): - chain.seed(1) - chain(int(length)) - runtime[name][index] = 1e3 * chain.elapsed_ + sc_chain.seed(1) + sc_chain(int(length)) + runtime["single carrier"][index] = 1e3 * sc_chain.elapsed_ + ofdm_chain.seed(1) + ofdm_chain(int(length)) + runtime["OFDM"][index] = 1e3 * ofdm_chain.elapsed_ # --- results: table and figure --------------------------------------- diff --git a/examples/optical/gn_model.py b/examples/optical/gn_model.py index 1ca5c27..c0d70e8 100644 --- a/examples/optical/gn_model.py +++ b/examples/optical/gn_model.py @@ -201,9 +201,3 @@ def measure(chain, powers_W, seed=3): noiseless.set_params(fibre__noise_scaling=0.0) value = measure(noiseless, [1e-3])[0] print(f"{label:12s} {value:10.2f} dB {value - nli_only_dB:+10.2f} dB") - -mermaid_dir = "../../docs/tutorials/mermaid/" -with open(f"{mermaid_dir}/gn_model.mmd", "w") as stream: - stream.write(get_chain().to_mermaid()) - -plt.show() diff --git a/examples/simple/channel_coding.py b/examples/simple/channel_coding.py index 274d1ab..6675770 100644 --- a/examples/simple/channel_coding.py +++ b/examples/simple/channel_coding.py @@ -1,9 +1,7 @@ -import time - import numpy as np import matplotlib.pyplot as plt -from comnumpy import monte_carlo +from comnumpy import monte_carlo, print_data from comnumpy.core import Sequential from comnumpy.core.channels import AWGN from comnumpy.core.generators import SymbolGenerator @@ -42,57 +40,63 @@ def snr_dB(ebn0, code_rate): spectrum = distance_spectrum((0o133, 0o171), n_terms=4) print(f"free distance d_free = {spectrum.d_free}") -def row(label, values): - """One line of the distance spectrum, six terms wide.""" - line = label - for value in values[:6]: - line += f"{value:7d}" - return line - - -print(row("d ", spectrum.distances)) -print(row("a_d ", spectrum.a_d)) -print(row("beta_d ", spectrum.beta_d)) - - -def get_uncoded_chain(): - return Sequential([ - SymbolGenerator(2, name="tx"), - SymbolMapper(BPSK), - AWGN(snr_dB=0.0, name="noise"), - SymbolDemapper(BPSK), - ], taps=["tx"], name="uncoded BPSK") - - -def get_coded_chain(soft): - return Sequential([ - SymbolGenerator(2, name="tx"), - ConvolutionalEncoder((0o133, 0o171)), - SymbolMapper(BPSK), - AWGN(snr_dB=0.0, name="noise"), - SymbolDemapper(BPSK, soft=soft), - ViterbiDecoder((0o133, 0o171), soft=soft), - ], taps=["tx"], name=f"K=7 convolutional, {'soft' if soft else 'hard'}") - +print_data({"x": spectrum.distances[:6], + "curves": {"a_d": spectrum.a_d[:6], + "beta_d": spectrum.beta_d[:6]}}, + xlabel="d", ylabel="distance spectrum") + + +# --- three chains, three sweeps -------------------------------------- +# One Sequential per scheme, written out; monte_carlo moves the noise. +uncoded = Sequential([ + SymbolGenerator(2, name="tx"), + SymbolMapper(BPSK), + AWGN(snr_dB=0.0, name="noise"), + SymbolDemapper(BPSK), +], taps=["tx"], name="uncoded BPSK") + +hard = Sequential([ + SymbolGenerator(2, name="tx"), + ConvolutionalEncoder((0o133, 0o171)), + SymbolMapper(BPSK), + AWGN(snr_dB=0.0, name="noise"), + SymbolDemapper(BPSK), + ViterbiDecoder((0o133, 0o171)), +], taps=["tx"], name="K=7 convolutional, hard") + +soft = Sequential([ + SymbolGenerator(2, name="tx"), + ConvolutionalEncoder((0o133, 0o171)), + SymbolMapper(BPSK), + AWGN(snr_dB=0.0, name="noise"), + SymbolDemapper(BPSK, soft=True), + ViterbiDecoder((0o133, 0o171), soft=True), +], taps=["tx"], name="K=7 convolutional, soft") curves = {} -for label, chain, code_rate in (("uncoded", get_uncoded_chain(), 1.0), - ("hard-decision Viterbi", get_coded_chain(False), rate), - ("soft-decision Viterbi", get_coded_chain(True), rate)): - start = time.perf_counter() - results = monte_carlo(chain, "noise.snr_dB", snr_dB(ebn0_dB, code_rate), - {"ber": compute_ser}, n_bits, reference="tx", seed=4) - curves[label] = results["ber"] - line = f"{label:24s} " - for value in results["ber"]: - line += f"{value:.2e} " - print(line + f" ({time.perf_counter() - start:.1f} s)") - +curves["uncoded"] = monte_carlo( + uncoded, "noise.snr_dB", snr_dB(ebn0_dB, 1.0), {"ber": compute_ser}, + n_bits, reference="tx", seed=4)["ber"] +curves["hard-decision Viterbi"] = monte_carlo( + hard, "noise.snr_dB", snr_dB(ebn0_dB, rate), {"ber": compute_ser}, + n_bits, reference="tx", seed=4)["ber"] +curves["soft-decision Viterbi"] = monte_carlo( + soft, "noise.snr_dB", snr_dB(ebn0_dB, rate), {"ber": compute_ser}, + n_bits, reference="tx", seed=4)["ber"] + +# --- results: table --------------------------------------------------- +print() +print_data({"x": ebn0_dB, "curves": curves}, + xlabel="Eb/N0 [dB]", ylabel="BER") + +# --- the union bound, against the soft measurement ------------------- bound = union_bound_ber(distance_spectrum((0o133, 0o171), n_terms=8), ebn0_dB) -line = "union bound " -for value in bound: - line += f"{value:.2e} " -print(line) +print() +print_data({"x": ebn0_dB, + "curves": {"soft-decision Viterbi": curves[ + "soft-decision Viterbi"], + "union bound": bound}}, + xlabel="Eb/N0 [dB]", ylabel="BER") ax = plot_error_rate(ebn0_dB, curves, theory={"soft-decision Viterbi": bound}, x_theory=ebn0_dB, xlabel="Eb/N0 [dB]", ylabel="BER", @@ -107,25 +111,32 @@ def get_coded_chain(soft): f"information bits, rate = {ldpc_encoder.rate:.3f}, " f"column weight {H.sum(axis=0)[0]}, row weight {H.sum(axis=1)[0]}") +# One chain; the iteration count is a parameter, so the second sweep +# reconfigures the decoder with set_params instead of rebuilding. n_frames = 40 +ldpc = Sequential([ + SymbolGenerator(2, name="tx"), + LDPCEncoder(H), + SymbolMapper(BPSK), + AWGN(snr_dB=0.0, name="noise"), + SymbolDemapper(BPSK, soft=True), + LDPCDecoder(H, n_iter=5, name="decoder"), +], taps=["tx"], name="LDPC") + ldpc_curves = {} -for n_iter in (5, 25): - link = Sequential([ - SymbolGenerator(2, name="tx"), - LDPCEncoder(H), - SymbolMapper(BPSK), - AWGN(snr_dB=0.0, name="noise"), - SymbolDemapper(BPSK, soft=True), - LDPCDecoder(H, n_iter=n_iter), - ], taps=["tx"], name=f"LDPC, {n_iter} iterations") - results = monte_carlo(link, "noise.snr_dB", snr_dB(ebn0_dB, ldpc_encoder.rate), - {"ber": compute_ser}, (n_frames, ldpc_encoder.k), - reference="tx", seed=6) - ldpc_curves[f"LDPC (2040, {ldpc_encoder.k}), {n_iter} iterations"] = results["ber"] - line = f"LDPC {n_iter:2d} iterations " - for value in results["ber"]: - line += f"{value:.2e} " - print(line) +ldpc_curves[f"LDPC (2040, {ldpc_encoder.k}), 5 iterations"] = monte_carlo( + ldpc, "noise.snr_dB", snr_dB(ebn0_dB, ldpc_encoder.rate), + {"ber": compute_ser}, (n_frames, ldpc_encoder.k), + reference="tx", seed=6)["ber"] +ldpc.set_params(decoder__n_iter=25) +ldpc_curves[f"LDPC (2040, {ldpc_encoder.k}), 25 iterations"] = monte_carlo( + ldpc, "noise.snr_dB", snr_dB(ebn0_dB, ldpc_encoder.rate), + {"ber": compute_ser}, (n_frames, ldpc_encoder.k), + reference="tx", seed=6)["ber"] + +print() +print_data({"x": ebn0_dB, "curves": ldpc_curves}, + xlabel="Eb/N0 [dB]", ylabel="BER") comparison = {"uncoded": curves["uncoded"], "soft-decision Viterbi": curves["soft-decision Viterbi"]} @@ -136,9 +147,4 @@ def get_coded_chain(soft): plt.tight_layout() plt.savefig(f"{img_dir}/channel_coding_fig2.png") -mermaid_dir = "../../docs/tutorials/mermaid/" -for diagram_name, diagram_chain in [("channel_coding", get_coded_chain(True))]: - with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: - stream.write(diagram_chain.to_mermaid()) -plt.show() diff --git a/examples/simple/monte_carlo_awgn.py b/examples/simple/monte_carlo_awgn.py index 13445c2..89e8e1f 100644 --- a/examples/simple/monte_carlo_awgn.py +++ b/examples/simple/monte_carlo_awgn.py @@ -22,24 +22,15 @@ snr_dB_list = np.arange(0, 22) -def get_chain(constellation): - """Symbols in, decided symbols out, over an AWGN channel. - - The last section of this page runs one of these per constellation - order, so the constellation is an argument rather than a global. The - transmitter is named because a sweep has to compare against what it - produced, and the channel is named because a sweep has to reconfigure - it. - """ - return Sequential([ - SymbolGenerator(constellation.order, name="tx"), - SymbolMapper(constellation), - AWGN(snr_dB=0, name="awgn_channel"), - SymbolDemapper(constellation), - ], taps=["tx"]) - - -chain = get_chain(constellation) +# One chain, transmitter to decision. Every block that a sweep must +# reconfigure is named: the channel for the SNR, the generator, mapper +# and demapper for the constellation the last section changes. +chain = Sequential([ + SymbolGenerator(constellation.order, name="tx"), + SymbolMapper(constellation, name="mapper"), + AWGN(snr_dB=0, name="awgn_channel"), + SymbolDemapper(constellation, name="demapper"), +], taps=["tx"]) # A Monte Carlo simulation is three parts, always in this order: the # metric pre-allocated, the loop that fills it, the display at the end. @@ -110,7 +101,11 @@ def qfunc(x): other = Constellation("QAM", int(order)) bits = int(bits_per_symbol[index]) snr_dB_values = ebn0_to_snr_dB(ebn0_dB_list, bits_per_symbol=bits) - collected = monte_carlo(get_chain(other), "awgn_channel.snr_dB", + # the same chain, reconfigured (D50): set_params re-runs the + # coercions a bare assignment would skip + chain.set_params(tx__M=other.order, mapper__alphabet=other, + demapper__alphabet=other) + collected = monte_carlo(chain, "awgn_channel.snr_dB", snr_dB_values, {"ber": partial(compute_ber, width=bits)}, N_compare, reference="tx", seed=1) @@ -126,10 +121,3 @@ def qfunc(x): ax.set_ylim(1e-6, 1) plt.savefig(f"{img_dir}/monte_carlo_awgn_orders.png") -# The chain diagrams this tutorial shows are exported from the chains -# themselves (D33c), so the picture cannot drift from the code -- the -# smoke test compares what a run writes with what the page displays. -mermaid_dir = "../../docs/tutorials/mermaid/" -for diagram_name, diagram_chain in [("awgn_chain", chain)]: - with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: - stream.write(diagram_chain.to_mermaid()) diff --git a/examples/simple/one_shot_cfo_iq_compensation.py b/examples/simple/one_shot_cfo_iq_compensation.py index eea7c85..7f7f36f 100644 --- a/examples/simple/one_shot_cfo_iq_compensation.py +++ b/examples/simple/one_shot_cfo_iq_compensation.py @@ -55,7 +55,7 @@ plot_iq(chain.tap(tap), title=name) # show evolution of the angular frequency estimate -w0_history = chain["cfo_comp"].history +w0_history = chain["cfo_comp"].history_ plt.figure() plt.plot(w0_history) plt.xlabel("number of iteration") diff --git a/examples/simple/probabilistic_shaping.py b/examples/simple/probabilistic_shaping.py index 4e78f4e..b6de16f 100644 --- a/examples/simple/probabilistic_shaping.py +++ b/examples/simple/probabilistic_shaping.py @@ -20,7 +20,6 @@ style.use() img_dir = "../../docs/tutorials/img/" -mermaid_dir = "../../docs/tutorials/mermaid/" snr_dB = np.arange(0, 27, 2) snr = 10 ** (snr_dB / 10) @@ -158,8 +157,6 @@ def shaped_rate(lam, value): AmplitudeMapper(amplitudes, name="signs"), ], name="probabilistic amplitude shaping") chain.seed(42) -with open(f"{mermaid_dir}/shaping_pas.mmd", "w") as stream: - stream.write(chain.to_mermaid()) # Run it and count. The signs are equiprobable, so the law carried by the # signed 8-PAM is P(+/- a_i) = P_A(a_i) / 2, which is Maxwell-Boltzmann diff --git a/examples/simple/profiling_awgn_ofdm.py b/examples/simple/profiling_awgn_ofdm.py index acafa31..e1bb0b7 100644 --- a/examples/simple/profiling_awgn_ofdm.py +++ b/examples/simple/profiling_awgn_ofdm.py @@ -60,5 +60,3 @@ for diagram_name, diagram_chain in [("profiling_chain", chain)]: with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: stream.write(diagram_chain.to_mermaid()) - -plt.show() diff --git a/src/comnumpy/core/channels.py b/src/comnumpy/core/channels.py index b31c4a0..562d594 100644 --- a/src/comnumpy/core/channels.py +++ b/src/comnumpy/core/channels.py @@ -32,7 +32,11 @@ class AWGN(Processor): For real-valued inputs, a real Gaussian noise of variance :math:`\sigma^2` is applied instead of a circular complex one. - Axes: *element-wise* -- applied pointwise, shape-agnostic. + Axes: *element-wise* -- applied pointwise, shape-agnostic. One + caveat under a batch (D51): with ``snr_dB=`` the signal power is + measured over the **whole** input, batch rows included, so rows of + unequal power share one variance rather than each receiving the + requested SNR; ``sigma2=`` is batch-neutral. Parameters ---------- @@ -236,8 +240,12 @@ class TappedDelayLineChannel(Processor): .. math:: - a_0[n] = \sqrt{\frac{K}{K+1}} \, e^{j \phi} - + \sqrt{\frac{1}{K+1}} \, g[n] + a_0[n] = \sqrt{\gamma_0}\left(\sqrt{\frac{K}{K+1}} \, + e^{j \phi} + \sqrt{\frac{1}{K+1}} \, g[n]\right) + + where :math:`g[n]` is the unit-power Rayleigh process of the first + path, so that :math:`\mathbb{E}[|a_0[n]|^2] = \gamma_0` still + holds. The output keeps the input length: the delayed copies are zero-padded at the start, i.e. the channel is causal and the transient is kept. diff --git a/src/comnumpy/core/compensators.py b/src/comnumpy/core/compensators.py index f7f9527..e987df3 100644 --- a/src/comnumpy/core/compensators.py +++ b/src/comnumpy/core/compensators.py @@ -237,7 +237,9 @@ class Normalizer(Amplifier): Axes: *element-wise* -- the gain :math:`\alpha` is a single scalar measured over the whole array in ``prepare()``, then applied - pointwise. + pointwise. On a batched input the rows are therefore **pooled into + one estimate** -- normalize per trial by looping when each row must + reach the target on its own. Unlike ``Amplifier``, the gain is *not* a parameter here: it is measured, so it is exposed as ``gain_`` (decision D23) and the @@ -409,11 +411,11 @@ class BlindIQCompensator(Processor): Attributes ---------- - alpha_ : complex + alpha_ : complex or np.ndarray Estimated coefficient :math:`\alpha` applied to the in-phase component (data-dependent, hence the trailing underscore, decision D23; initialized to 1). - beta_ : complex + beta_ : complex or np.ndarray Estimated coefficient :math:`\beta` applied to the quadrature component (data-dependent, hence the trailing underscore, decision D23; initialized to 0). @@ -527,12 +529,15 @@ class BlindCFOCompensator(Processor): Newton step uses the exact first and second derivatives of the periodogram with respect to :math:`\omega`. - Axes: *declared axis* -- accepts ``(N,)`` or ``(..., P, N)``; the - correction :math:`e^{-j\widehat{\omega}_0 n}` uses the sample index - :math:`n` along the last axis. The estimand is **shared** (D49): a - frequency offset is one laser beating against one local oscillator, - so the periodogram is accumulated over every path and a single - :math:`\widehat{\omega}_0` is applied to all of them. That is not + Axes: *declared axis* -- accepts ``(N,)`` or one polarization pair + ``(2, N)``; the correction :math:`e^{-j\widehat{\omega}_0 n}` uses + the sample index :math:`n` along the last axis. The estimand is + **shared** (D49): a frequency offset is one laser beating against + one local oscillator, so the periodogram is accumulated over the + pair and a single :math:`\widehat{\omega}_0` is applied to both + rows. A wider batch is refused with ``ShapeError`` -- one scalar + estimate cannot serve independent trials (D51); loop over them, one + compensator each. That is not only tidier than one estimate per path -- it is the better estimator, since it sees all the data. @@ -552,7 +557,7 @@ class BlindCFOCompensator(Processor): refinement. save_history : bool, optional, keyword-only If True, every intermediate value of :math:`\omega_0` is appended - to ``history``. Default is False. + to ``history_``. Default is False. method : {"grad", "newton"}, optional, keyword-only Local optimizer. Default is ``"newton"``. step_size : float, optional, keyword-only @@ -570,7 +575,7 @@ class BlindCFOCompensator(Processor): Estimated frequency offset :math:`\widehat{\omega}_0` in rad/sample (data-dependent, hence the trailing underscore, decision D23; None before the first ``fit``). - history : list + history_ : list Successive iterates of :math:`\omega_0` when ``save_history=True``. @@ -616,13 +621,13 @@ class BlindCFOCompensator(Processor): name: str = field(default="cfo_compensator", kw_only=True) # internal state (declared for slots, D40a) grid_search_array: np.ndarray = field(init=False, repr=False) - history: list[float] = field(init=False, repr=False, default_factory=list) + history_: list[float] = field(init=False, repr=False, default_factory=list) # estimated quantity (D23), declared for slots (D40a) w0_: Optional[float] = field(init=False, repr=False, default_factory=lambda: None) def __post_init__(self): self.grid_search_array = np.arange(self.grid_search_tuple[0], self.grid_search_tuple[1], self.grid_search_tuple[2]) - self.history = [] + self.history_ = [] def prepare(self, x: np.ndarray) -> None: from comnumpy.exceptions import ShapeError # local import (D36) @@ -653,7 +658,7 @@ def compute_dtft(self, x: np.ndarray, w: float) -> complex: def callback(self, intermediate_result: float) -> None: if self.save_history: - self.history.append(intermediate_result) + self.history_.append(intermediate_result) def fit(self, x: np.ndarray, w0: float) -> None: w = 4*w0 @@ -671,7 +676,7 @@ def fit(self, x: np.ndarray, w0: float) -> None: index_max = np.argmax(cost_vect) w = w_vect[index_max] - self.callback(4*w) + self.callback(float(np.real(w) / 4)) for _ in range(self.N_iter): @@ -695,7 +700,7 @@ def fit(self, x: np.ndarray, w0: float) -> None: f"{self.method!r}, expected 'grad' or 'newton'.") w = w + h - self.callback(4*w) + self.callback(float(np.real(w) / 4)) self.w0_ = float(np.real(w)/4) @@ -768,7 +773,7 @@ class BlindPhaseCompensation(Processor): Attributes ---------- - theta_ : float + theta_ : float or np.ndarray Estimated phase :math:`\widehat{\theta}` in radians (data-dependent, hence the trailing underscore, decision D23). @@ -1099,7 +1104,7 @@ class DataAidedPhaseCompensator(DataAidedMixin, Processor): Attributes ---------- - theta_ : float + theta_ : float or np.ndarray Estimated phase :math:`\widehat{\theta}` in radians (data-dependent, hence the trailing underscore, decision D23). Re-estimated at every call (per-block regime, D22). @@ -1221,7 +1226,7 @@ class DataAidedComplexGainCompensator(DataAidedMixin, Processor): Attributes ---------- - gain_ : complex + gain_ : complex or np.ndarray Estimated complex gain :math:`\widehat{g}` (data-dependent, hence the trailing underscore, decision D23). @@ -1534,7 +1539,7 @@ class DataAidedFineSynchronizer(DataAidedMixin, Processor): ``Sequential(wiring={"synchronizer.reference": "source"})`` instead of freezing an array. scale_correction : bool, optional - If True (default), the output is multiplied by the complex peak + If True (default), the output is divided by the complex peak :math:`\widehat{a}`; otherwise only the delay is corrected. up_factor : int, optional, keyword-only Interpolation factor :math:`U` applied before the correlation; @@ -1556,7 +1561,8 @@ class DataAidedFineSynchronizer(DataAidedMixin, Processor): delay_ : int Estimated delay :math:`\widehat{m}_0`, expressed in *upsampled* samples (data-dependent, hence the trailing underscore, decision - D23; divide by ``up_factor`` to get it in symbol periods). + D23; divide by ``up_factor`` to get it in input sample + periods). scale_ : complex Compensation gain :math:`1/\widehat{a}` actually applied, or 1 when ``scale_correction`` is False. @@ -1738,6 +1744,13 @@ class BlindPhaseSearchCompensator(Processor): name : str, optional, keyword-only Block name. Default ``"bps"``. + Attributes + ---------- + phase_ : np.ndarray + Estimated unwrapped phase trajectory of the last call, of the + same shape as the input -- one trajectory per leading-axis row + (data-dependent, hence the trailing underscore, D23). + References ---------- T. Pfau, S. Hoffmann, R. Noe, "Hardware-Efficient Coherent Digital diff --git a/src/comnumpy/core/filters.py b/src/comnumpy/core/filters.py index 428f542..44f7d52 100644 --- a/src/comnumpy/core/filters.py +++ b/src/comnumpy/core/filters.py @@ -84,17 +84,16 @@ class SRRCFilter(Processor): If True, normalize the taps to unit energy :math:`\sum_m h^2[m] = 1`. Default is True. scale : float, optional, keyword-only - Amplitude scaling factor applied to the output. Default is 1.0. - Only the ``"fft"`` method applies it. - method : {"lfilter", "time", "fft"}, optional, keyword-only + Amplitude scaling factor applied to the output. Default is + 1.0, applied by both methods. + method : {"lfilter", "fft"}, optional, keyword-only Filtering method: ``"lfilter"`` for the causal FIR convolution (``scipy.signal.lfilter``), ``"fft"`` for the circular frequency-domain product with the delay-compensated response - returned by :meth:`H`. Default is ``"lfilter"``. The ``"time"`` - value is declared in the type but not implemented. + returned by :meth:`H`. Default is ``"lfilter"``. axis : int, optional, keyword-only - Declared filtering axis. Default is -1; the implementation always - filters along the last axis. + Filtering axis. Default is -1. The ``"lfilter"`` method honours + it; the ``"fft"`` method always filters along the last axis. name : str, optional, keyword-only Name of the filter instance. Default is ``"SRRCFilter"``. @@ -231,8 +230,8 @@ class BWFilter(Processor): H[k] = \begin{cases} - 1 & \text{if } |f[k]| \leq w_n \\ - 0 & \text{if } |f[k]| > w_n + 1 & \text{if } |f[k]| \leq w_n / 2 \\ + 0 & \text{if } |f[k]| > w_n / 2 \end{cases} where :math:`f[k]` is the frequency of bin :math:`k` in cycles per @@ -259,8 +258,6 @@ class BWFilter(Processor): ------ ValueError If ``wn`` is outside :math:`(0, 1]`. - NotImplementedError - If the input has more than one dimension. References ---------- diff --git a/src/comnumpy/core/impairments.py b/src/comnumpy/core/impairments.py index 94afa37..55a2337 100644 --- a/src/comnumpy/core/impairments.py +++ b/src/comnumpy/core/impairments.py @@ -119,6 +119,10 @@ class Delay(Processor): Axes: *axis -1* -- samples are dropped along the last axis; leading axes are batch. + The name is historical: the block models the receiver-side *removal* + of a known delay -- it advances the signal -- which is why its model + reads :math:`y[n] = x[n + \tau]`. + Parameters ---------- tau : int diff --git a/src/comnumpy/core/mappers.py b/src/comnumpy/core/mappers.py index d38a1aa..99637df 100644 --- a/src/comnumpy/core/mappers.py +++ b/src/comnumpy/core/mappers.py @@ -114,6 +114,13 @@ class SymbolDemapper(Processor): name : str, optional, keyword-only Name of the symbol demapper instance. Default is ``"Symbol Demapper"``. + Raises + ------ + ValueError + If the alphabet size is not a power of two -- the bit table + behind the LLRs is built at construction, in hard mode too, so + the constraint applies regardless of ``soft``. + References ---------- J. G. Proakis, M. Salehi, *Digital Communications*, 5th ed., diff --git a/src/comnumpy/core/processors.py b/src/comnumpy/core/processors.py index fc60a49..d6304e2 100644 --- a/src/comnumpy/core/processors.py +++ b/src/comnumpy/core/processors.py @@ -68,15 +68,17 @@ class Upsampler(Processor): Axes: *declared axis* -- the zero stuffing runs along ``axis`` (default -1, the sample axis of the Serial layout ``(..., N)``) and - the remaining axes are left untouched. + the remaining axes are left untouched. The anti-imaging filter of + ``use_filter=True`` always acts along the **last** axis, whatever + ``axis`` says. Parameters ---------- L : int Upsampling factor :math:`L`. Must be a positive integer. phase : int, optional, keyword-only - Phase offset :math:`\tau \in \mathbb{N}` of the retained samples. - Default is 0. + Phase offset :math:`\tau` of the retained samples, with + :math:`0 \le \tau < L`. Default is 0. scale : float, optional, keyword-only Output scaling factor :math:`\alpha`. Default is 1.0. axis : int, optional, keyword-only @@ -168,7 +170,9 @@ class Downsampler(Processor): Axes: *declared axis* -- the decimation runs along ``axis`` (default -1, the sample axis of the Serial layout ``(..., N)``) and the - remaining axes are left untouched. + remaining axes are left untouched. The anti-aliasing filter of + ``use_filter=True`` always acts along the **last** axis, whatever + ``axis`` says. Parameters ---------- @@ -183,9 +187,8 @@ class Downsampler(Processor): Axis along which the downsampling is performed. Default is -1. use_filter : bool, optional, keyword-only If True, apply the anti-aliasing lowpass filter of normalized - cutoff :math:`1/L` before the decimation. Default is False. - The filter is a :class:`comnumpy.core.filters.BWFilter`, which - only accepts a 1D input. + cutoff :math:`1/L` (a :class:`comnumpy.core.filters.BWFilter`) + before the decimation. Default is False. name : str, optional, keyword-only Name of the processor instance. Default is ``"downsampler"``. @@ -261,7 +264,11 @@ class Serial2Parallel(Processor): .. math:: y[t, f] = x\!\left[t F + f\right], \qquad - 0 \le t < T, \quad 0 \le f < F, \quad T = \left\lceil \frac{N}{F} \right\rceil + 0 \le t < T, \quad 0 \le f < F, \quad + T = \left\lceil \frac{N}{F} \right\rceil + \; \text{(zero-padding)}, \qquad + T = \left\lfloor \frac{N}{F} \right\rfloor + \; \text{(truncate)} Equivalently, :math:`t = \lfloor n / F \rfloor` and :math:`f = n \bmod F`: the **serial index runs fastest on the last @@ -800,8 +807,9 @@ class SampleRemover(Processor): and :math:`y[n]` has length :math:`N - D`. It is the exact inverse of :class:`DataAdder` when the same :math:`n_0` is used. - Axes: *declared axis* -- operates on a 1D serial signal ``(N,)``; the - splice uses flat indexing and does not broadcast. + Axes: *1D only* -- operates on a serial signal ``(N,)``; the + splice uses flat indexing and a batch is refused with + ``ShapeError``. Parameters ---------- @@ -924,8 +932,9 @@ class DataAdder(Processor): :math:`N + D`. Applying :class:`SampleRemover` with the same :math:`n_0` and :math:`D` restores :math:`x[n]`. - Axes: *declared axis* -- operates on a 1D serial signal ``(N,)``; the - splice is a concatenation along the first axis. + Axes: *1D only* -- operates on a serial signal ``(N,)``; the + splice is a concatenation, and a batch is refused with + ``ShapeError``. Parameters ---------- @@ -1041,7 +1050,7 @@ def forward(self, x: np.ndarray) -> np.ndarray: if self.selector is None: return x - # transformer un tuple en slice si c'est un tuple simple + # turn a plain tuple into the slice it describes if isinstance(self.selector, tuple): self.selector = slice(*self.selector) @@ -1076,8 +1085,9 @@ class Resampler(Processor): a polyphase implementation (``scipy.signal.resample_poly``), which never computes the discarded samples. - Axes: *declared axis* -- the resampling runs along the last axis, the - sample axis of the Serial layout ``(..., N)``. + Axes: *axis -1* -- the resampling runs along the last axis, the + sample axis of the Serial layout ``(..., N)``; leading axes are + batch, every row resampled by the same factors. Parameters ---------- @@ -1121,7 +1131,9 @@ def __post_init__(self): def forward(self, x: np.ndarray) -> np.ndarray: from scipy import signal # local import (D36) - y = signal.resample_poly(x, self.up, self.down) + # axis=-1 explicitly: scipy's default is axis 0, which on a + # batched (..., N) input would resample the batch axis (D51) + y = signal.resample_poly(x, self.up, self.down, axis=-1) return y @@ -1237,10 +1249,11 @@ class BlindPhaseTracker(Processor): phases give mathematically identical costs and which one wins is arbitrary. That is a property of the criterion, not of the search. - Axes: *declared axis* -- operates on a 1D serial signal ``(N,)``; the - decision of a sample under a test phase does not depend on any other - sample, so the whole search is a table of decision errors plus a - moving sum over it, and nothing is done in order. + Axes: *axis -1* -- one independent trajectory per leading-axis row + (D49): a ``(..., P, N)`` signal gets one phase estimate per path. + The decision of a sample under a test phase does not depend on any + other sample, so the whole search is a table of decision errors plus + a moving sum over it, and nothing is done in order. Parameters ---------- @@ -1262,8 +1275,8 @@ class BlindPhaseTracker(Processor): underscore). theta_ : np.ndarray The estimated phase trajectory :math:`\hat{\theta}[n]` of the last - call, of length :math:`N` (data-derived, hence the trailing - underscore, D23). A block never draws (D25, D42): plot it from the + call, of the same shape as the input (data-derived, hence the + trailing underscore, D23). A block never draws (D25, D42): plot it from the caller, for instance with ``matplotlib.pyplot.plot(tracker.theta_)``. References diff --git a/src/comnumpy/core/shaping.py b/src/comnumpy/core/shaping.py index 2e860a3..506bcf5 100644 --- a/src/comnumpy/core/shaping.py +++ b/src/comnumpy/core/shaping.py @@ -1153,8 +1153,10 @@ class DistributionMatcher(Processor): Its output is what ``SymbolMapper`` consumes, so a shaped transmitter is a chain like any other. - Axes: *element-wise* -- the bit stream is read along the last axis, - and leading axes are carried through unchanged. + Axes: *axis -1* -- the bit stream is read in blocks of ``n_bits`` + along the last axis, whose length changes (:math:`k` bits in, + :math:`n` amplitudes out); leading axes are carried through + unchanged. Parameters ---------- @@ -1214,7 +1216,9 @@ class DistributionDematcher(Processor): identity on any sequence the matcher produced -- and an error on any other, which is the honest answer after a detection mistake. - Axes: *element-wise* -- indices are read along the last axis. + Axes: *axis -1* -- indices are read in blocks of ``length`` along + the last axis, whose length changes (:math:`n` amplitudes in, + :math:`k` bits out); leading axes are carried through unchanged. Parameters ---------- @@ -1354,13 +1358,15 @@ class AmplitudeDemapper(Processor): .. math:: - \hat{x}[n] = \arg\min_i \left| \left|y[n]\right| - a_i \right| + \hat{x}[n] = \arg\min_i + \left| \left|\Re e\, y[n]\right| - a_i \right| - Taking the magnitude first is not a shortcut, it is the maximum - likelihood decision: on a constellation symmetric about the origin, - the nearest point to :math:`y` always has the amplitude nearest to - :math:`|y|`, so deciding on :math:`|y|` and deciding on :math:`y` - give the same amplitude. The sign is dropped because in PAS it + The real part is taken first because the amplitudes live on one + quadrature lane of the constellation -- each lane of a QAM symbol + is demapped on its own. On that lane, deciding on the magnitude is + not a shortcut but the maximum likelihood decision: the ASK points + are symmetric about the origin, so the nearest point always has the + amplitude nearest to :math:`|\Re e\, y|`. The sign is dropped because in PAS it carries parity, not data -- a real receiver feeds it to the FEC decoder instead. diff --git a/src/comnumpy/fec/convolutional.py b/src/comnumpy/fec/convolutional.py index 9c96962..fe2c9fa 100644 --- a/src/comnumpy/fec/convolutional.py +++ b/src/comnumpy/fec/convolutional.py @@ -48,7 +48,7 @@ class ConvolutionalEncoder(Processor): Signal Model ------------ With input bits :math:`x[m]` and generator polynomials - :math:`g_1, \dots, g_n` of constraint length :math:`K`, the coded + :math:`g_0, \dots, g_{n-1}` of constraint length :math:`K`, the coded stream interleaves the :math:`n` parity outputs: .. math:: diff --git a/src/comnumpy/mimo/channels.py b/src/comnumpy/mimo/channels.py index 9a5d462..dc78a8d 100644 --- a/src/comnumpy/mimo/channels.py +++ b/src/comnumpy/mimo/channels.py @@ -198,7 +198,9 @@ class SelectiveMIMOChannel(BaseMIMOChannel): def forward(self, X: np.ndarray) -> np.ndarray: assert self.H is not None # validate_input rejects a None channel - validate_input(X, self.H.shape[1]) + # shape[2] is N_t of the (L, N_r, N_t) stack: validating against + # shape[1] rejected every non-square channel's correct input + validate_input(X, self.H.shape[2]) L, N_r, N_t = self.H.shape N = X.shape[1] diff --git a/src/comnumpy/mimo/coding.py b/src/comnumpy/mimo/coding.py index 835bbbd..c711546 100644 --- a/src/comnumpy/mimo/coding.py +++ b/src/comnumpy/mimo/coding.py @@ -130,7 +130,9 @@ class SpaceTimeCode: \in \mathbb{C}^{N_t \times T} and the receiver observes - :math:`\mathbf{Y} = \mathbf{H}\mathbf{G}(\mathbf{s}) + \mathbf{B}`. + :math:`\mathbf{Y} = \mathbf{H}\mathbf{G}(\mathbf{s}) + \mathbf{V}`, + with :math:`\mathbf{V}` the noise (:math:`\mathbf{B}_k` being taken + by the conjugate dispersion matrices above). The **rate** is :math:`K/T` symbols per channel use, the **diversity order** is the minimum rank of :math:`\mathbf{G}(\mathbf{s}) - \mathbf{G}(\mathbf{s}')` over @@ -777,6 +779,9 @@ class SpaceTimeDecoder(Processor): ------ ValueError If the code is not orthogonal, or if ``H`` was not set. + ShapeError + If the number of received samples is not a multiple of the + code's ``n_slots``. References ---------- @@ -820,9 +825,29 @@ def forward(self, X: np.ndarray) -> np.ndarray: f"{code.name!r} spans {code.n_slots} channel uses per " f"codeword and got {received.shape[-1]} samples, which " f"is not a multiple of it.") - matrix = code.equivalent_channel(self.H) - energy = (code.orthogonality_gain - * float(np.sum(np.abs(np.asarray(self.H)) ** 2))) + H = np.asarray(self.H) + if H.ndim > 2: + # stacked channel (D51): one equivalent matrix and one + # energy per draw -- channel k decodes frame k + if received.shape[:-2] != H.shape[:-2]: + raise ShapeError( + f"a stacked channel with leading shape " + f"{H.shape[:-2]} decodes one frame per draw and " + f"needs observations (..., n_rx, N) with the same " + f"leading shape, got {received.shape}.") + flat_H = H.reshape((-1,) + H.shape[-2:]) + matrices = [] + for k in range(flat_H.shape[0]): + matrices.append(code.equivalent_channel(flat_H[k])) + matrix = np.stack(matrices).reshape( + H.shape[:-2] + matrices[0].shape) + energy = (code.orthogonality_gain + * np.sum(np.abs(H) ** 2, + axis=(-2, -1)))[..., None, None] + else: + matrix = code.equivalent_channel(H) + energy = (code.orthogonality_gain + * float(np.sum(np.abs(H) ** 2))) n_words = received.shape[-1] // code.n_slots blocks = received.reshape(received.shape[:-1] + (n_words, code.n_slots)) diff --git a/src/comnumpy/mimo/compensators.py b/src/comnumpy/mimo/compensators.py index ce7de2a..14b435e 100644 --- a/src/comnumpy/mimo/compensators.py +++ b/src/comnumpy/mimo/compensators.py @@ -14,20 +14,24 @@ class BlindDualMIMOCompensator(Processor): Signal Model ------------ At each step :math:`n`, the two output polarizations are computed - from the estimated equalizer matrix :math:`\mathbf{H}` of size - :math:`2 \times 2(2L+1)`: + from the estimated equalizer matrix :math:`\mathbf{W}` of size + :math:`2 \times 2(2L+1)` -- written :math:`\mathbf{W}`, not + :math:`\mathbf{H}`, because everywhere else in this package + :math:`\mathbf{H}` is the channel; the attribute is ``H_`` for + historical reasons, the trailing underscore marking it as estimated + (D23): .. math:: - y_i[n] = \mathbf{h}_i^H \tilde{\mathbf{x}}[n], \qquad + y_i[n] = \mathbf{w}_i^H \tilde{\mathbf{x}}[n], \qquad \tilde{\mathbf{x}}[n] = \begin{bmatrix} x_0[n] & \cdots & x_0[n-2L] & x_1[n] & \cdots & x_1[n-2L] \end{bmatrix}^T - where :math:`\mathbf{h}_i^T` is the :math:`i`-th row of - :math:`\mathbf{H}` and :math:`x_0`, :math:`x_1` are the two received + where :math:`\mathbf{w}_i^T` is the :math:`i`-th row of + :math:`\mathbf{W}` and :math:`x_0`, :math:`x_1` are the two received polarizations. The equalizer is updated by stochastic gradient - descent, :math:`\mathbf{H} \leftarrow \mathbf{H} + \mu \mathbf{g}[n]`, + descent, :math:`\mathbf{W} \leftarrow \mathbf{W} + \mu \mathbf{g}[n]`, minimizing one of the following blind losses: * Constant Modulus Algorithm (CMA): @@ -136,8 +140,11 @@ class BlindDualMIMOCompensator(Processor): Raises ------ - ValueError - If the input is not a dual-polarization signal of shape ``(2, N)``. + ShapeError + If the input does not carry a polarization pair ``(..., 2, N)``, + or if ``H_`` carries the batch shape of a previous pass and the + new input carries a different one -- call :meth:`initialize_H` + to restart from the identity equalizer. References ---------- diff --git a/src/comnumpy/mimo/detectors.py b/src/comnumpy/mimo/detectors.py index 98349a4..9760886 100644 --- a/src/comnumpy/mimo/detectors.py +++ b/src/comnumpy/mimo/detectors.py @@ -100,7 +100,7 @@ class MaximumLikelihoodDetector(Processor): and the detector returns the integer indices of :math:`\widehat{\mathbf{x}}[n]` in the alphabet - (:math:`\widehat{\mathbf{x}}[n]` = ``alphabet[S[:, n]]``). + (:math:`\widehat{\mathbf{x}}[n]` = ``alphabet[S_[:, n]]``). .. WARNING:: @@ -123,6 +123,12 @@ class MaximumLikelihoodDetector(Processor): name : str, optional, keyword-only Name of the detector. Default is ``"ML Detector"``. + Attributes + ---------- + S_ : np.ndarray + Decisions of the last call, as alphabet indices + (data-dependent, hence the trailing underscore, D23). + Raises ------ ValueError @@ -148,8 +154,8 @@ class MaximumLikelihoodDetector(Processor): alphabet: np.ndarray H: Optional[np.ndarray] = field(default=None, kw_only=True) name: str = field(default="ML Detector", kw_only=True) - # internal state (declared for slots, D40a) - S: Optional[np.ndarray] = field(init=False, repr=False, default_factory=lambda: None) + # estimated state (D23), declared for slots (D40a) + S_: Optional[np.ndarray] = field(init=False, repr=False, default_factory=lambda: None) def __post_init__(self): # accept anything array-like, a Constellation included @@ -180,10 +186,10 @@ def forward(self, Y: np.ndarray) -> np.ndarray: H = _required_channel(self.H, type(self).__name__) if H.ndim > 2: S = _detect_stacked(self._detect, H, Y) - self.S = S + self.S_ = S return S S = self._detect(H, Y) - self.S = S + self.S_ = S return S def _detect(self, H: np.ndarray, Y: np.ndarray) -> np.ndarray: @@ -306,6 +312,8 @@ class SphereDecoder(Processor): Average number of tree nodes visited per detected vector in the last call -- data-dependent, hence the trailing underscore (D23). Compare it with :math:`|\mathcal{M}|^{N_t}`. + S_ : np.ndarray + Decisions of the last call, as alphabet indices (D23). Raises ------ @@ -360,8 +368,8 @@ class SphereDecoder(Processor): H: Optional[np.ndarray] = field(default=None, kw_only=True) name: str = field(default="Sphere Decoder", kw_only=True) # internal state (declared for slots, D40a) - S: Optional[np.ndarray] = field(init=False, repr=False, - default_factory=lambda: None) + S_: Optional[np.ndarray] = field(init=False, repr=False, + default_factory=lambda: None) nodes_: float = field(init=False, repr=False, default=0.0) # running counters behind nodes_, declared for slots (D40a) _visited: int = field(init=False, repr=False, default=0) @@ -434,12 +442,12 @@ def forward(self, Y: np.ndarray) -> np.ndarray: self._visited, self._samples = 0, 0 S = _detect_stacked(self._detect, H, Y) self.nodes_ = self._visited / max(1, self._samples) - self.S = S + self.S_ = S return S self._visited, self._samples = 0, 0 S = self._detect(H, Y) self.nodes_ = self._visited / max(1, self._samples) - self.S = S + self.S_ = S return S def _detect(self, H: np.ndarray, Y: np.ndarray) -> np.ndarray: diff --git a/src/comnumpy/ofdm/compensators.py b/src/comnumpy/ofdm/compensators.py index 49bc4b9..53c6fe9 100644 --- a/src/comnumpy/ofdm/compensators.py +++ b/src/comnumpy/ofdm/compensators.py @@ -1,5 +1,5 @@ import numpy as np -from typing import Literal, Optional +from typing import Optional from dataclasses import dataclass, field from comnumpy._backend import fft, fftshift # cupy-compatible (D3) from comnumpy.core.processors import WeightAmplifier @@ -43,13 +43,13 @@ class FrequencyDomainEqualizer(WeightAmplifier): The impulse response :math:`h[l]` of the channel to be equalized. Required (a ValueError is raised when it is missing). axis : int, optional, keyword-only - The axis along which to compute the DFT and apply the weights. - Default is -1. + The axis along which the weights are applied. Default is -1, + the only value currently supported: the channel DFT is computed + on the 1-D ``h`` along its own axis, so any other value raises + at ``prepare``. shift : bool, optional, keyword-only If True, applies an fftshift to the weights, for inputs whose zero-frequency subcarrier is centered. Default is False. - norm : {"ortho", "backward", "forward"}, optional, keyword-only - FFT normalization mode. Default is "ortho". name : str, optional, keyword-only Name of the equalizer instance. Default is ``"frequency domain equalizer"``. @@ -70,7 +70,6 @@ class FrequencyDomainEqualizer(WeightAmplifier): h: Optional[np.ndarray] = None axis: int = field(default=-1, kw_only=True) shift: bool = field(default=False, kw_only=True) - norm: Literal["ortho", "backward", "forward"] = field(default="ortho", kw_only=True) name: str = field(default="frequency domain equalizer", kw_only=True) # the base class takes its weight at construction; here it is # computed from h in prepare(), so it is init=False and starts empty diff --git a/src/comnumpy/ofdm/predistorders.py b/src/comnumpy/ofdm/predistorders.py index 6fc2123..99b3314 100644 --- a/src/comnumpy/ofdm/predistorders.py +++ b/src/comnumpy/ofdm/predistorders.py @@ -35,7 +35,7 @@ class HardClipper(Processor): For real-valued inputs, :math:`T_m \, \mathrm{sign}(x[n])` replaces the polar form. - Axes: *element-wise* -- the clipping is applied pointwise; the + Axes: *axis -1* -- the clipping is applied pointwise, but the threshold :math:`T_m` is computed from the mean power **along the last axis**, so each row of a batch is clipped against its own power -- one transmitter per trial (D51). @@ -91,7 +91,7 @@ class IctPaprReductor(Processor): ------------ The frequency-domain input is iteratively transformed to the time domain, clipped, transformed back and filtered. Starting from - :math:`X^{(0)}[k] = x[k]`, each iteration :math:`i = 1, \dots, N_{it}` + :math:`X^{(0)}[k] = X[k]`, each iteration :math:`i = 1, \dots, N_{it}` computes: .. math:: @@ -107,7 +107,7 @@ class IctPaprReductor(Processor): and where: - * :math:`x[k]` is the frequency-domain input at subcarrier :math:`k`, + * :math:`X[k]` is the frequency-domain input at subcarrier :math:`k`, * :math:`P_s` is the mean power of the current time-domain block :math:`s^{(i)}[n]`, * :math:`\mathrm{PAPR_{max}}` is the target maximum PAPR in dB, * :math:`w` is the filtering weight applied in the frequency domain. @@ -190,14 +190,14 @@ class PtsPaprReductor(Processor): Signal Model ------------ - Each frequency-domain block :math:`x[k]` of length :math:`N` is + Each frequency-domain block :math:`X[k]` of length :math:`N` is partitioned into :math:`M` disjoint adjacent sub-blocks - :math:`x_m[k]` (each keeping its :math:`N/M` contiguous subcarriers + :math:`X_m[k]` (each keeping its :math:`N/M` contiguous subcarriers and zeros elsewhere). The transmitted time-domain block is the phase-rotated combination: .. math:: - y[n] = \sum_{m=1}^{M} \hat{b}_m \, \mathrm{IDFT}\{x_m\}[n] + y[n] = \sum_{m=1}^{M} \hat{b}_m \, \mathrm{IDFT}\{X_m\}[n] where the phase factors are selected by exhaustive search over the phase alphabet :math:`\mathcal{B}` to minimize the PAPR: @@ -205,24 +205,24 @@ class PtsPaprReductor(Processor): .. math:: (\hat{b}_1, \dots, \hat{b}_M) = \arg\min_{b_m \in \mathcal{B}} - \mathrm{PAPR}\left(\sum_{m=1}^{M} b_m \, \mathrm{IDFT}\{x_m\}\right) + \mathrm{PAPR}\left(\sum_{m=1}^{M} b_m \, \mathrm{IDFT}\{X_m\}\right) and where: - * :math:`x[k]` is the frequency-domain input at subcarrier :math:`k`, + * :math:`X[k]` is the frequency-domain input at subcarrier :math:`k`, * :math:`M` is the number of sub-blocks, * :math:`\mathcal{B}` is the phase factor alphabet, * :math:`y[n]` is the time-domain output block. - Axes: *declared axis* -- expects the Block layout ``(T, F)``; one - phase combination is optimized per block and the IDFT runs along the - block content axis. + Axes: *declared axis* -- expects the 2-D Block layout ``(T, F)``, + validated in ``prepare()``; one phase combination is optimized per + block and the IDFT runs along the block content axis. Parameters ---------- - phase_alphabet : list - Phase factor alphabet :math:`\mathcal{B}` used in the PTS method - (e.g. ``[1, -1]``). + phase_alphabet : list, optional + Phase factor alphabet :math:`\mathcal{B}` used in the PTS + method. Default is ``[1, -1]``. N_sub : int, optional, keyword-only Number of sub-blocks :math:`M` the OFDM signal is divided into. Default is 16. name : str, optional, keyword-only @@ -259,6 +259,14 @@ def __post_init__(self) -> None: self.combinations = np.array( list(itertools.product(self.phase_alphabet, repeat=self.N_sub))) + def prepare(self, X: np.ndarray) -> None: + from comnumpy.exceptions import ShapeError # local import (D36) + if np.ndim(X) != 2: + raise ShapeError( + f"PtsPaprReductor expects the 2-D Block layout (T, F), " + f"got shape {np.shape(X)} -- reshape with " + f"Serial2Parallel first.") + def get_subblocks(self, X: np.ndarray) -> np.ndarray: # Adjacent partition: blocks consist of a contiguous set of subcarriers and are of equal size N = len(X) diff --git a/src/comnumpy/ofdm/processors.py b/src/comnumpy/ofdm/processors.py index da7450c..af9331b 100644 --- a/src/comnumpy/ofdm/processors.py +++ b/src/comnumpy/ofdm/processors.py @@ -424,8 +424,8 @@ class CarrierAllocator(Processor): array is taken as an FFT-order mask. pilots : np.ndarray or scalar, optional, keyword-only Pilot values :math:`p[j]`, one per PILOT subcarrier of an OFDM - symbol (or a scalar broadcast to all of them). Default is an - empty array. + symbol (or a scalar broadcast to all of them). Default is + None, treated as no pilots. name : str, optional, keyword-only Name of the processor instance. Default is ``"carrier allocator"``. diff --git a/src/comnumpy/optical/channels.py b/src/comnumpy/optical/channels.py index 9e8f81f..1f1cefb 100644 --- a/src/comnumpy/optical/channels.py +++ b/src/comnumpy/optical/channels.py @@ -143,7 +143,10 @@ class ChromaticDispersion(Processor): With ``direction=-1`` the phase shift and the attenuation are both inverted (backward propagation, used for CD compensation). - Axes: *declared axis* -- requires a full-field 1D signal (N,). + Axes: *axis -1* -- the transfer function is applied along the last + axis; leading axes are batch, every row through the same fibre. + :math:`\beta_2` is given in ps^2/km and converted internally to + s^2/km (factor :math:`10^{-24}`), so the exponent is dimensionless. Parameters ---------- @@ -228,9 +231,10 @@ class KerrNonLinearity(Processor): an optional amplitude gain. With ``direction=-1`` the phase rotation is conjugated (backward propagation, used in DBP). - Axes: *declared axis* -- requires a full-field 1D signal (N,). A - pointwise Kerr step applied to separate channels would model SPM only - (no XPM, no FWM). + Axes: *element-wise* -- the rotation is pointwise, leading axes + are batch, and the shape is not validated. The signal must be the + **full field**: a pointwise Kerr step applied to separate channels + would model SPM only (no XPM, no FWM). Parameters ---------- diff --git a/src/comnumpy/optical/dbp.py b/src/comnumpy/optical/dbp.py index e2d1511..5b98c82 100644 --- a/src/comnumpy/optical/dbp.py +++ b/src/comnumpy/optical/dbp.py @@ -51,13 +51,15 @@ class DBP(Processor): which is 1.8x faster than the same ``B`` calls at 12288 samples and 3.1x at 1024 -- the split-step loop, the per-span amplifier and the parameter precomputation are paid once instead of ``B`` times. A - batch of single-polarization fields is ``(B, 1, N)``; ``(B, N)`` - reads ``B`` as polarizations and is refused. + batch of single-polarization fields is ``(B, 1, N)``; a bare + ``(B, N)`` reads ``B`` as polarizations -- refused for + :math:`B \ge 3`, while ``(2, N)`` is read as one polarization pair + by convention. Parameters ---------- N_spans : int - Number of spans to back-propagate. + Number of spans to back-propagate. Default is 1. L_span : float, keyword-only Span length :math:`L_{\mathrm{span}}` in km. Default is 80. StPS : int, keyword-only @@ -71,24 +73,10 @@ class DBP(Processor): use_only_linear : bool, keyword-only If True, only the linear effects (CD and attenuation) are compensated, in a single exact step per span. Default is False. - c : float, keyword-only - Speed of light :math:`c` in m/s. Default is ``SPEED_OF_LIGHT``. - h : float, keyword-only - Planck constant :math:`h` in J.s (kept for symmetry with - :class:`FiberLink`; unused by the deterministic inverse). - gamma : float, keyword-only - Kerr coefficient :math:`\gamma` in rad/W/km. Default is - ``KERR_COEFFICIENT``. - lamb : float, keyword-only - Wavelength :math:`\lambda` in nm. Default is ``WAVELENGTH``. - alpha_dB : float, keyword-only - Fiber loss :math:`\alpha_{dB}` in dB/km. Default is ``FIBER_LOSS``. - cd_coefficient : float, keyword-only - Dispersion coefficient :math:`D` in ps/nm/km. Default is - ``CD_COEFFICIENT``. - nu : float, keyword-only - Optical carrier frequency :math:`\nu` in Hz (kept for symmetry - with :class:`FiberLink`; unused by the deterministic inverse). + fiber : FiberSpec, keyword-only + The physical fibre to invert (D46) -- the same object handed to + :class:`FiberLink`, so the two cannot disagree. Default is a + standard single-mode fibre, ``FiberSpec()``. step_log_factor : float, keyword-only Adjustment factor of the logarithmic step-size distribution. Default is 0.4. diff --git a/src/comnumpy/optical/devices.py b/src/comnumpy/optical/devices.py index bc3b4f5..d5b48fe 100644 --- a/src/comnumpy/optical/devices.py +++ b/src/comnumpy/optical/devices.py @@ -42,7 +42,8 @@ class Laser(Processor): Laser linewidth :math:`\Delta\nu` in Hz. Default is 0 (no phase noise). theta0 : float, optional, keyword-only - Initial phase :math:`\theta_0` in rad. Drawn at random when None. + Initial phase :math:`\theta_0` in rad. When None, drawn + uniformly on :math:`[-\pi, \pi)` from the local generator. seed : int, optional, keyword-only Local RNG seed (for :math:`\theta_0` and the phase noise). fs : float, keyword-only @@ -85,7 +86,9 @@ def forward(self, nb_samples: object) -> np.ndarray: assert isinstance(nb_samples, int) assert self.rng is not None # set in __post_init__ if not isinstance(self.theta0, (float, int)): - theta_0 = 2*np.pi*self.rng.standard_normal() - np.pi + # an optical phase with no prior is uniform on [-pi, pi); + # the Gaussian this replaced had std 2*pi and no wrapping + theta_0 = float(self.rng.uniform(-np.pi, np.pi)) else: theta_0 = self.theta0 @@ -121,8 +124,8 @@ class Optical90HybridCircuit(Processor): where :math:`E_{lo}[n]` is generated by ``laser_in`` (unit amplitude when no laser is provided). - Axes: *element-wise* -- applied pointwise; the local-oscillator - realization is drawn per sample of a 1D signal (N,). + Axes: *element-wise* -- applied pointwise; one local-oscillator + draw per array element, whatever the shape. Parameters ---------- @@ -241,8 +244,10 @@ class ErbiumDopedFiberAmplifier(Processor): :math:`N_{\mathrm{ase}}` the variance of the circular amplified spontaneous emission (ASE) noise, per state of polarization. - Axes: *element-wise* -- pointwise gain plus additive noise; the ASE - realization is drawn per sample of a 1D signal (N,). + Axes: *element-wise* -- pointwise gain plus additive noise; one + independent complex ASE draw per array element, whatever the shape + -- a polarization pair ``(..., 2, N)`` gets independent noise per + polarization. Parameters ---------- @@ -329,8 +334,8 @@ class MachZehnderModulator(Processor): amplitude when no laser is provided). Axes: *element-wise* -- applied pointwise; the normalization - :math:`m` is estimated from the whole array and the laser realization - is drawn per sample of a 1D signal (N,). + :math:`m` is estimated from the whole array (batch rows pooled) and + one laser draw is made per array element. Parameters ---------- diff --git a/src/comnumpy/optical/links.py b/src/comnumpy/optical/links.py index afe623e..3f6ca89 100644 --- a/src/comnumpy/optical/links.py +++ b/src/comnumpy/optical/links.py @@ -80,13 +80,15 @@ class FiberLink(Processor): which is 1.8x faster than the same ``B`` calls at 12288 samples and 3.1x at 1024 -- the split-step loop, the per-span amplifier and the parameter precomputation are paid once instead of ``B`` times. A - batch of single-polarization fields is ``(B, 1, N)``; ``(B, N)`` - reads ``B`` as polarizations and is refused. + batch of single-polarization fields is ``(B, 1, N)``; a bare + ``(B, N)`` reads ``B`` as polarizations -- refused for + :math:`B \ge 3`, while ``(2, N)`` is read as one polarization pair + by convention. Parameters ---------- N_spans : int - Number of spans in the fiber link. + Number of spans in the fiber link. Default is 1. L_span : float, keyword-only Span length :math:`L_{\mathrm{span}}` in km. Default is 80. StPS : int, keyword-only @@ -105,23 +107,11 @@ class FiberLink(Processor): use_only_linear : bool, keyword-only If True, only the linear effects (CD and attenuation) are simulated, in a single exact step per span. Default is False. - c : float, keyword-only - Speed of light :math:`c` in m/s. Default is ``SPEED_OF_LIGHT``. - h : float, keyword-only - Planck constant :math:`h` in J.s. Default is ``PLANCK_CONSTANT``. - gamma : float, keyword-only - Kerr coefficient :math:`\gamma` in rad/W/km. Default is - ``KERR_COEFFICIENT``. - lamb : float, keyword-only - Wavelength :math:`\lambda` in nm. Default is ``WAVELENGTH``. - alpha_dB : float, keyword-only - Fiber loss :math:`\alpha_{dB}` in dB/km. Default is ``FIBER_LOSS``. - cd_coefficient : float, keyword-only - Dispersion coefficient :math:`D` in ps/nm/km. Default is - ``CD_COEFFICIENT``. - nu : float, keyword-only - Optical carrier frequency :math:`\nu` in Hz. Default is - ``OPTICAL_CARRIER_FREQUENCY``. + fiber : FiberSpec, keyword-only + The physical fibre (D46): loss :math:`\alpha_{dB}`, Kerr + coefficient :math:`\gamma`, dispersion :math:`D` and the + wavelength the carrier frequency is derived from. Default is a + standard single-mode fibre, ``FiberSpec()``. step_log_factor : float, keyword-only Adjustment factor of the logarithmic step-size distribution. Default is 0.4. @@ -339,13 +329,15 @@ def budget(self, bandwidth_Hz: float, *, .. math:: - P_{\mathrm{ASE}} = P \, N_s \, N_{\mathrm{ASE}} \, B, + P_{\mathrm{ASE}} = n_{\mathrm{pol}} \, N_s \, + N_{\mathrm{ASE}} \, B, \qquad N_{\mathrm{ASE}} = (G - 1) h \nu \, n_{sp}, \qquad n_{sp} = \frac{\mathrm{NF}/2}{1 - 1/G} - with :math:`P` the number of polarizations. That factor is the + with :math:`n_{\mathrm{pol}}` the number of polarizations + (:math:`P` denotes power everywhere else in this module). That factor is the whole reason this method exists. :math:`N_{\mathrm{ASE}}` is a density **per polarization**, and a dual-polarization link carries two independent copies of it; quoting a channel power diff --git a/src/comnumpy/optical/utils.py b/src/comnumpy/optical/utils.py index ea62349..2d17548 100644 --- a/src/comnumpy/optical/utils.py +++ b/src/comnumpy/optical/utils.py @@ -110,11 +110,11 @@ def apply_chromatic_dispersion(x: np.ndarray, z: float, beta2: float, x: numpy array Complex signal z : float - Step length in meters (km). + Step length in km. beta2: float coefficient in ps**2/km alpha_dB: float, optional - gain in dB / km (defaut: None) + attenuation in dB/km (default: None) fs : float, optional Sampling frequency in hertz (Hz). direction : int, optional @@ -455,7 +455,7 @@ def compute_erbium_doped_fiber_N_ase(alpha_dB: float, L_span: float, .. math :: - \alpha =\alpha_{dB}/10 \log_{10}(e) + \alpha = \frac{\alpha_{dB}}{10 \log_{10} e} = \frac{\ln 10}{10}\,\alpha_{dB} References ---------- diff --git a/tests/mimo/test_stacked_channel.py b/tests/mimo/test_stacked_channel.py index d9da437..ce2eb42 100644 --- a/tests/mimo/test_stacked_channel.py +++ b/tests/mimo/test_stacked_channel.py @@ -85,6 +85,21 @@ def test_a_mismatched_batch_is_refused(self): with self.assertRaises(ShapeError): detector(np.ones((K + 1, N_R, N), dtype=complex)) + def test_the_alamouti_decoder_decodes_frame_k_against_channel_k(self): + from comnumpy.mimo.coding import (SpaceTimeDecoder, + SpaceTimeEncoder, get_code) + code = get_code("alamouti") + rng = np.random.default_rng(2) + s = (rng.standard_normal((K, 40)) + + 1j * rng.standard_normal((K, 40))) + stacked = rayleigh_channel(1, 2, seed=6, size=K) + received = FlatMIMOChannel(stacked)(SpaceTimeEncoder(code)(s)) + decoded = SpaceTimeDecoder(code, H=stacked)(received) + single = SpaceTimeDecoder(code, H=stacked[2])(received[2]) + np.testing.assert_allclose(decoded[2], single) + # the orthogonal design is lossless without noise + np.testing.assert_allclose(decoded, s) + def test_the_sphere_decoder_still_counts_its_nodes(self): decoder = SphereDecoder(CONSTELLATION, H=self.H) decoder(self.Y) diff --git a/tests/test_batch_contract.py b/tests/test_batch_contract.py index 9778f92..0dcfd41 100644 --- a/tests/test_batch_contract.py +++ b/tests/test_batch_contract.py @@ -42,6 +42,25 @@ PAIR = RNG.standard_normal((3, 2, N)) + 1j * RNG.standard_normal((3, 2, N)) +def _alamouti_pieces(): + # ONE 2-D channel broadcast over the batch: the row-equality check + # runs the same factory on the batch and on one row, so the config + # must not carry the batch itself (the stacked-H path is locked in + # tests/mimo/test_stacked_channel.py) + from comnumpy.mimo.channels import FlatMIMOChannel + from comnumpy.mimo.coding import SpaceTimeEncoder, get_code + from comnumpy.mimo.utils import rayleigh_channel + channel = rayleigh_channel(1, 2, seed=3) + encoded = SpaceTimeEncoder(get_code("alamouti"))(SERIAL) + return channel, FlatMIMOChannel(channel)(encoded) + + +ALAMOUTI_H, ALAMOUTI_RX = _alamouti_pieces() +# a rotated QPSK frame set, for the blind trackers +QPSK_POINTS = np.asarray(QPSK.alphabet) +ROTATED_QPSK = QPSK_POINTS[RNG.integers(0, 4, (3, 200))] * np.exp(1j * 0.3) + + def _fs(): return {"fs": 32e9} @@ -59,8 +78,6 @@ def _fs(): np.array([1.0, 0.4])), SERIAL), "core.compensators.BlindIQCompensator": (lambda: _cls("core.compensators", "BlindIQCompensator")(), SERIAL), - "core.compensators.Normalizer": - (lambda: _cls("core.compensators", "Normalizer")("max"), SERIAL), "core.compensators.DCCorrector": (lambda: _cls("core.compensators", "DCCorrector")(), SERIAL), "core.devices.RappAmplifier": @@ -91,6 +108,11 @@ def _fs(): (lambda: _cls("core.processors", "Downsampler")(2), SERIAL), "core.processors.Upsampler": (lambda: _cls("core.processors", "Upsampler")(2), SERIAL), + "core.processors.DelayRemover": + (lambda: _cls("core.processors", "DelayRemover")(3), SERIAL), + "core.processors.BlindPhaseTracker": + (lambda: _cls("core.processors", "BlindPhaseTracker")( + 16, QPSK_POINTS), ROTATED_QPSK), "core.processors.Serial2Parallel": (lambda: _cls("core.processors", "Serial2Parallel")(8), SERIAL), "core.processors.Parallel2Serial": @@ -101,6 +123,13 @@ def _fs(): "mimo.channels.FlatMIMOChannel": (lambda: _cls("mimo.channels", "FlatMIMOChannel")( _rayleigh(2, 2, seed=1)), PAIR), + "mimo.coding.SpaceTimeEncoder": + (lambda: _cls("mimo.coding", "SpaceTimeEncoder")( + _cls("mimo.coding", "get_code")("alamouti")), SERIAL), + "mimo.coding.SpaceTimeDecoder": + (lambda: _cls("mimo.coding", "SpaceTimeDecoder")( + _cls("mimo.coding", "get_code")("alamouti"), H=ALAMOUTI_H), + ALAMOUTI_RX), "mimo.detectors.LinearDetector": (lambda: _cls("mimo.detectors", "LinearDetector")( QPSK, H=_rayleigh(2, 2, seed=1), method="zf"), PAIR), @@ -173,23 +202,26 @@ def _fs(): EXEMPT = { # abstract bases and containers "mimo.channels.BaseMIMOChannel": "abstract base, forward raises", + "core.compensators.Normalizer": "the gain is ONE scalar measured " + "over the whole array -- a documented estimand, batch rows " + "pooled by design; the docstring says so", # generic event = the whole 1-D stream by declared design; a batch # has no meaning the docstring does not already refuse "core.compensators.DataAidedFineSynchronizer": "estimates one delay against one 1-D reference", "core.compensators.DataAidedSimpleSynchronizer": "estimates one delay against one 1-D reference", - "core.frames.Framer": "frames one 1-D stream against a FrameStructure", - "core.frames.Deframer": "deframes one 1-D stream", - "core.processors.AutoConcatenator": "concatenates along the last axis " - "with internal history: one stream per instance", - "core.processors.BlindPhaseTracker": "sequential per-sample tracker, " - "one trajectory per instance (superseded by " - "BlindPhaseSearchCompensator for pairs)", + "core.frames.Framer": "frames (..., T, N) against a " + "FrameStructure; needs one, not generically constructible here", + "core.frames.Deframer": "same FrameStructure dependency", + "core.processors.AutoConcatenator": "its masks are configured by a " + "companion block (HermitianPrefixer); not generically " + "constructible", + "core.processors.DataAdder": "adds a stored companion signal of the " "un-batched shape", "core.processors.DataExtractor": "extracts against stored 1-D markers", - "core.processors.DelayRemover": "drops a scalar delay from one stream", + "core.processors.SampleRemover": "drops N samples of one stream", "core.processors.Resampler": "scipy.signal.resample along a declared " "axis; not swept here because its output length depends on the " @@ -202,10 +234,11 @@ def _fs(): "amplitude stream", "core.shaping.AmplitudeDemapper": "shaping operates on one 1-D " "amplitude stream", - "core.shaping.DistributionMatcher": "arithmetic-coding state over one " - "1-D bit stream", - "core.shaping.DistributionDematcher": "arithmetic-coding state over " - "one 1-D bit stream", + "core.shaping.DistributionMatcher": "blocks of n_bits along the " + "last axis, leading axes carried through; not swept for want of " + "a canonical probe distribution", + "core.shaping.DistributionDematcher": "blocks of `length` along the " + "last axis, leading axes carried through; same", "core.channels.TappedDelayLineChannel": "frozen-realization channel: " "one seeded draw is the configuration (D51), sounded via " "impulse_response", @@ -218,10 +251,6 @@ def _fs(): "in tests/fec", "mimo.channels.SelectiveMIMOChannel": "builds an (L N_t, N) stacked " "convolution matrix for one frame", - "mimo.coding.SpaceTimeEncoder": "Alamouti layout (N_t, N) per code " - "block, one frame per call", - "mimo.coding.SpaceTimeDecoder": "Alamouti layout per code block, one " - "frame per call", "mimo.compensators.BlindDualMIMOCompensator": "adaptive per pair; its " "batch (one butterfly per pair) is locked in " "tests/core/test_batch_axes.py",