Enhance visualizations and animations for FEC and LDPC decoding examples

- Improved OFDM signal analysis with comprehensive plots for I/Q components, instantaneous power, constellation, and power distribution.
- Added detailed MIMO constraint analysis with per-antenna power and PAPR comparisons.
- Updated FEC decoder tutorial to clarify concepts and improve output formatting.
- Enhanced LDPC visualization with improved Tanner graph and added GIF animations for belief propagation iterations.
- Introduced GIF animations for successive cancellation decoding process in polar codes.

Author: selimfirat

examples/benchmarks/plot_visualization_example.py

@@ -84,32 +84,49 @@ def run_visualization_example():
     # Create comparison plot if we have multiple results
     print("\n5. Running parameter comparison...")
 
-    # Compare different modulation schemes
-    modulations = ["bpsk", "qpsk", "16qam"]
+    # Compare different modulation schemes using appropriate benchmarks
     comparison_results = []
-
-    for mod in modulations:
-        print(f"   Running {mod.upper()} simulation...")
-        mod_benchmark = get_benchmark("ber_simulation")(modulation=mod)
-        mod_result = runner.run_benchmark(mod_benchmark, snr_range=list(range(0, 16, 2)), num_bits=20000)
-        comparison_results.append(mod_result.metrics)
+    modulation_labels = []
+
+    # BPSK using BER simulation benchmark
+    print("   Running BPSK simulation...")
+    bpsk_benchmark = get_benchmark("ber_simulation")(modulation="bpsk")
+    bpsk_result = runner.run_benchmark(bpsk_benchmark, snr_range=list(range(0, 16, 2)), num_bits=20000)
+    comparison_results.append(bpsk_result.metrics)
+    modulation_labels.append("BPSK")
+
+    # 4-QAM (QPSK) using QAM benchmark
+    print("   Running QPSK simulation...")
+    qpsk_benchmark = get_benchmark("qam_ber")(constellation_size=4)
+    qpsk_result = runner.run_benchmark(qpsk_benchmark, snr_range=list(range(0, 16, 2)), num_symbols=10000)
+    comparison_results.append(qpsk_result.metrics)
+    modulation_labels.append("QPSK")
+
+    # 16-QAM using QAM benchmark
+    print("   Running 16-QAM simulation...")
+    qam16_benchmark = get_benchmark("qam_ber")(constellation_size=16)
+    qam16_result = runner.run_benchmark(qam16_benchmark, snr_range=list(range(0, 16, 2)), num_symbols=10000)
+    comparison_results.append(qam16_result.metrics)
+    modulation_labels.append("16-QAM")
 
     # Plot comparison
     print("\n6. Creating modulation comparison plot...")
     # Create individual BER plots for each modulation scheme
-    for i, (mod, result_metrics) in enumerate(zip(modulations, comparison_results)):
-        plot_name = f"ber_curve_{mod}.png"
+    for i, (mod_label, result_metrics) in enumerate(zip(modulation_labels, comparison_results)):
+        plot_name = f"ber_curve_{mod_label.lower().replace('-', '')}.png"
         visualizer.plot_ber_curve(result_metrics, save_path=str(output_dir / plot_name))
-        print(f"✓ {mod.upper()} BER curve saved to visualization_results/{plot_name}")
+        print(f"✓ {mod_label} BER curve saved to visualization_results/{plot_name}")
 
     # Create a combined comparison plot manually using matplotlib
     plt.figure(figsize=(12, 8))
-    for mod, result_metrics in zip(modulations, comparison_results):
+    for mod_label, result_metrics in zip(modulation_labels, comparison_results):
         snr_range = result_metrics.get("snr_range", [])
         if "ber_simulated" in result_metrics:
