Skip to content

Commit 4cce684

Browse files
ehsun-shclaude
andcommitted
Add wavelength-selective filtering and an optical spectrum analyser
The signal model has carried independently sampled bands and spectral noise bins since the first commit and nothing could select between them or look at them. OpticalFilter does the first, OpticalSpectrumAnalyzer the second. A filter has two jobs and the second is the one that matters here. Selecting a channel is what it looks like: because every band knows its own centre frequency this is genuinely wavelength-selective, so one tuned between two channels attenuates both. Gating ASE is the other. An amplifier emits across four terahertz and all of it reaches the diode and beats there. Eighty times less ASE gets through a 50 GHz demux, and the link then lands at 0.89x its own OSNR limit instead of 0.19x - while the OSNR figure barely moves, since it is quoted in a fixed 12.5 GHz band and cannot see what was removed outside it. PINPhotodiode's optical_bandwidth parameter is gone. It was introduced one commit ago as an explicit stand-in for this component and its own docstring said it should move here; two ways to express one piece of hardware only drift apart. The component reproduces the parameter's numbers to two decimal places, which is how I know the move was faithful. Filter skirts are floored at a declared extinction. A super-Gaussian's are not: a third-order 50 GHz passband is exp(-2838) one channel spacing off, which is not small but exactly zero, so rejection reported as infinite and a chain of filters would accumulate no crosstalk at all. Real parts specify 30-50 dB and that floor, not the shape, is what leaks. Four bugs, all mine, three of them found by the new multi-channel example: si("insertion_loss") already returns a linear ratio, so converting again turned a declared 0 dB into 0.79x. A default of "no loss" was quietly attenuating everything. The OSA stored a density and multiplied by the resolution bandwidth, which scales a carrier along with the noise floor. It must not: a carrier narrower than the setting is already inside the window. Now it stores power per bin and integrates over a Gaussian resolution filter - Gaussian and not rectangular, because a boxcar makes every peak a flat plateau and argmax then reports it half a resolution bandwidth off centre. osnr() summed every band's power over every band's noise. Right only while the channels are equal; after a demux suppressed three of four it read the survivor 6 dB low. Now per-channel, on the strongest band. The detector took its reference frequency from bands[0]. After a demux that is a suppressed neighbour a channel spacing away, whose frequency falls outside the noise bins the filter just clipped - so the ASE density read zero, the beat terms contributed nothing, and the link looked four times better than its OSNR allowed. Silently. Now the strongest band, matching the meter. Nine sabotages, nine caught, two only after tests were added: nothing asserted the OSA's absolute ASE level (every other spectrum test compares two traces, and a constant error survives all of them), and nothing asserted the extinction floor. 552 tests, 35 components. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
1 parent f98f1fb commit 4cce684

10 files changed

Lines changed: 1191 additions & 49 deletions

File tree

README.md

