download

Noise Under Large-Signal Modulation

4f60d9f8327f46c2abbe3a745b379337

In this notebook, we will explore how to simulate and analyze the noise performance of an Intensity-Modulated Direct-Detection (IMDD) microwave photonic link.

While basic mathematical models often assume noise behaves as a constant, stationary background hiss, real-world large-signal modulation causes the noise to pulse alongside the signal. This is known as cyclostationary noise.

To bridge the gap between simple small-signal theory and physical reality, we will leverage PhotonForge’s abstract components to run comprehensive system-level time-domain simulations. By wiring together behavioral black-box models, such as ideal lasers, electro-optic phase shifters, and noisy photodiodes, we can easily simulate the entire link architecture and capture complex, dynamic physical effects that simple analytical formulas miss.

What You Will Learn

  • Link Construction: Setting up a Mach-Zehnder Modulator (MZM) based IMDD circuit using fundamental components.

  • Metric Extraction: Calculating signal gain (\(G_c\)), noise gain (\(G_n\)), and overall Noise Figure (NF) directly from physical time-domain waveforms.

  • Large-Signal Dynamics: Understanding and visualizing how a strong RF drive artificially inflates baseband noise beyond standard mathematical predictions.

Core PhotonForge Tools

Reference

      1. Hosseini and A. Banai, “Noise figure of microwave photonic links operating under large-signal modulation and its application to optoelectronic oscillators,” Applied Optics 2014 53 (28), 6414–6421, doi: 10.1364/AO.53.006414

[1]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import rfft
import photonforge as pf
from photonforge.live_viewer import LiveViewer
import warnings

# Ignore the specific unconnected port warning from PhotonForge
warnings.filterwarnings(
    action='ignore',
    category=RuntimeWarning,
    module='photonforge.circuit_base'
)

viewer = LiveViewer()

LiveViewer started at http://localhost:45459

Calculating Equivalent Noise Bandwidth

When dealing with noise in any physical or simulated system, the shape of your measurement filter matters. A real-world bandpass filter does not have perfectly sharp “brick-wall” edges; it has roll-off.

To accurately compare our simulated noise to theoretical formulas, we must calculate the Equivalent Noise Bandwidth (ENBW) of our digital filter. The ENBW is the width of an ideal, perfectly rectangular filter that would pass the exact same amount of white noise power as our actual filter.

The mathematical definition is:

\[\text{ENBW} = \frac{\int_0^\infty |H(f)|^2 \, df}{\max(|H(f)|^2)}\]

where \(H(f)\) is the frequency response of the filter.

[4]:

def calculate_enbw(bpf, fs):
    """
    Calculates the exact Equivalent Noise Bandwidth of a pre-existing PhotonForge filter component.
    """
    time_step = 1.0 / fs
    steps = 200000

    # Inject a discrete impulse
    ts = bpf.setup_time_stepper(time_step=time_step, carrier_frequency=193e12)
    ts.reset()

    impulse_arr = np.zeros(steps)
    impulse_arr[0] = 1.0  # The impulse

    # Map the impulse to the input port
    in_port_name = "E0@0"
    out_port_name = "E1@0"
    inputs = pf.TimeSeries({in_port_name:impulse_arr}, time_step)

    outputs = ts.step(inputs, show_progress=False)

    # Extract the impulse response and calculate the Power Transfer Function
    h = np.real(outputs[out_port_name])
    H = rfft(h)
    H_power = np.abs(H)**2

    # ENBW is the integral of the power response divided by the peak
    df = fs / steps
    enbw_hz = np.sum(H_power) * df / np.max(H_power)

    return float(enbw_hz)

Before wiring up the main optical circuit, we must establish the underlying simulation environment and define our custom electrical filters.

  • Technology & Ports: We initialize PhotonForge’s base technology and set up “virtual ports.” These act as the standardized optical and electrical interfaces connecting our components.

  • RF Bandpass Filter (BPF): Positioned at the output of our link, this filter isolates our 100 MHz signal and the surrounding 20 MHz noise band, throwing away out-of-band noise to allow for clean metric extraction.

  • Input Low-Pass Filter (LPF): This scrubs high-frequency wideband noise from the input RF drive while perfectly preserving our crucial DC bias voltage.

