Lab 2: Cognitive Radio -- Spectrum Sensing and Opportunistic Access Pipeline

Week: 3 -- Cognitive Radio
Points: 20
Time estimate: 90 min lab + 2.5 hr independent
Deliverable: lab-2-report.md + Python scripts

Objectives

1. Implement an energy detector with configurable threshold; verify detection probability and false-alarm rate against the theoretical model.
2. Build an occupancy map across M simulated channels from wideband spectrum data.
3. Implement an opportunistic channel selector that finds and uses available channels.
4. Demonstrate the cognitive cycle: sense → decide → transmit → detect-primary-return → vacate.

Part A: Energy Detector Calibration (30 min)

#!/usr/bin/env python3
"""Lab 2 Part A: Energy detector calibration."""
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

np.random.seed(42)

def energy_detector(x, N_window, tau):
    """
    Sliding window energy detector.
    x:        input signal (complex)
    N_window: sensing window length (samples)
    tau:      detection threshold
    Returns:  detection decision (True/False) for each window
    """
    n_windows = len(x) // N_window
    decisions = []
    for i in range(n_windows):
        window = x[i*N_window:(i+1)*N_window]
        energy = np.mean(np.abs(window)**2)
        decisions.append(energy > tau)
    return decisions

def theoretical_threshold(N, noise_var, P_fa):
    """Threshold for P_fa false alarm rate, sensing window N."""
    # Under H0: 2*E/noise_var ~ chi2(2N); Gaussian approx for large N
    sigma_E = noise_var / np.sqrt(N)
    return noise_var + stats.norm.ppf(1 - P_fa) * sigma_E

# Simulation parameters
noise_var = 1.0
N_window = 1000      # sensing window samples
P_fa_target = 0.01

# Compute threshold for target false alarm rate
tau = theoretical_threshold(N_window, noise_var, P_fa_target)
print(f"Threshold for P_fa={P_fa_target}: τ = {tau:.4f} (noise σ² = {noise_var})")

# Monte Carlo evaluation: measure P_fa (under H0)
n_trials = 10000
noise_only = np.random.randn(n_trials * N_window) + 1j * np.random.randn(n_trials * N_window)
noise_only *= np.sqrt(noise_var / 2)  # normalize to variance noise_var

decisions_H0 = energy_detector(noise_only, N_window, tau)
P_fa_measured = np.mean(decisions_H0)
print(f"Measured P_fa: {P_fa_measured:.4f} (target: {P_fa_target})")

# ROC curve: P_d vs P_fa at various SNR levels
snr_db_levels = [-10, -5, 0, 5, 10]
tau_range = np.linspace(noise_var * 0.5, noise_var * 2.0, 100)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for snr_db in snr_db_levels:
    snr_linear = 10**(snr_db/10)
    signal_var = noise_var * snr_linear
    
    # Generate signal + noise
    signal = np.sqrt(signal_var/2) * (np.random.randn(n_trials * N_window) +
                                       1j * np.random.randn(n_trials * N_window))
    x_signal = signal + noise_only
    
    P_fa_roc, P_d_roc = [], []
    for tau_roc in tau_range:
        d_H0 = energy_detector(noise_only, N_window, tau_roc)
        d_H1 = energy_detector(x_signal, N_window, tau_roc)
        P_fa_roc.append(np.mean(d_H0))
        P_d_roc.append(np.mean(d_H1))
    
    axes[0].plot(P_fa_roc, P_d_roc, label=f'SNR={snr_db} dB')

axes[0].set_xlabel('P_fa (False Alarm Rate)')
axes[0].set_ylabel('P_d (Detection Probability)')
axes[0].set_title('ROC Curves - Energy Detector')
axes[0].legend()
axes[0].grid(True)
axes[0].plot([0, 1], [0, 1], 'k--', alpha=0.3)

# P_d vs SNR at fixed P_fa
P_d_vs_snr = []
snr_range = np.linspace(-15, 10, 30)
tau_fixed = theoretical_threshold(N_window, noise_var, P_fa_target)
for snr_db in snr_range:
    snr_linear = 10**(snr_db/10)
    signal_var = noise_var * snr_linear
    signal = np.sqrt(signal_var/2) * (np.random.randn(n_trials * N_window) +
                                       1j * np.random.randn(n_trials * N_window))
    x_signal = signal + noise_only
    decisions = energy_detector(x_signal, N_window, tau_fixed)
    P_d_vs_snr.append(np.mean(decisions))