-            plt.semilogy(snr_range, result_metrics["ber_simulated"], "o-", label=f"{mod.upper()} (Simulated)", linewidth=2, markersize=6)
+            plt.semilogy(snr_range, result_metrics["ber_simulated"], "o-", label=f"{mod_label} (Simulated)", linewidth=2, markersize=6)
+        elif "ber_results" in result_metrics:
+            plt.semilogy(snr_range, result_metrics["ber_results"], "o-", label=f"{mod_label} (Simulated)", linewidth=2, markersize=6)
         if "ber_theoretical" in result_metrics:
-            plt.semilogy(snr_range, result_metrics["ber_theoretical"], "--", label=f"{mod.upper()} (Theoretical)", linewidth=2)
+            plt.semilogy(snr_range, result_metrics["ber_theoretical"], "--", label=f"{mod_label} (Theoretical)", linewidth=2)
 
     plt.xlabel("SNR (dB)", fontsize=12)
     plt.ylabel("Bit Error Rate", fontsize=12)

examples/channels/plot_binary_channels.py

@@ -307,9 +307,11 @@
 
 # Calculate capacities for each channel type
 def calculate_bsc_capacity(p):
-    """Calculate BSC capacity: C = 1 - H(p)"""
+    """Calculate BSC capacity: C = 1 - H(p) where H(p) is binary entropy"""
     if p == 0 or p == 1:
         return 1.0 if p == 0 else 0.0
+    # Binary entropy H(p) = -p*log2(p) - (1-p)*log2(1-p)
+    # BSC capacity C = 1 - H(p)
     return 1 + p * np.log2(p) + (1 - p) * np.log2(1 - p)
 
 
@@ -323,9 +325,18 @@ def calculate_z_capacity(p):
     if p == 0:
         return 1.0
     if p == 1:
+        return np.log2(2) - 1  # log2(2) = 1, so this is 0
+
+    # Z-channel capacity: C = log2(1 + (1-p)^(1-p) * p^p)
+    # This is the correct formula for Z-channel capacity
+    try:
+        term1 = (1 - p) ** (1 - p) if (1 - p) > 0 else 0
+        term2 = p**p if p > 0 else 1
+        capacity = np.log2(1 + term1 * term2)
+        return capacity
+    except (ValueError, ZeroDivisionError):
+        # Handle edge cases
         return 0.0
-    # Z-channel capacity formula
-    return 1 + (1 - p) * np.log2(1 - p) + p * np.log2(p)
 
 
 # Calculate capacities

examples/channels/plot_fading_channels.py

@@ -78,10 +78,32 @@
 # First 5 complex symbols: {qpsk_symbols[:5]}
 
 # Show the QPSK constellation diagram
-fig, ax = plt.subplots(figsize=(6, 6))
-qpsk_modulator.plot_constellation()
-ax.set_title("QPSK Constellation")
-ax.grid(True)
+fig, ax = plt.subplots(figsize=(8, 6))
+
+# Create QPSK constellation points manually
+qpsk_points = np.array([1 + 1j, -1 + 1j, -1 - 1j, 1 - 1j]) / np.sqrt(2)  # Normalized QPSK
+labels = ["00", "10", "11", "01"]
+
+# Plot constellation points
+ax.scatter(qpsk_points.real, qpsk_points.imag, s=100, c="red", marker="o", edgecolors="black", linewidth=2)
+
+# Add labels for each point
+for point, label in zip(qpsk_points, labels):
+    ax.annotate(label, (point.real, point.imag), xytext=(10, 10), textcoords="offset points", fontsize=12, fontweight="bold")
+
+# Set up the plot
+ax.set_xlabel("In-phase (I)", fontweight="bold")
+ax.set_ylabel("Quadrature (Q)", fontweight="bold")
+ax.set_title("QPSK Constellation", fontweight="bold", fontsize=14)
+ax.grid(True, alpha=0.3)
+ax.axis("equal")
+ax.set_xlim(-1.5, 1.5)
+ax.set_ylim(-1.5, 1.5)
+
+# Add axes through origin
+ax.axhline(y=0, color="k", linewidth=0.5)
+ax.axvline(x=0, color="k", linewidth=0.5)
+
 plt.tight_layout()
 plt.show()
 