[5]:

# Technology + virtual port specs
pf.config.default_technology = pf.basic_technology()
opt_vir = pf.virtual_port_spec(classification="optical")
elec_vir = pf.virtual_port_spec(classification="electrical")

# Create a highly selective, lossless rectangular bandpass filter
bpf_model = pf.TwoPortModel()
bpf_model.time_stepper = pf.FilterTimeStepper(
    family="rectangular",
    shape="bp",
    f_cutoff=(f_rf - B / 2, f_rf + B / 2),
    taps=4001,
    insertion_loss=0.0,
    z0=Z0,
    verbose=False,
)
bpf = bpf_model.black_box_component(elec_vir, name="RF_BPF")

# Create an Input Low-Pass Filter to scrub high-frequency noise but preserve DC Bias
in_filter_model = pf.TwoPortModel()
in_filter_model.time_stepper = pf.FilterTimeStepper(
    family="rectangular",
    shape="lp",
    f_cutoff=f_rf + B,
    taps=4001,
    insertion_loss=0.0,
    z0=Z0
)
in_filter = in_filter_model.black_box_component(elec_vir, name="Input_LPF")

Executing the Time-Domain Simulation

With our physical parameters and circuit layout defined, we need a way to easily run the simulation across different RF power levels. This cell defines run_link_once, an execution wrapper that dynamically builds the link and runs the time-stepping engine for a single set of conditions.

Here is what happens under the hood during execution:

  1. Circuit Instantiation: It calls our build_imdd_link function to construct a fresh circuit tailored to the current RF input power and noise settings.

  2. Placing Probes (Monitors): Just like hooking up an oscilloscope in the lab, we place specific monitors in the circuit. We tap into the RF input drive and the optical signal right before it hits the photodiode. (The final RF output is automatically monitored because we exposed it as a main circuit port).

  3. Simulation Engine: We initialize PhotonForge’s time_stepper, defining the central optical carrier frequency (\(f_{opt}\)) and the frequency span to track. We then step the simulation forward in time.

  4. Data Formatting: PhotonForge natively tracks electrical signals as traveling wave amplitudes to accurately model physical power flow. Before returning the data, we convert these wave amplitudes back into standard absolute voltages using the relationship \(V = A_{wave} \sqrt{Z_0}\).

[7]:

def run_link_once(
    *,
    Pin_dBm: float,
    include_rin: bool,
    include_thermal: bool,
    seed: int,
    time_step: float,
    steps: int,
):
    Pin_W = 1e-3 * 10 ** (Pin_dBm / 10)
    f_opt = pf.C_0 / 1.55

    link, ref_indices = build_imdd_link(
        f_opt=f_opt,
        f_rf=f_rf,
        B=B,
        Pin_W=Pin_W,
        include_rin=include_rin,
        include_thermal=include_thermal,
        seed=seed,
    )

    # locate references for monitors (by instance insertion order)
    rf_ref = link.references[ref_indices["rf_drive"]]
    outdc_ref = link.references[ref_indices["out_dc"]]

    # setup and run
    freq_span = 2e9
    frequencies = np.linspace(f_opt - freq_span, f_opt + freq_span, 101)

    ts = link.setup_time_stepper(
        time_step=time_step,
        carrier_frequency=f_opt,
        time_stepper_kwargs={
            "frequencies": frequencies,
            "monitors": {
                "rf_drive": rf_ref["E1"],
                # Optical monitor just before PD (output directional coupler port feeding PD)
                "opt_pre_pd": outdc_ref["P2"],
            },
        },
    )

    out = ts.step(steps=steps, time_step=time_step)

    t = out.times

    drive_wave = np.real(out["rf_drive@0+"])
    vout_wave = np.real(out["E0@0"])  # exposed port: (bpf, E1)

    # Optical field monitor (complex envelope)
    Aopt_pre_pd = out["opt_pre_pd@0+"]

    v_in = drive_wave * np.sqrt(Z0)
    v_out = vout_wave * np.sqrt(Z0)

    return t, v_in, v_out, Aopt_pre_pd

