Lab 8: FHSS Analysis and Chirp Spread Spectrum

Week: 10 -- Anti-Jamming, LPI/LPD
Points: 20
Time estimate: 90 min lab + 1 hr independent
Deliverable: lab-8-report.md + Python scripts

Objectives

1. Simulate a 16-channel FHSS signal and recover the hop sequence from a spectrogram.
2. Generate a chirp spread spectrum signal and verify the matched-filter processing gain formula.
3. Calculate an FHSS anti-jamming margin budget from a given system specification.
4. Analyze why the Week 7 SIGINT pipeline fails on LPI/LPD waveforms.

Setup

mkdir -p lab8/plots
cd lab8

Required packages: numpy, scipy, matplotlib. All are available in the standard RF-301 environment.

Part A: FHSS Simulation and Hop Sequence Recovery (7 points)

A.1 Simulate the FHSS Signal

#!/usr/bin/env python3
"""Lab 8 Part A: FHSS simulation and hop sequence recovery."""
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import spectrogram

def simulate_fhss(n_channels=16, hop_rate=1600, fs=4e6, duration=0.1, seed=42):
    """
    Simulate a 16-channel FHSS signal.
    
    Parameters:
        n_channels: number of hop channels
        hop_rate: hops per second
        fs: sample rate (Hz)
        duration: signal duration (seconds)
        seed: RNG seed for reproducibility
    
    Returns:
        iq: complex IQ samples (complex64)
        hop_sequence: array of channel indices used
        hop_freqs: channel center frequencies (Hz)
    """
    f_min = -0.8e6
    f_max = 0.8e6
    hop_freqs = np.linspace(f_min, f_max, n_channels)
    
    samples_per_hop = int(fs / hop_rate)
    total_samples = int(fs * duration)
    
    rng = np.random.default_rng(seed=seed)
    n_hops = int(duration * hop_rate)
    hop_sequence = rng.integers(0, n_channels, size=n_hops)
    
    iq = np.zeros(total_samples, dtype=np.complex64)
    for i, hop_idx in enumerate(hop_sequence):
        start = i * samples_per_hop
        end = min(start + samples_per_hop, total_samples)
        n = end - start
        t = np.arange(n) / fs
        iq[start:end] = np.exp(1j * 2 * np.pi * hop_freqs[hop_idx] * t).astype(np.complex64)
    
    return iq, hop_sequence, hop_freqs

# Generate the signal
iq, hop_seq_true, hop_freqs = simulate_fhss()
print(f"Generated {len(iq)/4e6*1e3:.1f} ms of FHSS signal")
print(f"First 10 hops (channel indices): {hop_seq_true[:10]}")
print(f"Channel frequencies: {hop_freqs/1e3} kHz")

A.2 Visualize with Spectrogram

def plot_fhss_spectrogram(iq, fs, hop_rate, n_channels=16, title="FHSS Spectrogram"):
    """Generate and save FHSS spectrogram."""
    f, t, Sxx = spectrogram(iq, fs=fs, nperseg=256, noverlap=200, return_onesided=False)
    f_shifted = np.fft.fftshift(f) / 1e3  # kHz
    Sxx_shifted = np.fft.fftshift(Sxx, axes=0)
    
    plt.figure(figsize=(14, 5))
    plt.pcolormesh(t * 1e3, f_shifted,
                   10 * np.log10(Sxx_shifted + 1e-10),
                   shading='gouraud', cmap='hot', vmin=-40, vmax=0)
    plt.xlabel("Time (ms)")
    plt.ylabel("Frequency (kHz from center)")
    plt.colorbar(label="Power (dB)")
    plt.title(title)
    plt.tight_layout()
    plt.savefig('lab8/plots/fhss_spectrogram.png', dpi=150)
    print("Saved: lab8/plots/fhss_spectrogram.png")
    return f_shifted, t, Sxx_shifted

f_kHz, t_s, Sxx = plot_fhss_spectrogram(iq, fs=4e6, hop_rate=1600)

A.3 Recover the Hop Sequence

def recover_hop_sequence(Sxx, f_kHz, t_s, hop_freqs_hz, fs):
    """
    Recover the hop sequence from a spectrogram by finding the peak frequency
    at each time slice.
    
    Returns:
        recovered_hops: array of recovered channel indices (one per time slice)
        hop_times: center time of each spectrogram time slice
    """
    hop_freqs_kHz = hop_freqs_hz / 1e3
    
    recovered_hops = []
    for t_idx in range(Sxx.shape[1]):
        col = Sxx[:, t_idx]
        peak_f = f_kHz[np.argmax(col)]
        
        # Find nearest hop channel
        dists = np.abs(hop_freqs_kHz - peak_f)
        recovered_hops.append(np.argmin(dists))
    
    return np.array(recovered_hops), t_s