axes[1].plot(snr_range, P_d_vs_snr, 'b-', label=f'N={N_window}, P_fa={P_fa_target}')
axes[1].axhline(0.9, color='r', linestyle='--', label='P_d=0.9 target')
axes[1].set_xlabel('SNR (dB)')
axes[1].set_ylabel('P_d')
axes[1].set_title('Detection Probability vs. SNR')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.savefig('lab2/roc_curves.png', dpi=150)

Record: At what SNR (dB) does the energy detector achieve P_d ≥ 0.9 at P_fa = 0.01 with N=1000?

Part B: Occupancy Map Across M Channels (30 min)

Simulate a wideband spectrum with M=16 channels. Some channels have primary users (randomly assigned). Build the occupancy map using the energy detector.

#!/usr/bin/env python3
"""Lab 2 Part B: Multi-channel occupancy map."""
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)
M = 16             # number of channels
N_sensing = 500    # sensing samples per channel
fs = 16e6          # total bandwidth: 16 MHz (1 MHz per channel)
noise_var = 1.0
P_fa = 0.01
tau = theoretical_threshold(N_sensing, noise_var, P_fa)

# Ground truth occupancy (random; 40% channels occupied)
true_occupancy = np.random.rand(M) < 0.4
print(f"True occupied channels: {np.where(true_occupancy)[0].tolist()}")
print(f"True vacant channels: {np.where(~true_occupancy)[0].tolist()}")

# Simulate sensing of each channel
SNR_primary_dB = 3.0   # primary users at low SNR (challenging scenario)
snr_linear = 10**(SNR_primary_dB/10)

occupancy_map = np.zeros(M, dtype=bool)
energy_per_channel = np.zeros(M)

for ch in range(M):
    # Generate samples for this channel
    noise = np.sqrt(noise_var/2) * (np.random.randn(N_sensing) + 1j*np.random.randn(N_sensing))
    
    if true_occupancy[ch]:
        # Primary user present
        signal_var = noise_var * snr_linear
        signal = np.sqrt(signal_var/2) * (np.random.randn(N_sensing) + 1j*np.random.randn(N_sensing))
        x = signal + noise
    else:
        x = noise
    
    energy = np.mean(np.abs(x)**2)
    energy_per_channel[ch] = energy
    occupancy_map[ch] = energy > tau

# Evaluate detection performance
true_positive  = np.sum(true_occupancy & occupancy_map)
false_positive = np.sum(~true_occupancy & occupancy_map)
true_negative  = np.sum(~true_occupancy & ~occupancy_map)
false_negative = np.sum(true_occupancy & ~occupancy_map)

print(f"\nOccupancy detection results:")
print(f"  True positives (correctly detected occupied): {true_positive}")
print(f"  False positives (incorrectly flagged vacant): {false_positive}")
print(f"  True negatives (correctly identified vacant): {true_negative}")
print(f"  False negatives (missed primary users): {false_negative}")
print(f"  Detection rate: {true_positive/np.sum(true_occupancy):.2f}")
print(f"  False alarm rate: {false_positive/np.sum(~true_occupancy):.2f}")

# Visualize occupancy map
fig, ax = plt.subplots(figsize=(12, 4))
bar_colors = []
for ch in range(M):
    if true_occupancy[ch] and occupancy_map[ch]:
        bar_colors.append('red')      # correctly detected occupied
    elif not true_occupancy[ch] and not occupancy_map[ch]:
        bar_colors.append('green')    # correctly detected vacant
    elif not true_occupancy[ch] and occupancy_map[ch]:
        bar_colors.append('orange')   # false positive
    else:
        bar_colors.append('blue')     # missed detection

ax.bar(range(M), energy_per_channel, color=bar_colors)
ax.axhline(tau, color='black', linestyle='--', linewidth=2, label=f'Threshold τ={tau:.3f}')
ax.set_xlabel('Channel index')
ax.set_ylabel('Energy')
ax.set_title(f'Occupancy Map (M={M} channels, SNR={SNR_primary_dB} dB)')
ax.legend()
plt.tight_layout()
plt.savefig('lab2/occupancy_map.png', dpi=150)

Part C: Opportunistic Channel Selector and Cognitive Cycle (30 min)

#!/usr/bin/env python3
"""Lab 2 Part C: Opportunistic access and cognitive cycle demonstration."""
import numpy as np
import matplotlib.pyplot as plt