Metric Extraction from Waveforms

This section processes our raw time-domain waveforms to compute Signal Gain (\(G_c\)), Noise Gain (\(G_n\)), and Noise Figure (NF).

Because we are working with raw signal arrays containing both the powerful RF tone and the microscopic noise fluctuations, we must carefully separate them:

  1. Tone Stripping: For each measurement point, we fit an ideal sine wave to the fundamental RF carrier (\(100\) MHz) and subtract it from the raw waveform. What remains is pure residual noise.

  2. Band-limited Noise Power:

    • For the Output Noise (:math:`N_{out}`): Because our simulated circuit already includes a highly selective physical brick-wall Bandpass Filter (BPF), all out-of-band noise is gone. We can instantly find the in-band noise power simply by taking the statistical variance of the residual waveform (np.var() / R).

    • For the Input/Modulator Noise (:math:`N_{in}` and :math:`G_n`): We use a scipy.signal.periodogram to calculate the Power Spectral Density (PSD) and numerically integrate it exactly over our \(20\) MHz observation window.

  3. Isolating :math:`G_n`: To find the noise gain strictly from the modulator (without the photodiode’s thermal/shot noise), we tap into the optical monitor just before the PD, mathematically convert its complex envelope into an ideal equivalent voltage, and measure its noise power.

[8]:

def brickwall_band_power_W(v: np.ndarray, fs: float, f_lo: float, f_hi: float, R: float) -> float:
    """Calculates noise power (W) inside a strict [f_lo, f_hi] frequency band
    by integrating a one-sided periodogram (PSD)."""
    # Remove DC offset
    v = v - np.mean(v)

    # Calculate Power Spectral Density (V^2/Hz)
    f, Pxx = signal.periodogram(
        v, fs=fs, window="boxcar", detrend="constant",
        scaling="density", return_onesided=True,
    )

    # Create a boolean mask to isolate our exact frequency band
    m = (f >= f_lo) & (f <= f_hi)
    if not np.any(m):
        return 0.0

    # Numerically integrate the PSD over the band and convert to Watts
    v2 = float(np.trapezoid(Pxx[m], f[m]))
    return v2 / R

def measure_Gc_Gn_NF_photonforge(*, t: np.ndarray, v_in: np.ndarray, v_out: np.ndarray,
                                 Aopt_pre_pd: np.ndarray, Pin_W_cmd: float,
                                 steady_state: bool = True) -> dict:
    """Extracts Gc, Gn, and NF directly from PhotonForge time-domain waveforms."""

    # Drop initial simulation transients to ensure we only measure steady-state
    drop = 10000 if steady_state else 0
    tt = t[drop:]
    vin = v_in[drop:]
    vout = v_out[drop:]
    fs = 1.0 / (t[1] - t[0])
    f_lo, f_hi = f_rf - B / 2, f_rf + B / 2

    # Remove DC bias from the waveforms
    vin_ac = vin - np.mean(vin)
    vout_ac = vout - np.mean(vout)

    # 1. Extract the signal tones
    fit_in = fit_single_tone(vin_ac, tt, f_rf)
    fit_out = fit_single_tone(vout_ac, tt, f_rf)
    Pin_W_meas = power_of_tone(fit_in.amp_peak, RL)
    Pout_sig_W = power_of_tone(fit_out.amp_peak, RL)

    # 2. Strip the fundamental RF tones to leave only the residual noise
    vin_res = vin_ac - (fit_in.a_cos * np.cos(2 * np.pi * f_rf * tt) + fit_in.b_sin * np.sin(2 * np.pi * f_rf * tt))
    vout_res = vout_ac - (fit_out.a_cos * np.cos(2 * np.pi * f_rf * tt) + fit_out.b_sin * np.sin(2 * np.pi * f_rf * tt))

    # Calculate input noise power via periodogram integration
    Nin_W_meas = brickwall_band_power_W(vin_res, fs, f_lo, f_hi, RL)

    # Since v_out already passed through our physical brick-wall BPF in the simulation,
    # we can bypass the periodogram and calculate total in-band noise power directly via variance!
    Nout_W_meas = float(np.var(vout_res)) / RL

    # Calculate Signal Gain (Gc) and true physical Noise Figure (NF)
    Gc = Pout_sig_W / (Pin_W_cmd + 1e-300)
    NF_meas = (Pin_W_meas / (Nin_W_meas + 1e-300)) / (Pout_sig_W / (Nout_W_meas + 1e-300))

    # 3. Calculate Modulator Noise Gain (Gn)
    # Convert pre-PD optical envelope to an ideal equivalent voltage
    v_mod = rho * (np.abs(Aopt_pre_pd) ** 2) * TIA_gain
    vmod = v_mod[drop:]
    vmod_ac = vmod - np.mean(vmod)

    # Strip the tone and measure the modulator's specific noise contribution
    fit_mod = fit_single_tone(vmod_ac, tt, f_rf)
    vmod_res = vmod_ac - (fit_mod.a_cos * np.cos(2 * np.pi * f_rf * tt) + fit_mod.b_sin * np.sin(2 * np.pi * f_rf * tt))

    Nout_mod_W = brickwall_band_power_W(vmod_res, fs, f_lo, f_hi, RL)
    Gn_mod_meas = Nout_mod_W / (Nin_W_meas + 1e-300)

    # Return only the essential metrics needed for plotting
    return {
        "Pin_W_meas": float(Pin_W_meas),
        "Gc": float(Gc),
        "Gn_mod_meas": float(Gn_mod_meas),
        "NF_meas": float(NF_meas),
    }