# Note: the spectrogram has more time slices than hops (256-sample window at 4MSPS)
# Each hop is 2500 samples long; each spectrogram slice is 256-200=56 samples
recovered, times = recover_hop_sequence(Sxx, f_kHz, t_s, hop_freqs, fs=4e6)
print(f"Recovered {len(recovered)} spectrogram slices from {len(hop_seq_true)} true hops")

# Downsample recovered sequence to one-per-hop for comparison
slices_per_hop = len(recovered) // len(hop_seq_true)
recovered_downsampled = recovered[slices_per_hop//2::slices_per_hop][:len(hop_seq_true)]

accuracy = np.mean(recovered_downsampled == hop_seq_true)
print(f"Hop sequence recovery accuracy: {accuracy:.1%}")
print(f"True hops    (first 10): {hop_seq_true[:10]}")
print(f"Recovered    (first 10): {recovered_downsampled[:10]}")

Part A deliverable: The spectrogram plot showing the FHSS pattern, and your measured hop sequence recovery accuracy. Include the printed output from A.3.

Part B: Chirp Spread Spectrum and Processing Gain Verification (7 points)

B.1 Generate a Chirp Signal

#!/usr/bin/env python3
"""Lab 8 Part B: Chirp spread spectrum and matched filter."""
import numpy as np
import matplotlib.pyplot as plt

def generate_chirp(f_start, f_end, T, fs):
    """Generate a linear frequency-modulated (LFM) chirp."""
    t = np.arange(int(T * fs)) / fs
    k = (f_end - f_start) / T
    phase = 2 * np.pi * (f_start * t + 0.5 * k * t**2)
    return np.exp(1j * phase).astype(np.complex64)

# Parameters matching LoRa SF7, BW=125kHz
BW = 125e3    # chirp bandwidth (Hz)
T = 2.048e-3  # chirp duration (seconds) -- LoRa SF7 at 125 kHz
fs = 1e6      # sample rate

chirp = generate_chirp(f_start=-BW/2, f_end=BW/2, T=T, fs=fs)
print(f"Chirp duration: {T*1e3:.3f} ms")
print(f"Chirp bandwidth: {BW/1e3:.0f} kHz")
print(f"Time-bandwidth product TW: {T * BW:.1f}")
print(f"Expected processing gain: {10 * np.log10(T * BW):.1f} dB")

B.2 Matched Filter Reception

def matched_filter_receive(received, chirp_template):
    """
    Apply matched filter by correlating received signal with conjugated template.
    Returns the correlation magnitude; the peak marks the chirp's start.
    """
    template_conj = np.conj(chirp_template[::-1])
    correlation = np.convolve(received, template_conj, mode='full')
    return np.abs(correlation)

def measure_processing_gain(chirp, snr_input_db, n_trials=100, fs=1e6):
    """
    Measure the SNR improvement from matched filter reception.
    
    Add AWGN at snr_input_db, apply matched filter, measure output SNR.
    Average over n_trials for reliability.
    """
    # Signal power
    sig_power = np.mean(np.abs(chirp)**2)
    
    # Noise power that gives the target input SNR
    snr_linear = 10 ** (snr_input_db / 10)
    noise_sigma = np.sqrt(sig_power / snr_linear / 2)
    
    snr_outputs = []
    for _ in range(n_trials):
        # Add AWGN to the chirp
        noise = noise_sigma * (np.random.randn(len(chirp)) + 1j * np.random.randn(len(chirp)))
        received = chirp + noise
        
        # Apply matched filter
        mf_out = matched_filter_receive(received, chirp)
        
        # Measure peak SNR
        peak_val = np.max(mf_out)
        # Background = noise power in matched filter output (away from peak)
        n_background = int(0.1 * len(mf_out))
        background = mf_out[:n_background]
        noise_rms = np.std(background)
        
        snr_out = 20 * np.log10(peak_val / (noise_rms + 1e-12))
        snr_outputs.append(snr_out)
    
    return np.mean(snr_outputs)

# Measure processing gain at -10 dB input SNR
snr_in_db = -10
snr_out_db = measure_processing_gain(chirp, snr_in_db, n_trials=50, fs=fs)
pg_measured = snr_out_db - snr_in_db
pg_theoretical = 10 * np.log10(T * BW)

print(f"\nProcessing gain verification:")
print(f"  Input SNR: {snr_in_db} dB")
print(f"  Output SNR (measured): {snr_out_db:.1f} dB")
print(f"  Processing gain (measured): {pg_measured:.1f} dB")
print(f"  Processing gain (formula 10·log10(TW)): {pg_theoretical:.1f} dB")
print(f"  Discrepancy: {abs(pg_measured - pg_theoretical):.1f} dB")

B.3 Spectrogram of Chirp Signal

from scipy.signal import spectrogram as sp_spectrogram

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Time-domain amplitude
ax = axes[0]
t_ax = np.arange(len(chirp)) / fs * 1e3
ax.plot(t_ax, chirp.real[:int(0.2 * len(chirp))])
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Amplitude (real part)")
ax.set_title("Chirp Signal (first 20%)")
ax.grid(True)

# Spectrogram
ax = axes[1]
f_sg, t_sg, Sxx = sp_spectrogram(chirp, fs=fs, nperseg=64, noverlap=60, return_onesided=False)
f_sg_shift = np.fft.fftshift(f_sg) / 1e3
Sxx_shift = np.fft.fftshift(Sxx, axes=0)
ax.pcolormesh(t_sg * 1e3, f_sg_shift, 10 * np.log10(Sxx_shift + 1e-10),
              shading='gouraud', cmap='hot')
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Frequency (kHz)")
ax.set_title("Chirp Spectrogram (linear sweep visible)")
ax.grid(True)

plt.tight_layout()
plt.savefig('lab8/plots/chirp_analysis.png', dpi=150)
print("Saved: lab8/plots/chirp_analysis.png")

Part B deliverable: The chirp spectrogram plot, and the printed processing gain verification output showing measured vs. theoretical gain. If the discrepancy exceeds 3 dB, describe why.

Part C: Anti-Jamming Margin Budget (4 points)

A communications system uses FHSS with the following parameters:

Calculate the following:

def jamming_margin_budget(N_h, P_jammer_dBm, P_signal_dBm, EbN0_required_dB):
    """
    Calculate FHSS anti-jamming margin.
    
    J/S = P_jammer - P_signal  (jammer-to-signal ratio in dB)
    PG  = 10 * log10(N_h)      (processing gain)
    Margin = PG - EbN0_required - J/S
    
    Positive margin: system survives the jammer.
    Negative margin: jammer succeeds.
    """
    J_over_S_dB = P_jammer_dBm - P_signal_dBm
    PG_dB = 10 * np.log10(N_h)
    margin_dB = PG_dB - EbN0_required_dB - J_over_S_dB
    
    print(f"Jamming Margin Budget:")
    print(f"  J/S ratio: {J_over_S_dB:.1f} dB")
    print(f"  Processing gain (N_h={N_h} channels): {PG_dB:.1f} dB")
    print(f"  Required Eb/N0: {EbN0_required_dB:.1f} dB")
    print(f"  Anti-jamming margin: {margin_dB:.1f} dB")
    print(f"  Result: {'System SURVIVES jammer' if margin_dB > 0 else 'Jammer SUCCEEDS'}")
    return margin_dB

# Compute for the given scenario
margin = jamming_margin_budget(
    N_h=32,
    P_jammer_dBm=10,
    P_signal_dBm=-20,
    EbN0_required_dB=8
)

# Now compute the N_h required for a positive margin in this scenario
print("\n--- Minimum channels for positive margin ---")
for n in [32, 64, 128, 256, 512]:
    m = jamming_margin_budget(n, P_jammer_dBm=10, P_signal_dBm=-20, EbN0_required_dB=8)
    if m > 0:
        print(f">>> N_h = {n} channels is sufficient <<<")
        break

Part C deliverable: The printed budget output for N_h = 32, and the minimum N_h required for a positive margin. Answer: at what power ratio (J/S) does a 32-channel FHSS system break even?

Part D: Reflection (2 points)

Write 150-200 words in your lab report answering:

"Your SIGINT pipeline from Lab 7 would classify an FHSS signal as 'frequency hopping' at Stage 3 (via the spectrogram) and a chirp signal as 'unknown wideband noise' at Stage 2. Explain which stage of the pipeline fails first for each waveform, and what additional measurement beyond the current Stage 1-4 toolset would recover the modulation type for each."

Your answer should reference specific stages (1-4) and specific observations from this lab.

Lab Report

Create lab-8-report.md with:

1. plots/fhss_spectrogram.png -- FHSS spectrogram with hop pattern visible
2. plots/chirp_analysis.png -- Chirp time-domain and spectrogram
3. Hop sequence recovery accuracy and a brief explanation of why it is not 100%
4. Processing gain: measured value, theoretical value, and discrepancy explanation (if any)
5. Jamming margin budget: full printed output for N_h = 32 and minimum N_h result
6. Reflection (Part D)

Grading