@@ -182,29 +204,30 @@
 # --------------------------------------------------
 # Let's analyze how fading affects the amplitude distribution of the symbols.
 
-# Calculate amplitudes - ensure we're using real values
-perfect_amp = np.abs(perfect_complex).real
-awgn_amp = np.abs(awgn_complex).real
-fading_amp = np.abs(fading_complex).real
+# Calculate amplitudes correctly
+perfect_amp = np.abs(perfect_complex)
+awgn_amp = np.abs(awgn_complex)
+fading_amp = np.abs(fading_complex)
+
 plt.figure(figsize=(12, 5))
 
 # Histogram of amplitudes
 plt.subplot(1, 2, 1)
-plt.hist(perfect_amp, bins=30, alpha=0.3, label="Perfect Channel", density=True)
-plt.hist(awgn_amp, bins=30, alpha=0.3, label="AWGN Channel", density=True)
-plt.hist(fading_amp, bins=30, alpha=0.3, label="Fading Channel", density=True)
-plt.grid(True)
-plt.xlabel("Symbol Amplitude")
-plt.ylabel("Probability Density")
-plt.title("Symbol Amplitude Distribution")
+plt.hist(perfect_amp, bins=30, alpha=0.6, label="Perfect Channel", density=True, color="blue")
+plt.hist(awgn_amp, bins=30, alpha=0.6, label="AWGN Channel", density=True, color="green")
+plt.hist(fading_amp, bins=30, alpha=0.6, label="Fading Channel", density=True, color="red")
+plt.grid(True, alpha=0.3)
+plt.xlabel("Symbol Amplitude", fontweight="bold")
+plt.ylabel("Probability Density", fontweight="bold")
+plt.title("Symbol Amplitude Distribution", fontweight="bold")
 plt.legend()
 
 # Theoretical vs. Empirical Rayleigh Distribution
+plt.subplot(1, 2, 2)
 x = np.linspace(0, 3, 1000)
 # Rayleigh PDF: (x/σ²) * exp(-x²/(2σ²))
 # For unit variance Rayleigh, σ² = 1/2
 rayleigh_pdf = x * np.exp(-(x**2) / 2)
-plt.subplot(1, 2, 2)
 plt.hist(fading_amp, bins=30, alpha=0.5, density=True, label="Empirical (Fading Channel)")
 plt.plot(x, rayleigh_pdf, "r-", linewidth=2, label="Theoretical Rayleigh")
 plt.grid(True)
@@ -362,18 +385,20 @@ def calculate_ser(received, original_indices):
 # -------------------------------------------------------------------------
 # Let's analyze the frequency characteristics of the fading process.
 