Noise and Signal Gain Measurement

With our measurement functions in place, we are ready to run the core simulation sweep. We will iterate over a range of RF input powers from \(-20\) dBm to \(+15\) dBm. This sweep allows us to observe how the link transitions from a linear, small-signal regime into nonlinear, large-signal operation.

Key Simulation Parameters:

  • Sampling Rate (:math:`f_s`): We set the sampling frequency to 40 times the RF carrier frequency. This heavy oversampling is crucial to prevent aliasing and accurately capture the high-frequency harmonics generated during large-signal modulation.

  • Duration: We run each power level for 20,000 time steps. Noise is inherently random, so we need a sufficiently long time record for our statistical periodogram to yield a stable, accurate measurement of the average noise power. Increase this number for more accurate results.

  • Isolating Noise Sources: Notice that we set include_rin=False in this specific loop. By turning off the laser’s Relative Intensity Noise, we isolate the baseline thermal and shot noise behaviors. We will run a separate sweep later with RIN enabled to see its massive cyclostationary impact!

[9]:

# Sweep RF input power from -20 dBm to +15 dBm in 5 dB increments
Pin_dBm_sweep = np.arange(-20.0, 20.0, 5.0)

steps = 20000
seed0 = 100

# Set a high sampling frequency (40x the RF tone) to capture harmonics and noise
fs = 40 * f_rf
time_step = 1 / fs

# Initialize lists to store the extracted metrics for each power level
Gc_dB = []
Gn_mod_meas_dB = []

for i, Pin_dBm in enumerate(Pin_dBm_sweep):
    # Convert commanded dBm to linear Watts for the gain calculations
    Pin_W_cmd = 1e-3 * 10**(Pin_dBm / 10)

    # Run the physical simulation for this specific power level.
    # Note: RIN is explicitly turned OFF here to establish a thermal/shot noise baseline.
    t, v_in, v_out, Aopt_pre_pd = run_link_once(
        Pin_dBm=float(Pin_dBm),
        include_rin=False,
        include_thermal=True,
        seed=seed0 + i,
        time_step=time_step,
        steps=steps,
    )

    # Pass the raw waveforms to our DSP function to extract the metrics
    m = measure_Gc_Gn_NF_photonforge(
        t=t,
        v_in=v_in,
        v_out=v_out,
        Aopt_pre_pd=Aopt_pre_pd,
        Pin_W_cmd=Pin_W_cmd,
        steady_state=True,
    )

    # Convert the linear power and gain ratios back into decibels (dB / dBm)
    # A tiny offset (1e-300) is added to prevent log10(0) errors if a value drops to zero
    Gc_dB.append(10 * np.log10(m['Gc'] + 1e-300))
    Gn_mod_meas_dB.append(10 * np.log10(m['Gn_mod_meas'] + 1e-300))

Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000
Progress: 100%
Progress: 20000/20000

