From 42452f32bd2b97cdc0dccadc8ecacdd97d8c4c67 Mon Sep 17 00:00:00 2001 From: vishwaksen-1 Date: Thu, 5 Feb 2026 21:46:11 +0530 Subject: [PATCH 1/2] "Add real-time packet detection section to detection chapter" --- _images/detection_realtime.svg | 95912 ++++++++++++++++ content/detection.rst | 560 + .../detection_realtime.py | 154 + 3 files changed, 96626 insertions(+) create mode 100644 _images/detection_realtime.svg create mode 100644 figure-generating-scripts/detection_realtime.py diff --git a/_images/detection_realtime.svg b/_images/detection_realtime.svg new file mode 100644 index 00000000..cd1c80d9 --- /dev/null +++ b/_images/detection_realtime.svg @@ -0,0 +1,95912 @@ + + + + + + + + 2026-02-05T20:17:37.262558 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/content/detection.rst b/content/detection.rst index b9cf5eff..c8b6e591 100644 --- a/content/detection.rst +++ b/content/detection.rst @@ -410,3 +410,563 @@ In a DSSS system, the receiver's ability to recover data is entirely dependent o :alt: DSSS The peak occurs at zero offset as expected, and it drops linearly, i.e. it gets to half the peak value at a half-chip offset. After more than one chip offset the correlation might seem like it's going back up, but the actual peak will be low because it's not aligned to the sequence anymore. + +**************************************************** +Real-Time Packet Detection in Continuous IQ Streams +**************************************************** + +So far we have explored the theoretical foundations of signal detection, from correlators to CFAR detectors to spread spectrum systems. Now we bring it all together to solve a common practical problem: **detecting intermittent packets in a continuous stream of IQ samples from an SDR**. + +The Challenge +############# + +Consider this scenario: You have a modem or IoT device that transmits a data packet once per second (or at irregular intervals). Your SDR is continuously receiving samples at, say, 1 MHz. The packets arrive at unpredictable times, buried in noise and interference. You need to: + +1. Detect when a packet arrives +2. Determine the exact sample index where it starts +3. Extract the packet for further processing (demodulation, decoding, etc.) +4. Do this in real-time without missing packets + +This is fundamentally different from processing a pre-recorded IQ file where you can analyze the entire signal at once. Here, samples arrive continuously, and you must make decisions in real-time with limited computational resources. + +Real-World Applications +####################### + +This pattern appears everywhere in wireless systems: + +- **IoT Sensor Networks**: LoRa, Sigfox, and NB-IoT devices that wake up periodically to transmit sensor data +- **Amateur Radio Packet Networks**: APRS beacons transmitting GPS position reports +- **Satellite Communications**: LEO satellites passing overhead transmitting telemetry bursts +- **Industrial Wireless**: Factory sensors reporting only when measurements exceed thresholds +- **Smart Home Devices**: Zigbee, Z-Wave devices transmitting when button pressed or state changes + +In all these cases, the receiver must remain vigilant, detecting packets when they arrive without knowing precisely when that will be. + +The Building Blocks +################### + +We will combine several techniques covered in this chapter: + +1. **Cross-Correlation**: To find the known preamble pattern +2. **CFAR Detection**: To adaptively set thresholds despite varying noise +3. **Buffer Management**: To handle continuous streaming data +4. **Peak Detection**: To extract precise packet timing + +The key insight is that we don't need to process every single sample individually. Instead, we: +- Accumulate samples in **buffers** (chunks of, say, 100,000 samples) +- Run our detector on each buffer +- Extract detected packets +- Maintain state across buffer boundaries to avoid missing packets that span buffers + +Implementation Strategy +####################### + +Our detector will follow this workflow: + +.. code-block:: text + + ┌─────────────────────────────────────────────────────────────┐ + │ Continuous IQ Stream from SDR (e.g., 1 Msps) │ + └────────────────────┬────────────────────────────────────────┘ + │ + ▼ + ┌─────────────────────────────────────────────────────────────┐ + │ Buffer Accumulation (e.g., 100k samples = 0.1 sec) │ + └────────────────────┬────────────────────────────────────────┘ + │ + ▼ + ┌─────────────────────────────────────────────────────────────┐ + │ Cross-Correlation with Known Preamble │ + │ → Produces correlation vs. sample index │ + └────────────────────┬────────────────────────────────────────┘ + │ + ▼ + ┌─────────────────────────────────────────────────────────────┐ + │ CFAR Threshold Computation │ + │ → Adaptive threshold that tracks noise floor │ + └────────────────────┬────────────────────────────────────────┘ + │ + ▼ + ┌─────────────────────────────────────────────────────────────┐ + │ Peak Detection (correlation > threshold) │ + │ → List of candidate packet start indices │ + └────────────────────┬────────────────────────────────────────┘ + │ + ▼ + ┌─────────────────────────────────────────────────────────────┐ + │ Packet Extraction & Validation │ + │ → Extract samples, pass to demodulator │ + └─────────────────────────────────────────────────────────────┘ + +**Buffer Overlap Strategy** + +To avoid missing packets that straddle buffer boundaries, we use an **overlap-save** approach: + +- Each buffer includes the last ``N_preamble`` samples from the previous buffer +- This ensures any packet starting near the end of buffer ``i`` will be fully contained in buffer ``i+1`` +- Trade-off: Small computational overhead vs. guaranteed detection + +Python Implementation +##################### + +Let's build a complete packet detector step by step. We'll use a Zadoff-Chu preamble (as introduced earlier) and implement adaptive CFAR detection. + +Step 1: Define the Preamble and Parameters +******************************************* + +.. code-block:: python + + import numpy as np + import matplotlib.pyplot as plt + from scipy.signal import correlate + + # Preamble: Zadoff-Chu sequence (excellent correlation properties) + N_zc = 63 # ZC sequence length (typically prime or power of 2 - 1) + u = 5 # ZC root + t = np.arange(N_zc) + preamble = np.exp(-1j * np.pi * u * t * (t + 1) / N_zc) + + # System parameters + sample_rate = 1e6 + buffer_size = 100000 + overlap_size = len(preamble) # Overlap to catch boundary packets + + # CFAR parameters + cfar_guard = 10 + cfar_train = 50 + pfa_target = 1e-6 + + # Packet parameters (for simulation) + packet_length = 500 # Total packet length in samples (preamble + data) + snr_db = -5 + +Step 2: CFAR Detector Function +******************************* + +We'll use the Cell-Averaging CFAR (CA-CFAR) from earlier, slightly optimized: + +.. code-block:: python + + def ca_cfar_1d(signal, num_train, num_guard, pfa): + """ + 1D Cell-Averaging CFAR detector. + + Args: + signal: Input signal (typically correlation magnitude) + num_train: Number of training cells (on each side) + num_guard: Number of guard cells (on each side) + pfa: Target probability of false alarm + + Returns: + threshold: Adaptive threshold array + """ + n = len(signal) + threshold = np.zeros(n) + alpha = num_train * (pfa**(-1/num_train) - 1) + + for i in range(n): + # Define training window indices + train_start_left = max(0, i - num_guard - num_train) + train_end_left = max(0, i - num_guard) + train_start_right = min(n, i + num_guard + 1) + train_end_right = min(n, i + num_guard + num_train + 1) + + # Collect training cells (avoid guard cells and CUT) + train_cells = np.concatenate([ + signal[train_start_left:train_end_left], + signal[train_start_right:train_end_right] + ]) + + if len(train_cells) > 0: + noise_est = np.mean(train_cells) + threshold[i] = alpha * noise_est + + return threshold + +Step 3: Packet Detection Function +********************************** + +.. code-block:: python + + def detect_packets(buffer, preamble, cfar_guard, cfar_train, pfa, + min_spacing=None): + """ + Detect packets in a buffer of IQ samples. + + Args: + buffer: Complex IQ samples + preamble: Known preamble sequence + cfar_guard: CFAR guard cells + cfar_train: CFAR training cells + pfa: Target false alarm probability + min_spacing: Minimum samples between detections (prevents duplicates) + + Returns: + detections: List of sample indices where packets start + """ + # Correlate buffer with preamble + corr = correlate(buffer, preamble, mode='same') + corr_power = np.abs(corr)**2 + + # Compute adaptive threshold + threshold = ca_cfar_1d(corr_power, cfar_train, cfar_guard, pfa) + + # Find peaks above threshold + detections_raw = np.where(corr_power > threshold)[0] + + # Compensate for correlation offset (peak occurs at len(preamble)//2 after true start) + half_preamble = len(preamble) // 2 + detections_raw = detections_raw - half_preamble + + # Remove edge detections (unreliable) + half_preamble = len(preamble) // 2 + detections_raw = detections_raw[ + (detections_raw > half_preamble) & + (detections_raw < len(buffer) - half_preamble) + ] + + # Remove duplicate detections (peaks close together) + if min_spacing is None: + min_spacing = len(preamble) + + detections = [] + if len(detections_raw) > 0: + detections.append(detections_raw[0]) + for det in detections_raw[1:]: + if det - detections[-1] > min_spacing: + detections.append(det) + + return detections, corr_power, threshold + +Step 4: Simulation - Generate Test Signal +****************************************** + +.. code-block:: python + + def generate_packet_stream(preamble, packet_length, num_packets, + sample_rate, snr_db): + """ + Generate a simulated IQ stream with intermittent packets. + + Returns: + signal: Complex IQ samples + true_starts: Ground truth packet start indices + """ + # Calculate noise power from SNR + signal_power = 1.0 # Normalized preamble power + noise_power = signal_power / (10**(snr_db/10)) + noise_std = np.sqrt(noise_power / 2) # Complex noise + + # Generate QPSK data (random payload after preamble) + qpsk_map = np.array([1+1j, -1+1j, -1-1j, 1-1j]) / np.sqrt(2) + + # Time between packets (1 second +/- 20% jitter) + packets_per_sec = 1 + avg_gap = int(sample_rate / packets_per_sec) + + signal = [] + true_starts = [] + + for i in range(num_packets): + # Add gap (noise only) + if i == 0: + gap_length = np.random.randint(avg_gap//2, avg_gap) + else: + gap_length = np.random.randint(int(avg_gap*0.8), int(avg_gap*1.2)) + + noise = noise_std * (np.random.randn(gap_length) + + 1j*np.random.randn(gap_length)) + signal.extend(noise) + + # Record true packet start + true_starts.append(len(signal)) + + # Add packet (preamble + data) + data_length = packet_length - len(preamble) + data = qpsk_map[np.random.randint(0, 4, data_length)] + packet = np.concatenate([preamble, data]) + + # Add noise to packet + packet_noisy = packet + noise_std * (np.random.randn(len(packet)) + + 1j*np.random.randn(len(packet))) + signal.extend(packet_noisy) + + # Add final gap + gap_length = np.random.randint(avg_gap//2, avg_gap) + noise = noise_std * (np.random.randn(gap_length) + + 1j*np.random.randn(gap_length)) + signal.extend(noise) + + return np.array(signal), true_starts + + # Generate 5 seconds of signal with ~5 packets + signal, true_starts = generate_packet_stream( + preamble, packet_length, num_packets=5, + sample_rate=sample_rate, snr_db=snr_db + ) + + print(f"Generated {len(signal)} samples ({len(signal)/sample_rate:.1f} sec)") + print(f"True packet starts: {true_starts}") + +Step 5: Run Detection in Streaming Mode +**************************************** + +Now we process the signal in chunks, simulating real-time streaming: + +.. code-block:: python + + def process_stream(signal, preamble, buffer_size, overlap_size, + cfar_guard, cfar_train, pfa): + """ + Process continuous IQ stream in buffers (simulates real-time). + + Returns: + all_detections: List of detected packet starts (global indices) + """ + all_detections = [] + n_samples = len(signal) + current_pos = 0 + + while current_pos < n_samples: + # Define buffer with overlap + buffer_start = max(0, current_pos - overlap_size) + buffer_end = min(n_samples, current_pos + buffer_size) + buffer = signal[buffer_start:buffer_end] + + # Detect packets in this buffer + detections, corr_power, threshold = detect_packets( + buffer, preamble, cfar_guard, cfar_train, pfa + ) + + # Convert buffer-relative indices to global indices + for det in detections: + global_idx = buffer_start + det + + # Avoid duplicate detections from overlap region + if len(all_detections) == 0 or \ + global_idx - all_detections[-1] > len(preamble): + all_detections.append(global_idx) + + current_pos += buffer_size + + return all_detections + + + detected_starts = process_stream( + signal, preamble, buffer_size, overlap_size, + cfar_guard, cfar_train, pfa_target + ) + + print(f"\nDetection Results:") + print(f"True packets: {len(true_starts)}") + print(f"Detected packets: {len(detected_starts)}") + print(f"Detected starts: {detected_starts}") + +Step 6: Evaluate Performance +***************************** + +.. code-block:: python + + # Calculate detection statistics + tolerance = len(preamble) + + matched_detections = [] + false_alarms = [] + + for det in detected_starts: + # Check if detection matches any true packet + matched = False + for true_start in true_starts: + if abs(det - true_start) <= tolerance: + matched_detections.append(det) + matched = True + break + if not matched: + false_alarms.append(det) + + missed_packets = len(true_starts) - len(matched_detections) + + print(f"\nPerformance Metrics:") + print(f" Correct detections: {len(matched_detections)}/{len(true_starts)}") + print(f" Missed packets: {missed_packets}") + print(f" False alarms: {len(false_alarms)}") + + # Calculate timing errors + timing_errors = [] + for det in matched_detections: + errors = [abs(det - ts) for ts in true_starts] + timing_errors.append(min(errors)) + + if len(timing_errors) > 0: + print(f" Timing error (avg): {np.mean(timing_errors):.1f} samples") + print(f" Timing error (max): {np.max(timing_errors):.1f} samples") + +Step 7: Visualize Results +************************** + +.. code-block:: python + + # Process one buffer for detailed visualization + buffer_start = max(0, true_starts[0] - 5000) + buffer_end = min(len(signal), true_starts[0] + 20000) + viz_buffer = signal[buffer_start:buffer_end] + + detections_viz, corr_viz, thresh_viz = detect_packets( + viz_buffer, preamble, cfar_guard, cfar_train, pfa_target + ) + + # Convert to global indices for plotting + detections_viz_global = [d + buffer_start for d in detections_viz] + + # Create visualization + fig, axes = plt.subplots(3, 1, figsize=(14, 10)) + time_axis = (np.arange(len(viz_buffer)) + buffer_start) / sample_rate * 1000 # ms + + # Subplot 1: Received signal power + axes[0].plot(time_axis, np.abs(viz_buffer)**2, 'gray', alpha=0.6, linewidth=0.5) + axes[0].set_ylabel('Power') + axes[0].set_title('Received IQ Signal Power') + axes[0].grid(True, alpha=0.3) + + # Mark true packet locations + for ts in true_starts: + if buffer_start <= ts <= buffer_end: + t_ms = ts / sample_rate * 1000 + axes[0].axvline(t_ms, color='green', linestyle='--', alpha=0.7, + label='True Packet' if ts == true_starts[0] else '') + axes[0].legend() + + # Subplot 2: Correlation output + axes[1].plot(time_axis, corr_viz, 'blue', linewidth=1, label='Correlation') + axes[1].plot(time_axis, thresh_viz, 'red', linestyle='--', linewidth=1.5, + label='CFAR Threshold') + axes[1].set_ylabel('Correlation Power') + axes[1].set_title('Preamble Correlation with Adaptive CFAR Threshold') + axes[1].grid(True, alpha=0.3) + axes[1].legend() + + # Subplot 3: Detections + detection_mask = np.zeros(len(viz_buffer)) + for det in detections_viz: + detection_mask[det] = corr_viz[det] + + axes[2].plot(time_axis, corr_viz, 'blue', alpha=0.4, linewidth=0.8) + axes[2].scatter(time_axis[detection_mask > 0], detection_mask[detection_mask > 0], + color='lime', edgecolors='black', s=100, zorder=5, + label='Detected Packets') + axes[2].set_xlabel('Time (ms)') + axes[2].set_ylabel('Correlation Power') + axes[2].set_title('Detected Packet Locations') + axes[2].grid(True, alpha=0.3) + axes[2].legend() + + plt.tight_layout() + plt.savefig('../_images/detection_realtime.svg', bbox_inches='tight') + plt.show() + +The visualization should show: + +1. **Top plot**: Raw signal power with true packet locations marked +2. **Middle plot**: Correlation output with adaptive CFAR threshold tracking the noise floor +3. **Bottom plot**: Detected packets highlighted as peaks above threshold + +.. image:: ../_images/detection_realtime.svg + :align: center + :target: ../_images/detection_realtime.svg + :alt: Real-time packet detection results + +Practical Considerations and Tuning +#################################### + +Buffer Size Trade-offs +*********************** + +**Larger buffers (e.g., 1M samples):** + +- ✅ Better CFAR noise estimation (more training cells) +- ✅ Lower computational overhead (fewer processing calls) +- ❌ Higher latency (must wait for buffer to fill) +- ❌ More memory required + +**Smaller buffers (e.g., 10k samples):** + +- ✅ Lower latency (faster response) +- ✅ Less memory usage +- ❌ CFAR performance degrades (fewer training cells) +- ❌ Higher CPU usage (more frequent processing) + +**Recommendation**: Start with buffer size = 10× to 100× your preamble length. For a 63-sample preamble at 1 Msps, try 10k-100k samples. + +CFAR Parameter Tuning +********************** + +The three CFAR parameters control detector behavior: + +**num_guard** (guard cells): + +- Purpose: Prevents signal leakage into noise estimate +- Too small: Signal leaks into training region → raised threshold → missed detections +- Too large: Fewer training cells → poor noise estimate +- **Rule of thumb**: Set to ~0.5 to 1.0× preamble length + +**num_train** (training cells): + +- Purpose: Estimates local noise floor +- Too small: Noisy threshold → false alarms or missed detections +- Too large: Threshold doesn't adapt quickly enough to noise changes +- **Rule of thumb**: Set to ~3 to 5× preamble length + +**pfa** (probability of false alarm): + +- Purpose: Controls detection sensitivity +- Too high (e.g., 1e-2): Many false alarms +- Too low (e.g., 1e-10): Misses weak packets +- **Rule of thumb**: Start with 1e-5 for per-lag PFA, then adjust based on system-level false alarm rate + +Remember the relationship between per-lag and system-level false alarm rates from earlier in the chapter! + +Handling Real SDR Hardware +*************************** + +When connecting to actual SDR hardware (RTL-SDR, PlutoSDR, USRP, etc.), you'll need to: + +1. **Replace simulation with SDR API**: + + .. code-block:: python + + # Example with RTL-SDR + from rtlsdr import RtlSdr + + sdr = RtlSdr() + sdr.sample_rate = 1e6 + sdr.center_freq = 915e6 # ISM band + sdr.gain = 'auto' + + # Read buffer + samples = sdr.read_samples(buffer_size) + +2. **Handle DC offset and IQ imbalance**: Real hardware often has DC spikes and phase imbalance. Pre-process with: + + .. code-block:: python + + # Remove DC offset + samples = samples - np.mean(samples) + + # Optional: High-pass filter to remove residual DC + from scipy.signal import butter, lfilter + b, a = butter(3, 0.01, btype='high') + samples = lfilter(b, a, samples) + +3. **Frequency offset correction**: If your modem's carrier frequency is slightly offset from your SDR's LO, the preamble correlation will be reduced. Consider adding coarse frequency offset estimation or using frequency-resilient correlation (covered earlier in this chapter). + +4. **Sample rate mismatch**: Ensure your SDR sample rate exactly matches the expected symbol rate of your packets. Even small mismatches accumulate over long preambles. + +Extending to Unknown Preambles +******************************* + +What if you don't know the exact preamble? Some options: + +- **Energy detection**: Use the simple ``np.abs()`` and ``np.diff()`` approach suggested earlier, but with CFAR thresholding +- **Blind synchronization**: Exploit statistical properties (cyclostationary features) covered in other chapters +- **Multi-hypothesis detection**: Run multiple correlators in parallel, one for each candidate preamble +- **Machine learning**: Train a neural network to recognize packet boundaries (advanced topic) \ No newline at end of file diff --git a/figure-generating-scripts/detection_realtime.py b/figure-generating-scripts/detection_realtime.py new file mode 100644 index 00000000..2b740301 --- /dev/null +++ b/figure-generating-scripts/detection_realtime.py @@ -0,0 +1,154 @@ +""" +Generate visualization for real-time packet detection section. +""" +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import correlate + +def ca_cfar_1d(signal, num_train, num_guard, pfa): + """Cell-Averaging CFAR detector.""" + n = len(signal) + threshold = np.zeros(n) + alpha = num_train * (pfa**(-1/num_train) - 1) + + for i in range(n): + train_start_left = max(0, i - num_guard - num_train) + train_end_left = max(0, i - num_guard) + train_start_right = min(n, i + num_guard + 1) + train_end_right = min(n, i + num_guard + num_train + 1) + + train_cells = np.concatenate([ + signal[train_start_left:train_end_left], + signal[train_start_right:train_end_right] + ]) + + if len(train_cells) > 0: + threshold[i] = alpha * np.mean(train_cells) + + return threshold + +def detect_packets(buffer, preamble, cfar_guard, cfar_train, pfa): + """Detect packets in IQ buffer.""" + corr = correlate(buffer, preamble, mode='same') + corr_power = np.abs(corr)**2 + threshold = ca_cfar_1d(corr_power, cfar_train, cfar_guard, pfa) + + detections_raw = np.where(corr_power > threshold)[0] + half_preamble = len(preamble) // 2 + detections_raw = detections_raw - half_preamble + detections_raw = detections_raw[ + (detections_raw > half_preamble) & + (detections_raw < len(buffer) - half_preamble) + ] + + detections = [] + if len(detections_raw) > 0: + detections.append(detections_raw[0]) + for det in detections_raw[1:]: + if det - detections[-1] > len(preamble): + detections.append(det) + + return detections, corr_power, threshold + +def generate_packet_stream(preamble, packet_length, num_packets, sample_rate, snr_db): + """Generate simulated IQ stream with packets.""" + signal_power = 1.0 + noise_power = signal_power / (10**(snr_db/10)) + noise_std = np.sqrt(noise_power / 2) + + qpsk_map = np.array([1+1j, -1+1j, -1-1j, 1-1j]) / np.sqrt(2) + avg_gap = int(sample_rate / 1) # 1 packet per second + + signal = [] + true_starts = [] + + for i in range(num_packets): + if i == 0: + gap_length = np.random.randint(avg_gap//2, avg_gap) + else: + gap_length = np.random.randint(int(avg_gap*0.8), int(avg_gap*1.2)) + + noise = noise_std * (np.random.randn(gap_length) + 1j*np.random.randn(gap_length)) + signal.extend(noise) + + true_starts.append(len(signal)) + + data_length = packet_length - len(preamble) + data = qpsk_map[np.random.randint(0, 4, data_length)] + packet = np.concatenate([preamble, data]) + packet_noisy = packet + noise_std * (np.random.randn(len(packet)) + + 1j*np.random.randn(len(packet))) + signal.extend(packet_noisy) + + gap_length = np.random.randint(avg_gap//2, avg_gap) + noise = noise_std * (np.random.randn(gap_length) + 1j*np.random.randn(gap_length)) + signal.extend(noise) + + return np.array(signal), true_starts + +# Generate test signal +np.random.seed(42) +N_zc, u = 63, 5 +t = np.arange(N_zc) +preamble = np.exp(-1j * np.pi * u * t * (t + 1) / N_zc) + +sample_rate = 1e6 +packet_length = 500 +snr_db = -5 + +signal, true_starts = generate_packet_stream( + preamble, packet_length, num_packets=5, + sample_rate=sample_rate, snr_db=snr_db +) + +# Process one buffer for visualization +buffer_start = max(0, true_starts[0] - 5000) +buffer_end = min(len(signal), true_starts[2] + 10000) +viz_buffer = signal[buffer_start:buffer_end] + +detections_viz, corr_viz, thresh_viz = detect_packets( + viz_buffer, preamble, cfar_guard=10, cfar_train=50, pfa=1e-5 +) + +# Create visualization +fig, axes = plt.subplots(3, 1, figsize=(14, 10)) +time_axis = (np.arange(len(viz_buffer)) + buffer_start) / sample_rate * 1000 + +# Subplot 1: Received signal power +axes[0].plot(time_axis, np.abs(viz_buffer)**2, 'gray', alpha=0.6, linewidth=0.5) +axes[0].set_ylabel('Power', fontsize=11) +axes[0].set_title('Received IQ Signal Power', fontsize=12, fontweight='bold') +axes[0].grid(True, alpha=0.3) +for i, ts in enumerate(true_starts): + if buffer_start <= ts <= buffer_end: + t_ms = ts / sample_rate * 1000 + axes[0].axvline(t_ms, color='green', linestyle='--', alpha=0.7, linewidth=1.5, + label='True Packet' if i == 0 else '') +axes[0].legend(fontsize=10) + +# Subplot 2: Correlation output +axes[1].plot(time_axis, corr_viz, 'blue', linewidth=1, label='Correlation') +axes[1].plot(time_axis, thresh_viz, 'red', linestyle='--', linewidth=1.5, label='CFAR Threshold') +axes[1].set_ylabel('Correlation Power', fontsize=11) +axes[1].set_title('Preamble Correlation with Adaptive CFAR Threshold', fontsize=12, fontweight='bold') +axes[1].grid(True, alpha=0.3) +axes[1].legend(fontsize=10) + +# Subplot 3: Detections +detection_mask = np.zeros(len(viz_buffer)) +for det in detections_viz: + detection_mask[det] = corr_viz[det] + +axes[2].plot(time_axis, corr_viz, 'blue', alpha=0.4, linewidth=0.8) +axes[2].scatter(time_axis[detection_mask > 0], detection_mask[detection_mask > 0], + color='lime', edgecolors='black', s=100, zorder=5, label='Detected Packets') +axes[2].set_xlabel('Time (ms)', fontsize=11) +axes[2].set_ylabel('Correlation Power', fontsize=11) +axes[2].set_title('Detected Packet Locations', fontsize=12, fontweight='bold') +axes[2].grid(True, alpha=0.3) +axes[2].legend(fontsize=10) + +plt.tight_layout() +plt.savefig('../_images/detection_realtime.svg', bbox_inches='tight', dpi=150) +print("Figure saved to ../_images/detection_realtime.svg") +plt.show() \ No newline at end of file From fd67f4b6b7f1f445eacd6f402340bb8d6c545640 Mon Sep 17 00:00:00 2001 From: vishwaksen-1 Date: Thu, 5 Feb 2026 22:48:06 +0530 Subject: [PATCH 2/2] Add full E2E QPSK simulation to Sync chapter --- content/sync.rst | 276 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 267 insertions(+), 9 deletions(-) diff --git a/content/sync.rst b/content/sync.rst index 52c92d54..30803934 100644 --- a/content/sync.rst +++ b/content/sync.rst @@ -484,13 +484,271 @@ You can see it's 11 (length of the sequence) in the center, and -1 or 0 for all Another sequence with great autocorrelation properties is Zadoff-Chu sequences, which are used in LTE. They have the benefit of being in sets; you can have multiple different sequences that all have good autocorrelation properties, but they won't trigger each other (i.e., also good cross-correlation properties, when you cross-correlate different sequences in the set). Thanks to that feature, different cell towers will be assigned different sequences so that a phone can not only find the start of the frame but also know which tower it is receiving from. +*********************************************** +Full End-to-End QPSK Communication Example +*********************************************** + +Now that we've learned about time and frequency synchronization, let's put it all together in a complete end-to-end simulation. This example demonstrates a full QPSK communication system that: + +1. Takes an ASCII text message as input +2. Modulates it using QPSK +3. Applies pulse shaping with a raised-cosine filter +4. Simulates a realistic wireless channel (frequency offset, timing delay, and noise) +5. Performs time synchronization using the Mueller & Muller algorithm +6. Performs frequency synchronization using a Costas Loop +7. Demodulates the QPSK symbols +8. Resolves phase ambiguity +9. Decodes the message and calculates Bit Error Rate (BER) - - - - - - - - - +.. code-block:: python + + import numpy as np + import matplotlib.pyplot as plt + from scipy import signal + + # Step 1: Generate Input Data (ASCII Message) + message = "Hello SDR!" + print(f"Original Message: {message}") + + # Convert ASCII to bits + bits = [] + for char in message: + byte = format(ord(char), '08b') # Convert to 8-bit binary + bits.extend([int(b) for b in byte]) + bits = np.array(bits) + print(f"Total bits: {len(bits)}") + + # Step 2: QPSK Modulation + # Map bits to QPSK symbols: 00->-1-1j, 01->-1+1j, 10->1-1j, 11->1+1j + num_symbols = len(bits) // 2 + symbols = np.zeros(num_symbols, dtype=complex) + + for i in range(num_symbols): + bit_pair = (bits[2*i], bits[2*i+1]) + if bit_pair == (0, 0): + symbols[i] = -1-1j + elif bit_pair == (0, 1): + symbols[i] = -1+1j + elif bit_pair == (1, 0): + symbols[i] = 1-1j + else: # (1, 1) + symbols[i] = 1+1j + + # Normalize to unit power + symbols = symbols / np.sqrt(2) + + # Step 3: Pulse Shaping (Raised Cosine Filter) + sps = 8 # samples per symbol + num_taps = 101 + beta = 0.35 + Ts = sps + t = np.arange(-51, 52) + h = np.sinc(t/Ts) * np.cos(np.pi*beta*t/Ts) / (1 - (2*beta*t/Ts)**2) + + # Upsample symbols + upsampled = np.zeros(len(symbols) * sps, dtype=complex) + upsampled[::sps] = symbols + + # Apply pulse shaping + tx_samples = np.convolve(upsampled, h, mode='same') + + # Step 4: Channel Simulation + # Add frequency offset + fs = 1e6 + fo = 800 # 800 Hz frequency offset + Ts_sample = 1/fs + t_vec = np.arange(len(tx_samples)) * Ts_sample + channel_samples = tx_samples * np.exp(1j*2*np.pi*fo*t_vec) + + # Add fractional delay + delay = 0.4 + N_delay = 21 + n = np.arange(-N_delay//2, N_delay//2) + h_delay = np.sinc(n - delay) + h_delay *= np.hamming(N_delay) + h_delay /= np.sum(h_delay) + channel_samples = np.convolve(channel_samples, h_delay, mode='same') + + # Add AWGN noise + SNR_dB = 15 + signal_power = np.mean(np.abs(channel_samples)**2) + noise_power = signal_power / (10**(SNR_dB/10)) + noise = np.sqrt(noise_power/2) * (np.random.randn(len(channel_samples)) + + 1j*np.random.randn(len(channel_samples))) + rx_samples = channel_samples + noise + + # Step 5: Time Synchronization (Mueller & Muller) + # Interpolate for fractional timing recovery + samples_interpolated = signal.resample_poly(rx_samples, 16, 1) + + mu = 0 + out = np.zeros(len(rx_samples) + 10, dtype=complex) + out_rail = np.zeros(len(rx_samples) + 10, dtype=complex) + i_in = 0 + i_out = 2 + + while i_out < len(rx_samples) and i_in*16+16 < len(samples_interpolated): + out[i_out] = samples_interpolated[i_in*16 + int(mu*16)] + out_rail[i_out] = (np.real(out[i_out])>0)*2-1 + 1j*((np.imag(out[i_out])>0)*2-1) + x = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1]) + y = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1]) + mm_val = np.real(y - x) + mu += sps + 0.01*mm_val # Reduced gain for stability + i_in += int(np.floor(mu)) + mu = mu - np.floor(mu) + i_out += 1 + + out = out[2:i_out] + time_synced = out + + # Step 6: Frequency Synchronization (Costas Loop) + N = len(time_synced) + phase = 0 + freq = 0 + alpha = 0.01 # Reduced proportional gain to prevent oscillations + beta = 0.001 # Reduced integral gain + freq_synced = np.zeros(N, dtype=complex) + + for i in range(N): + freq_synced[i] = time_synced[i] * np.exp(-1j*phase) + error = np.real(freq_synced[i]) * np.imag(freq_synced[i]) + + freq += beta * error + phase += freq + alpha * error + + while phase >= 2*np.pi: + phase -= 2*np.pi + while phase < 0: + phase += 2*np.pi + + # Step 7: QPSK Demodulation + # Skip first 35 symbols (allow sync algorithms to fully converge) + skip = min(35, len(freq_synced) // 2) # At least 35 symbols or half the data + demod_symbols = freq_synced[skip:] + + # Try all 4 phase rotations and pick the one with lowest error + best_ber = 1.0 + best_rotation = 0 + best_bits = None + + for phase_rotation in range(4): + # Rotate by 0, π/2, π, 3π/2 + rotation = np.exp(1j * np.pi / 2 * phase_rotation) + rotated_symbols = demod_symbols * rotation + + # Decode QPSK symbols back to bits + test_bits = [] + for sym in rotated_symbols: + I = np.real(sym) + Q = np.imag(sym) + + if I >= 0 and Q >= 0: + test_bits.extend([1, 1]) + elif I < 0 and Q >= 0: + test_bits.extend([0, 1]) + elif I >= 0 and Q < 0: + test_bits.extend([1, 0]) + else: + test_bits.extend([0, 0]) + + test_bits = np.array(test_bits) + if len(test_bits) > len(bits): + test_bits = test_bits[:len(bits)] + elif len(test_bits) < len(bits): + test_bits = np.pad(test_bits, (0, len(bits) - len(test_bits))) + + # Calculate BER for this rotation + ber_test = np.sum(bits != test_bits) / len(bits) + if ber_test < best_ber: + best_ber = ber_test + best_rotation = phase_rotation + best_bits = test_bits + + received_bits = best_bits + print(f"Best phase rotation: {best_rotation * 90}°") + + # Step 8: Decode ASCII and Calculate BER + received_bits = np.array(received_bits) + + # Trim to match original length + if len(received_bits) > len(bits): + received_bits = received_bits[:len(bits)] + elif len(received_bits) < len(bits): + # Pad with zeros if we have fewer bits + received_bits = np.pad(received_bits, (0, len(bits) - len(received_bits))) + + # Calculate BER + bit_errors = np.sum(bits != received_bits) + ber = bit_errors / len(bits) + + # Decode received message + received_message = "" + for i in range(len(received_bits) // 8): + byte_bits = received_bits[i*8:(i+1)*8] + byte_str = ''.join([str(b) for b in byte_bits]) + received_message += chr(int(byte_str, 2)) + + # Step 9: Display Results + print("\n" + "="*50) + print("E2E QPSK Communication System Results") + print("="*50) + print(f"Transmitted Message: {message}") + print(f"Received Message: {received_message}") + print(f"SNR: {SNR_dB} dB") + print(f"Frequency Offset: {fo} Hz") + print(f"Bit Error Rate: {ber:.6f} ({bit_errors}/{len(bits)} bit errors)") + print("="*50) + + # Plot constellation before and after sync + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) + ax1.plot(np.real(rx_samples[::sps]), np.imag(rx_samples[::sps]), '.') + ax1.set_title('Before Synchronization') + ax1.set_xlabel('I') + ax1.set_ylabel('Q') + ax1.grid(True) + ax1.axis('equal') + + ax2.plot(np.real(freq_synced[skip:]), np.imag(freq_synced[skip:]), '.') + ax2.set_title('After Synchronization') + ax2.set_xlabel('I') + ax2.set_ylabel('Q') + ax2.grid(True) + ax2.axis('equal') + + plt.tight_layout() + plt.savefig('sync_e2e_qpsk_constellation.png', dpi=150, bbox_inches='tight') + plt.show() + +When you run this simulation, you should see: + +- The original transmitted message printed to the console +- The received message (which should closely match the original, depending on SNR and channel conditions) +- The Bit Error Rate (BER), showing the ratio of incorrectly received bits +- Two constellation diagrams side-by-side: + + * **Left plot**: The received signal before synchronization, showing the effects of frequency offset, timing errors, and noise. The constellation points will appear rotated and spread out. + * **Right plot**: The signal after both time and frequency synchronization. The QPSK constellation should show four distinct clusters near the ideal symbol locations (±1±j)/√2. + +**Understanding the Results** + +This example demonstrates several important concepts: + +- **Phase Ambiguity**: QPSK (and most modulations) have inherent phase ambiguity. The Costas Loop can lock to any of four possible phase offsets (0°, 90°, 180°, 270°). The code resolves this by trying all four rotations and selecting the one with the lowest BER. + +- **Convergence Time**: The synchronization algorithms need time to converge. That's why we skip the first 35 symbols before demodulation. In a real system, you would use a preamble or training sequence for synchronization. + +- **SNR Impact**: Try adjusting the :code:`SNR_dB` parameter. Lower SNR values will increase the BER, while higher values should result in near-perfect reception. + +- **Frequency Offset**: The :code:`fo` parameter simulates oscillator drift or Doppler shift. Larger offsets make synchronization more challenging. + +**Exercises** + +Try modifying the code to explore different scenarios: + +1. Change the SNR to see how it affects BER (try values from 5 dB to 25 dB) +2. Increase the frequency offset to 2000 Hz or 5000 Hz and observe the impact +3. Modify the message to be longer and see if the BER changes +4. Replace QPSK with BPSK by changing the symbol mapping (you'll need to adjust the demodulator accordingly) +5. Add a plot showing BER vs SNR by running the simulation multiple times + +This simulation ties together all the concepts we've covered in this chapter: pulse shaping, channel effects, timing recovery, and frequency synchronization. It demonstrates that even with impairments like frequency offset, timing errors, and noise, we can successfully recover digital data using the synchronization techniques we've learned. \ No newline at end of file