Lines changed: 38 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@
2929
> is drawn as a random realisation with the right Maxwellian statistics.
3030
>
3131
> **Not implemented yet:** cross-phase modulation and four-wave mixing between channels,
32-
> wavelength-selective filtering, and the GUI. See the [roadmap](#roadmap).
32+
> and the GUI. See the [roadmap](#roadmap).
3333
>
3434
> This is not yet a useful simulator. It is a foundation with the expensive decisions made and
3535
> tested. Criticism of those decisions is worth more right now than any feature —
@@ -345,6 +345,37 @@ why a coherent front end needs no optical filter at all: it is filtered by its o
345345
bandwidth, so the ASE-ASE term that forces a direct-detection receiver to carry one is absent by
346346
construction.
347347

348+
### Wavelength selection, and what a filter is really for
349+
350+
`python examples/wdm_demux.py` puts four channels on a 100 GHz grid through four amplified spans
351+
and demultiplexes one. Because every band carries its own centre frequency, that is a real
352+
wavelength-selective operation and not a choice of array index — a filter tuned *between* two
353+
channels attenuates both.
354+
355+
The second job is the one that surprises people. An amplifier emits ASE across four terahertz and
356+
every hertz of it reaches the photodiode and beats there:
357+
358+
| link | OSNR | ASE power | Q | vs OSNR limit |
359+
| :--- | ---: | ---: | ---: | ---: |
360+
| no demultiplexer | 29.93 dB | 0.3250 mW | 7.43 | 0.19× |
361+
| 50 GHz demux | 26.92 dB | 0.0040 mW | 25.04 | **0.89×** |
362+
363+
Eighty times less ASE reaches the diode, and the demultiplexed link lands just under its own OSNR
364+
limit — which is where a real receiver sits. **The OSNR figure barely moves**: it is quoted in a
365+
fixed 12.5 GHz reference bandwidth, so it cannot see ASE removed outside that band. The cheapest
366+
improvement available to a receiver is invisible to the number everyone quotes.
367+
368+
The filter's skirts are floored at a declared `extinction`, because a super-Gaussian's are not. A
369+
third-order 50 GHz passband is `exp(-2838)` one channel spacing away — not a small number but
370+
exactly zero in double precision, which would make rejection infinite and a chain of filters
371+
accumulate no crosstalk at all. Real hardware specifies 30–50 dB and it is that floor, not the
372+
shape, that decides what leaks through a long line of them.
373+
374+
The **[OSA](src/oosim/components/filters.py)** finally makes the signal model visible: bands and
375+
noise bins rendered onto one grid, the way an instrument shows them. Its resolution bandwidth is
376+
not cosmetic — widening it raises the ASE trace decibel for decibel and leaves a carrier exactly
377+
where it is, which is the clearest demonstration of why OSNR needs a stated reference bandwidth.
378+
348379
## What this is
349380

350381
A block-diagram simulator for optical systems: drop components on a canvas, wire a link, run it,
@@ -479,7 +510,7 @@ time window, and results are reproducible.
479510
| **1 — MVP: linear link** *(essentially done)* | ✅ PRBS → NRZ → laser → MZM → fiber (α + CD) → PIN → filter → eye/Q/BER, validated end to end. **Python only, no GUI.** | ~2–3 months |
480511
| **1.5 — Nonlinear & amplified**| Adaptive-step SSFM, Kerr, EDFA with ASE, OSNR, PMD, APD | ~2 months |
481512
| **2 — Coherent transceiver**| Gray-coded M-QAM to 256, IQ modulator with bias and quadrature error, 90° hybrid, balanced detection, blind carrier phase recovery, dual polarization with a blind butterfly equaliser, root-raised-cosine shaping and matched filtering, differential quadrant encoding, receiver-side dispersion compensation over spans to 1000 km, EVM/MER, constellation diagram, validated against closed-form SER | ~3 months |
482-
| **3 — GUI & WDM** | Session server, React Flow graph editor, wavelength-selective filters and an OSA · DWDM + XPM/FWM crosstalk, 400G/800G references, CuPy back-end | ~6 months |
513+
| **3 — GUI & WDM** | ✅ Wavelength-selective filters and an OSA. Session server, React Flow graph editor · DWDM + XPM/FWM crosstalk, 400G/800G references, CuPy back-end | ~6 months |
483514
| **4 — PIC** | Waveguides, ring resonators, MMI, MZI via integration with an existing S-matrix solver; PDK import ||
484515

485516
¹ One developer, part-time. Estimates, not commitments.
@@ -522,6 +553,11 @@ Every physics block ships with a test against a closed-form result, run in CI
522553
| **Signal-ASE beat** | Q on an amplified link tracks `2√(B_ref/B_e)·OSNR/(1+√(1+4·OSNR))` to 15% ||
523554
| ASE beat, coherent | Electrical SNR converges on `2·OSNR·B_ref/R_s` as ASE dominates — 0.23 dB ||
524555
| Beat is polarization-selective | Co-polarized ASE beats; orthogonal ASE does not, on both detectors ||
556+
| Filter noise bandwidth | `B_n = B·Γ(1+1/2n)/ln2^(1/2n)`, against numerical integration; order 1 is the Gaussian ||
557+
| Filtered ASE power | Exactly density × `B_n`; a demux passes its own equivalent noise bandwidth ||
558+
| Wavelength selectivity | A filter between two channels attenuates both; rejection stops at `extinction` ||
559+
| OSA normalisation | Trace integrates back to an independent power meter; ASE reads density × RBW ||
560+
| Per-channel OSNR | Survives a demultiplexer that suppresses three channels of four ||
525561
| Matched filtering | Costs `10·log10(f_s/R_s)` to omit — the receiver integrates noise it cannot use ||
526562
| **PMD** | DGD Maxwellian: `⟨τ²⟩/⟨τ⟩² = 3π/8`, mean `∝√L`, spread `0.42·mean` ||
527563
| APD | `F(M) = kM + (2−1/M)(1−k)`; an **interior optimum gain** exists ||

examples/wdm_demux.py

Lines changed: 222 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,222 @@
1+
"""Four channels on a grid, amplified, then demultiplexed.
2+
3+
The signal model has carried independently sampled bands since the first commit
4+
precisely so that this would be possible without a sample rate no machine can
5+
afford: four carriers 100 GHz apart span 300 GHz, but each is sampled only across
6+
its own 160 GHz and never on a common grid. Put those two lasers 6 THz apart
7+
instead and nothing about the run changes.
8+
9+
Three things are worth reading off the output.
10+
11+
**The demultiplexer is wavelength-selective, not an index lookup.** Each filter
12+
is tuned to one channel and rejects its neighbours by the response at their
13+
offset — a number the model produces rather than an assumption it makes. Tune one
14+
between two channels and it attenuates both.
15+
16+
**A filter is also the ASE gate.** The amplifier emits spontaneous emission
17+
across four terahertz and every hertz of it reaches a photodiode and beats there.
18+
The final table is what that costs: the same link, the same OSNR, and a factor of
19+
three in Q depending on whether anything filtered it first.
20+
21+
**OSNR does not see it.** The reference bandwidth is fixed at 12.5 GHz, so
22+
cutting ASE outside that band changes the noise power by orders of magnitude and
23+
the OSNR figure not at all. That is worth seeing once: the number everyone quotes
24+
is blind to the single cheapest improvement available to a receiver.
25+
"""
26+
27+
from __future__ import annotations
28+
29+
import math
30+
31+
import numpy as np
32+
33+
from oosim import Graph, SimulationContext
34+
from oosim.component import Component
35+
from oosim.components import (
36+
EDFA,
37+
BERAnalyzer,
38+
Combiner,
39+
CWLaser,
40+
ElectricalFilter,
41+
Fiber,
42+
MachZehnderModulator,
43+
NRZDriver,
44+
OpticalFilter,
45+
OpticalSpectrumAnalyzer,
46+
OSNRMeter,
47+
PINPhotodiode,
48+
PowerMeter,
49+
PRBSGenerator,
50+
)
51+
from oosim.kernels import super_gaussian_response
52+
from oosim.units import C_LIGHT, wavelength_to_frequency
53+
54+
ANCHOR = 1550.0 # nm — channel 0
55+
SPACING = 100e9 # Hz
56+
CHANNELS = 4
57+
BIT_RATE = 10e9
58+
FILTER_WIDTH = 50.0 # GHz
59+
OSNR_REFERENCE = 12.5e9
60+
RECEIVER_BANDWIDTH = 7.0 # GHz
61+
62+
63+
def channel_wavelength(index: int) -> float:
64+
"""Wavelength of channel ``index`` on the grid [nm]."""
65+
return C_LIGHT / (wavelength_to_frequency(ANCHOR * 1e-9) + index * SPACING) * 1e9
66+
67+
68+
def q_from_osnr(osnr_db: float) -> float:
69+
"""Textbook Q for NRZ-OOK limited by ASE beat noise."""
70+
ratio = 10.0 ** (osnr_db / 10.0)
71+
noise_bandwidth = RECEIVER_BANDWIDTH * 1e9 * math.sqrt(math.pi / (4.0 * math.log(2.0)))
72+
return (
73+
2.0
74+
* math.sqrt(OSNR_REFERENCE / noise_bandwidth)
75+
* ratio
76+
/ (1.0 + math.sqrt(1.0 + 4.0 * ratio))
77+
)
78+
79+
80+
def build(select: int | None, *, filtered: bool = True) -> tuple[Graph, dict[str, object]]:
81+
"""A four-channel comb through one amplifier, optionally demultiplexed.
82+
83+
``select`` picks which channel the receiver is tuned to; ``None`` leaves the
84+
whole comb on the detector, which is what a link without a demultiplexer
85+
would do and is included because the difference is the point.
86+
"""
87+
ctx = SimulationContext(bit_rate=BIT_RATE, samples_per_symbol=16, sequence_length=2048, seed=17)
88+
graph = Graph(ctx)
89+
combiner = graph.add(Combiner(CHANNELS, label="mux"))
90+
91+
# Only the channel under test is modulated; the neighbours are unmodulated
92+
# carriers at the same power. That keeps the crosstalk number about the
93+
# filter's response rather than about someone else's data pattern.
94+
prbs = graph.add(PRBSGenerator(order=15.0, label="prbs"))
95+
driver = graph.add(NRZDriver(v_low=4.0, v_high=0.0, label="drv"))
96+
graph.connect(prbs["out"], driver["in"])
97+
98+
for index in range(CHANNELS):
99+
laser = graph.add(
100+
CWLaser(power=0.0, wavelength=channel_wavelength(index), label=f"ch{index}")
101+
)
102+
if index == (select or 0):
103+
modulator = graph.add(MachZehnderModulator(v_pi=4.0, label=f"mzm{index}"))
104+
graph.connect(laser, modulator["optical_in"])
105+
graph.connect(driver, modulator["electrical_in"])
106+
graph.connect(modulator, combiner[f"in{index}"])
107+
else:
108+
graph.connect(laser, combiner[f"in{index}"])
109+
110+
# Four spans, each amplified back to transparency. Without real loss to make
111+
# up there is no ASE worth speaking of and the last table has nothing to show.
112+
node: Component = combiner
113+
for span in range(4):
114+
fiber = graph.add(Fiber(length=80.0, attenuation=0.2, dispersion=0.0, label=f"f{span}"))
115+
amplifier = graph.add(EDFA(gain=16.0, noise_figure=6.0, label=f"edfa{span}"))
116+
graph.connect(node, fiber["in"])
117+
graph.connect(fiber, amplifier["in"])
118+
node = amplifier
119+
120+
if filtered and select is not None:
121+
demux = graph.add(
122+
OpticalFilter(
123+
center_wavelength=channel_wavelength(select),
124+
bandwidth=FILTER_WIDTH,
125+
order=3.0,
126+
label="demux",
127+
)
128+
)
129+
graph.connect(node, demux["in"])
130+
node = demux
131+
132+
osa = graph.add(
133+
OpticalSpectrumAnalyzer(
134+
center_wavelength=channel_wavelength(1),
135+
span=800.0,
136+
points=4096,
137+
label="osa",
138+
)
139+
)
140+
meter = graph.add(PowerMeter(label="pm"))
141+
osnr = graph.add(OSNRMeter(label="osnr"))
142+
detector = graph.add(PINPhotodiode(label="pin"))
143+
receiver = graph.add(ElectricalFilter(bandwidth=RECEIVER_BANDWIDTH, label="lpf"))
144+
analyzer = graph.add(BERAnalyzer(label="ber"))
145+
graph.connect(node, osa["in"])
146+
graph.connect(node, meter["in"])
147+
graph.connect(node, osnr["in"])
148+
graph.connect(node, detector["in"])
149+
graph.connect(detector, receiver["in"])
150+
graph.connect(receiver, analyzer["in"])
151+
graph.connect(prbs["out"], analyzer["reference"])
152+
153+
return graph, {"osa": osa, "pm": meter, "osnr": osnr, "ber": analyzer}
154+
155+
156+
def main() -> None:
157+
print(f"{CHANNELS} channels on a {SPACING / 1e9:.0f} GHz grid from {ANCHOR:.1f} nm")
158+
print(f"four amplified 80 km spans, then a {FILTER_WIDTH:.0f} GHz third-order demultiplexer\n")
159+
160+
# The shape itself, before any link is involved. At a 100 GHz spacing every
161+
# sensible order buries the neighbour far below any real component's
162+
# extinction, so what limits crosstalk on this grid is the floor and not the
163+
# skirt — which is why the floor is a parameter rather than an idealisation.
164+
print(" Passband shape: transmission against offset from centre")
165+
print(f" {'offset':>10} {'order 1':>10} {'order 3':>10} {'order 5':>10}")
166+
print(" " + "-" * 46)
167+
for offset_ghz in (0.0, 12.5, 25.0, 31.25, 37.5, 50.0):
168+
row = []
169+
for order in (1, 2 + 1, 5):
170+
t = super_gaussian_response(np.array([offset_ghz * 1e9]), FILTER_WIDTH * 1e9, order)[0]
171+
db = 10.0 * math.log10(max(t**2, 1e-30))
172+
row.append(f"{db:9.1f} " if db > -300 else " --- ")
173+
print(f" {offset_ghz:8.1f}GHz " + " ".join(row))
174+
print(f" a neighbour sits {SPACING / 1e9:.0f} GHz out, past every column above —")
175+
print(" so the extinction floor, not the shape, is what it lands on.")
176+
print()
177+
178+
print(" Demultiplexer selectivity — power in each channel after the filter")
179+
print(" tuned to ch0 ch1 ch2 ch3 worst rejection")
180+
print(" " + "-" * 74)
181+
for select in range(CHANNELS):
182+
graph, ports = build(select)
183+
reading = graph.run()[ports["pm"]] # type: ignore[index]
184+
by_frequency = sorted(reading.bands, key=lambda b: b.wavelength_nm, reverse=True)
185+
powers = [band.power_dbm for band in by_frequency]
186+
others = [p for i, p in enumerate(powers) if i != select]
187+
print(
188+
f" ch{select} "
189+
+ "".join(f"{p:9.1f} " for p in powers)
190+
+ f" {powers[select] - max(others):6.1f} dB"
191+
)
192+
193+
print("\n What the filter is worth at the receiver")
194+
print(f" {'link':>28} {'OSNR':>9} {'ASE power':>11} {'Q':>7} {'vs OSNR limit':>14}")
195+
print(" " + "-" * 78)
196+
for filtered, label in ((False, "no demultiplexer"), (True, f"{FILTER_WIDTH:.0f} GHz demux")):
197+
graph, ports = build(1, filtered=filtered)
198+
results = graph.run()
199+
osnr_db = float(results[ports["osnr"]]) # type: ignore[index]
200+
reading = results[ports["pm"]] # type: ignore[index]
201+
q = results[ports["ber"]].q_factor # type: ignore[index]
202+
print(
203+
f" {label:>28} {osnr_db:6.2f} dB {reading.noise_power_w * 1e3:8.4f} mW "
204+
f"{q:7.2f} {q / q_from_osnr(osnr_db):13.2f}x"
205+
)
206+
207+
print("\n The OSNR figure is identical either way. It is quoted in a fixed")
208+
print(" 12.5 GHz reference bandwidth, so it cannot see ASE removed outside it.")
209+
210+
graph, ports = build(1)
211+
spectrum = graph.run()[ports["osa"]] # type: ignore[index]
212+
peak_frequency, peak_power = spectrum.peak()
213+
print(
214+
f"\n OSA: {len(spectrum.frequencies)} points, peak at "
215+
f"{C_LIGHT / peak_frequency * 1e9:.3f} nm, "
216+
f"{10.0 * math.log10(peak_power * 1e3):.2f} dBm per "
217+
f"{spectrum.resolution_bandwidth / 1e9:.1f} GHz"
218+
)
219+
220+
221+
if __name__ == "__main__":
222+
main()

src/oosim/analysis.py

Lines changed: 20 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -83,20 +83,33 @@ def noise_psd_at(signal: OpticalSignal, frequency: float) -> float:
8383
def osnr(signal: OpticalSignal, *, reference_bandwidth: float = OSNR_REFERENCE_BANDWIDTH) -> float:
8484
"""Optical signal-to-noise ratio [dB], in a reference bandwidth.
8585
86-
Signal power is the total in the sampled bands; noise power is the noise
87-
PSD *at each band's own centre frequency*, integrated over the reference
88-
bandwidth. Quoting OSNR in a fixed reference bandwidth rather than the
89-
signal's own is the convention, and it is the reason a 10 G and a 100 G
90-
channel with the same OSNR do not have the same margin.
86+
**OSNR is per channel, and this reports the strongest one.** Summing the
87+
power of every band and dividing by the noise at every band centre gives the
88+
right answer only while all the channels are equal, and silently the wrong one
89+
the moment they are not: after a demultiplexer suppresses three channels of
90+
four, the numerator loses them but the denominator still counts their noise,
91+
and the surviving channel's OSNR reads 6 dB low. That is what a first version
92+
of this function did, and a four-channel example is what found it.
93+
94+
Taking the strongest band is what an instrument-based measurement does — you
95+
find the peak and read the shoulder beside it — and it leaves the
96+
single-channel and equal-channel cases exactly as they were.
97+
98+
Quoting OSNR in a fixed reference bandwidth rather than the signal's own is
99+
the convention, and it is the reason a 10 G and a 100 G channel with the same
100+
OSNR do not have the same margin.
91101
"""
92102
if reference_bandwidth <= 0.0:
93103
raise ValueError(f"reference_bandwidth must be positive, got {reference_bandwidth}")
104+
if not signal.bands:
105+
return -math.inf
94106

95-
signal_power = signal.signal_power()
107+
strongest = max(signal.bands, key=lambda band: band.average_power())
108+
signal_power = strongest.average_power()
96109
if signal_power <= 0.0:
97110
return -math.inf
98111

99-
noise_power = sum(noise_psd_at(signal, band.f0) * reference_bandwidth for band in signal.bands)
112+
noise_power = noise_psd_at(signal, strongest.f0) * reference_bandwidth
100113
if noise_power <= 0.0:
101114
return math.inf
102115
return 10.0 * math.log10(signal_power / noise_power)

src/oosim/components/__init__.py

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@
1919
from .dsp import ButterflyEqualizer, DispersionCompensator
2020
from .electrical import DCVoltage, IQDriver, NRZDriver, PRBSGenerator
2121
from .fiber import Fiber
22-
from .filters import ElectricalFilter
22+
from .filters import ElectricalFilter, OpticalFilter, OpticalSpectrumAnalyzer
2323
from .mapping import DifferentialDecoder, QAMMapper
2424
from .meters import OSNRMeter, PowerMeter
2525
from .modulators import IQModulator, MachZehnderModulator
@@ -58,6 +58,8 @@
5858
"MachZehnderModulator",
5959
"NRZDriver",
6060
"OSNRMeter",
61+
"OpticalFilter",
62+
"OpticalSpectrumAnalyzer",
6163
"PINPhotodiode",
6264
"PRBSGenerator",
6365
"PolarizationCombiner",

0 commit comments

Comments
 (0)