-# Calculate PSD using FFT
+# Calculate PSD using Welch's method
+# Ensure nperseg is not larger than input length
+input_length = len(fading_magnitude)
+nperseg = min(256, input_length // 4)  # Use at most 1/4 of input length
 
-# Use Welch's method to estimate PSD
-f, psd = signal.welch(fading_magnitude, fs=1.0, nperseg=256)
+f, psd = signal.welch(fading_magnitude, fs=1.0, nperseg=nperseg)
 
 plt.figure(figsize=(10, 5))
-plt.semilogy(f, psd, linewidth=2)
-plt.grid(True)
-plt.xlabel("Normalized Frequency")
-plt.ylabel("Power Spectral Density")
-plt.title("PSD of Rayleigh Fading Process")
-plt.axvline(x=0.05, color="r", linestyle="--", label="Doppler Frequency (0.05)")
+plt.semilogy(f, psd, linewidth=2, color="blue")
+plt.grid(True, alpha=0.3)
+plt.xlabel("Normalized Frequency", fontweight="bold")
+plt.ylabel("Power Spectral Density", fontweight="bold")
+plt.title("PSD of Rayleigh Fading Process", fontweight="bold")
+plt.axvline(x=0.05, color="r", linestyle="--", linewidth=2, label="Doppler Frequency (0.05)")
 plt.legend()
 plt.tight_layout()
 plt.show()

examples/channels/plot_poisson_channel.py

@@ -375,29 +375,47 @@ def compute_local_snr(signal, noise, signal_levels):
 poisson_snr = compute_local_snr(poisson_profile["ramp_signal"], poisson_profile["ramp_noise"], signal_levels)
 awgn_snr = compute_local_snr(awgn_profile["ramp_signal"], awgn_profile["ramp_noise"], signal_levels)
 
-# Plot SNR vs signal level
-if poisson_snr and awgn_snr:
-    p_centers, p_snr_db = zip(*poisson_snr)
-    a_centers, a_snr_db = zip(*awgn_snr)
+# Create the plot
+plt.figure(figsize=(10, 6))
 
-    plt.plot(p_centers, p_snr_db, "ro-", label="Poisson Channel")
-    plt.plot(a_centers, a_snr_db, "go-", label="AWGN Channel")
+# Initialize flag to track if any data was plotted
+data_plotted = False
 
-    # Plot theoretical curves
-    x = np.linspace(0.1, 1.0, 100)
-    # For Poisson: SNR = rate * x (signal ^ 2 / signal = signal)
-    poisson_theory_snr = 10 * np.log10(rate * x)
-    # For AWGN: SNR = x^2 / (1/rate) = x^2 * rate
-    awgn_theory_snr = 10 * np.log10(x**2 * rate)
+# Plot SNR vs signal level if we have data
+if poisson_snr:
+    p_centers, p_snr_db = zip(*poisson_snr)
+    plt.plot(p_centers, p_snr_db, "ro-", linewidth=2, markersize=6, label="Poisson Channel")
+    data_plotted = True
 
-    plt.plot(x, poisson_theory_snr, "r--", alpha=0.7, label="Poisson Theory")
-    plt.plot(x, awgn_theory_snr, "g--", alpha=0.7, label="AWGN Theory")
+if awgn_snr:
+    a_centers, a_snr_db = zip(*awgn_snr)
+    plt.plot(a_centers, a_snr_db, "go-", linewidth=2, markersize=6, label="AWGN Channel")
+    data_plotted = True
 
-plt.grid(True)
-plt.title(f"Local SNR vs. Signal Level (Rate = {rate})")
-plt.xlabel("Signal Intensity")
-plt.ylabel("SNR (dB)")
-plt.legend()
+# Always plot theoretical curves for reference
+x = np.linspace(0.1, 1.0, 100)
+# For Poisson: SNR = signal (since noise variance = signal for Poisson)
+poisson_theory_snr = 10 * np.log10(rate * x)
+# For AWGN: SNR = x^2 / noise_power = x^2 * rate (for equivalent SNR)
+awgn_theory_snr = 10 * np.log10(x**2 * rate)
+
+plt.plot(x, poisson_theory_snr, "r--", alpha=0.7, linewidth=2, label="Poisson Theory")
+plt.plot(x, awgn_theory_snr, "g--", alpha=0.7, linewidth=2, label="AWGN Theory")
+data_plotted = True
+
+plt.grid(True, alpha=0.3)
+plt.title(f"Local SNR vs Signal Level (Rate = {rate})", fontweight="bold", fontsize=14)
+plt.xlabel("Signal Intensity", fontweight="bold")
+plt.ylabel("SNR (dB)", fontweight="bold")
+
+# Only call legend if we actually plotted something
+if data_plotted:
+    plt.legend()
+else:
+    plt.text(0.5, 0.5, "No data available for plotting", transform=plt.gca().transAxes, ha="center", va="center", fontsize=12)
+
+plt.tight_layout()
+plt.show()
 plt.tight_layout()
 plt.show()