class CognitiveRadio:
    """Minimal cognitive radio implementing the sense-decide-transmit cycle."""
    
    def __init__(self, M_channels, noise_var, P_fa, N_sensing, vacate_time_slots=3):
        self.M = M_channels
        self.noise_var = noise_var
        self.tau = theoretical_threshold(N_sensing, noise_var, P_fa)
        self.N_sensing = N_sensing
        self.vacate_time_slots = vacate_time_slots
        self.current_channel = None
        self.state = 'SENSING'    # SENSING → TRANSMITTING → VACATING
        self.log = []
    
    def sense(self, channel_energies):
        """Stage 1: Sense all channels, build occupancy map."""
        occupancy = channel_energies > self.tau
        return occupancy
    
    def decide(self, occupancy):
        """Stage 2: Select best vacant channel."""
        vacant = np.where(~occupancy)[0]
        if len(vacant) == 0:
            return None
        # Simple strategy: pick the channel with lowest energy (most vacant)
        return vacant[0]
    
    def transmit(self, channel):
        """Stage 3: Transmit on selected channel."""
        self.current_channel = channel
        self.state = 'TRANSMITTING'
        self.log.append(f"TX on ch {channel}")
    
    def check_primary_return(self, channel_energy):
        """Stage 4: Check if primary returned on current channel."""
        return channel_energy > self.tau
    
    def vacate(self):
        """Vacate channel due to primary return."""
        old_ch = self.current_channel
        self.current_channel = None
        self.state = 'SENSING'
        self.log.append(f"VACATED ch {old_ch}: primary detected")

# Simulate 100 time slots
np.random.seed(42)
M = 8
N_sensing = 500
noise_var = 1.0
snr_primary_dB = 5.0
snr_linear = 10**(snr_primary_dB/10)

cr = CognitiveRadio(M, noise_var, P_fa=0.01, N_sensing=N_sensing)

# Primary user activity: channel 3 returns at slot 50
tx_log = []
channel_log = []

for slot in range(100):
    # Generate channel energies
    # Channel 2 always occupied; channel 3 occupied after slot 50
    occupancy_true = np.zeros(M, dtype=bool)
    occupancy_true[2] = True  # always occupied
    if slot >= 50:
        occupancy_true[3] = True  # primary returns at slot 50
    
    energies = np.zeros(M)
    for ch in range(M):
        noise = np.sqrt(noise_var/2) * (np.random.randn(N_sensing) + 1j*np.random.randn(N_sensing))
        if occupancy_true[ch]:
            sig_var = noise_var * snr_linear
            s = np.sqrt(sig_var/2) * (np.random.randn(N_sensing) + 1j*np.random.randn(N_sensing))
            energies[ch] = np.mean(np.abs(s + noise)**2)
        else:
            energies[ch] = np.mean(np.abs(noise)**2)
    
    if cr.state == 'SENSING':
        occ = cr.sense(energies)
        ch = cr.decide(occ)
        if ch is not None:
            cr.transmit(ch)
    elif cr.state == 'TRANSMITTING':
        if cr.current_channel is not None:
            if cr.check_primary_return(energies[cr.current_channel]):
                cr.vacate()
    
    channel_log.append(cr.current_channel if cr.current_channel is not None else -1)

print("\nCognitive Radio Event Log:")
for event in cr.log:
    print(f"  {event}")

# Plot channel assignment over time
plt.figure(figsize=(12, 4))
plt.step(range(100), channel_log, where='post')
plt.axvline(50, color='r', linestyle='--', label='Primary returns on ch 3 at slot 50')
plt.xlabel('Time slot')
plt.ylabel('CR channel in use (-1 = sensing)')
plt.title('Cognitive Radio Channel Assignment Over Time')
plt.legend()
plt.grid(True)
plt.savefig('lab2/cr_channel_log.png', dpi=150)

Record:

1. What channel does the CR initially select?
2. At slot 50 (primary returns), does the CR detect the return and vacate? Within how many slots?
3. What channel does the CR move to after vacating?

Lab Report

Create lab-2-report.md with:

1. ROC curve (P_d vs. P_fa) at 5 SNR levels; annotate the operating point for P_fa=0.01
2. Minimum SNR for P_d ≥ 0.9 at P_fa = 0.01 with N=1000 (read from curve)
3. Occupancy map figure with color-coded detection results; detection rate and false alarm rate
4. Cognitive radio event log showing initial channel selection, primary-return detection, and channel handoff
5. Analysis: "In a real cognitive radio implementation, what is the trade-off between sensing window length (N) and throughput? At N=1000 samples and fs=1 MHz, how many ms does one sensing window take, and what fraction of a 10 ms frame does it consume?"

Grading