[10]:

# Reference data (Gc)
pin_gc_ref = [-20.04236, -15.06066, -10.07911, -4.96967, 0.01160, 4.94935, 9.96967, 15.10764]
gc_ref = [10.01362, 10.00399, 9.96728, 9.95740, 9.86655, 9.59981, 8.84560, 6.19607]

# Reference data (Gn)
pin_gn_ref = [-20.04251, -15.01808, -9.99409, -5.01268, 0.01088, 4.99013, 9.96383, 15.09229]
gn_ref = [9.98655, 10.00390, 9.94004, 9.87626, 9.73119, 9.26132, 7.74916, 3.31285]

plt.figure()

# Current simulation plots
plt.plot(Pin_dBm_sweep, Gc_dB, 'o-', label='Gc (PhotonForge)')
plt.plot(Pin_dBm_sweep, Gn_mod_meas_dB, 'o-', label='Gn (PhotonForge)')

# Reference data plots
plt.plot(pin_gc_ref, gc_ref, 's--', color='blue', alpha=0.6, label='Gc (Reference)')
plt.plot(pin_gn_ref, gn_ref, 'd--', color='orange', alpha=0.6, label='Gn (Reference)')

plt.xlabel('Pin_cmd (dBm)')
plt.ylabel('Gain (dB)')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
../_images/examples_Large_Signal_Noise_18_0.png

The Cyclostationary Sweep (Enabling RIN)

Now for the main event. We run the exact same power sweep, but with two critical changes:

  1. ``include_rin = True``: We turn the laser’s Relative Intensity Noise back on. Because RIN scales with the square of the instantaneous optical power, the large RF swings will cause the noise to fold and multiply, creating the cyclostationary effect.

  2. ``steps = 200000``: Notice we increased the simulation duration by a factor of 10. Cyclostationary noise is highly statistical and dynamic. To get a clean, converged measurement of the total Noise Figure (NF), we must observe the system over a much longer time window.

[11]:

# Increase time steps tenfold to allow statistical noise variance to fully converge
steps = 200000

# Re-initialize arrays specifically for the total Noise Figure extraction
NF_meas_dB = []
Pin_meas_dBm = []

for i, Pin_dBm in enumerate(Pin_dBm_sweep):
    Pin_W_cmd = 1e-3 * 10 ** (Pin_dBm / 10)

    # Run the physical simulation.
    # CRITICAL CHANGE: include_rin is now True!
    t, v_in, v_out, Aopt_pre_pd = run_link_once(
        Pin_dBm=float(Pin_dBm),
        include_rin=True,
        include_thermal=True,
        seed=seed0 + i,
        time_step=time_step,
        steps=steps,
    )

    # Extract the metrics from the long-duration waveforms
    m = measure_Gc_Gn_NF_photonforge(
        t=t,
        v_in=v_in,
        v_out=v_out,
        Aopt_pre_pd=Aopt_pre_pd,
        Pin_W_cmd=Pin_W_cmd,
        steady_state=True,
    )

    # Store the measured Total Noise Figure and actual measured Input Power
    NF_meas_dB.append(10 * np.log10(m['NF_meas'] + 1e-300))
    Pin_meas_dBm.append(10 * np.log10(m['Pin_W_meas'] / 1e-3 + 1e-300))

Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000
Progress: 100%
Progress: 200000/200000

Stationary vs. Cyclostationary Noise

When comparing our time-domain simulation to the paper’s standard analytical formula, you will notice a growing divergence at high RF input powers (large-signal regime). The simulated Noise Figure (NF) climbs systematically higher than the mathematical prediction.

This mismatch is not a simulation error; rather, it highlights a physical limitation of the standard analytical formula:

  • The Stationary Assumption (Standard Formula): The paper’s baseline formula relies on a small-signal approximation:

    \[NF = 10 \log_{10} \left[ \frac{2G_n(P_{in})}{G_c(P_{in})} + \frac{1}{G_c(P_{in})} + \frac{\langle I_D \rangle^2 N_{rin} R_L + 2e\langle I_D \rangle R_L}{G_c(P_{in}) k_B T_0} \right]\]

    It assumes that Shot noise and Relative Intensity Noise (RIN) act as a constant background hiss dictated strictly by the time-averaged DC photocurrent (\(\langle I_D \rangle\)).

  • The Cyclostationary Reality: Under large-signal RF drive, the MZM optical output swings violently, causing the noise itself to pulse at the RF carrier frequency. Because RIN power is proportional to the square of the instantaneous current (\(I_D(t)^2\)), its true time-averaged power grows faster than the simple DC average predicts.

  • The Correction Factor: By applying a cyclostationary correction factor—specifically, scaling the RIN term by \(1.5 - 0.5 J_0(2m)\), where \(m\) is the modulation index—we can mathematically capture this large-signal noise inflation.

The figure below displays all three curves to demonstrate how the cyclostationary correction perfectly bridges the gap between the paper’s original math and the physical simulation.

[12]:

import scipy.special as sp
from scipy.constants import e

# 1. Standard Analytical Calculation (Stationary)
# Convert Gain from dB to Linear
gc_lin = 10 ** (np.array(gc_ref) / 10)
gn_lin = 10 ** (np.array(gn_ref) / 10)

# Calculate Theoretical Noise Figure Terms
term1 = 2 * (gn_lin / gc_lin)
term2 = 1 / gc_lin

thermal_noise_base = kB * T0
shot_noise = 2 * e * Id_avg * RL
rin_noise = (Id_avg ** 2) * RIN * RL

term3 = (rin_noise + shot_noise) / (gc_lin * thermal_noise_base)

# Total Noise Factor (Stationary)
enbw = calculate_enbw(bpf, fs)
F_linear = term1 + term2 + term3 * enbw / B
nf_calculated_dB = 10 * np.log10(F_linear)

# 2. Corrected Analytical Calculation (Cyclostationary)
# Calculate Peak RF Voltage and Modulation Index (m)
V_peak = np.sqrt(2 * RL * (1e-3 * 10**(Pin_dBm_sweep / 10)))
m_mod = np.pi * V_peak / Vpi

# Calculate Cyclostationary RIN Enhancement Factor
rin_enhancement = 1.5 - 0.5 * sp.j0(2 * m_mod)

# Apply to the Noise Terms (Shot noise remains unchanged)
rin_noise_cyclo = (Id_avg ** 2) * rin_enhancement * RIN * RL

# Total Noise Factor (Cyclostationary)
term3_cyclo = (rin_noise_cyclo + shot_noise) / (gc_lin * thermal_noise_base)
F_linear_cyclo = term1 + term2 + term3_cyclo * enbw / B
nf_calculated_cyclo_dB = 10 * np.log10(F_linear_cyclo)

# 3. Final Unified Plot
plt.figure()

# Simulated
plt.plot(Pin_dBm_sweep, NF_meas_dB, 'o-', linewidth=2, label='Simulated (PhotonForge)')

# Analytical Models
plt.plot(Pin_dBm_sweep, nf_calculated_dB, 's--', alpha=0.7, label='Analytical (Stationary / Original)')
plt.plot(Pin_dBm_sweep, nf_calculated_cyclo_dB, 'd--', color='purple', linewidth=2, label='Analytical (Cyclostationary / Corrected)')

plt.xlabel('RF Input Power, $P_{in}$ (dBm)')
plt.ylabel('Noise Figure, NF (dB)')
plt.title('Large-Signal Noise Figure: Simulation vs. Analytical Models')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
../_images/examples_Large_Signal_Noise_22_0.png