Lab 11: ML Modulation Classifier (Optional)

Week: 11c -- ML-Based Modulation Classification
Points: 20 (optional; see scoring note)
Time estimate: 90 min lab + 2 hr independent
Prerequisites: RadioML 2016.10A dataset (download separately, ~54 MB); torch, sklearn
Deliverable: lab-11-report.md + trained model + plots

Scoring Note

Lab 11 is optional. If completed, the 20 points can replace your lowest non-capstone lab score. If not completed, there is no penalty. The dataset download and GPU training requirement make this lab unsuitable for all lab environments; complete it if your workstation supports it.

Setup

# Install required packages (if not already installed)
pip install torch torchvision scikit-learn

# Download RadioML 2016.10A dataset
# Register at https://radioml.com and download RML2016.10a.tar.bz2 (~54 MB)
# Place the extracted file as: radioml/RML2016.10a_dict.pkl

mkdir -p lab11/plots lab11/models

Confirm setup:

import pickle, torch
import numpy as np

with open('radioml/RML2016.10a_dict.pkl', 'rb') as f:
    dataset = pickle.load(f, encoding='latin1')

modulations = sorted(set(m for m, _ in dataset.keys()))
snrs = sorted(set(s for _, s in dataset.keys()))
print(f"Modulations: {modulations}")
print(f"SNR range: {min(snrs)} to {max(snrs)} dB")
n_examples = sum(v.shape[0] for v in dataset.values())
print(f"Total examples: {n_examples:,}")
print(f"Sample shape: {list(dataset.values())[0].shape}")  # (1000, 2, 128)
print(f"PyTorch available: CUDA={torch.cuda.is_available()}, CPU=ok")

Part A: Dataset Loading and Exploration (3 points)

#!/usr/bin/env python3
"""Lab 11 Part A: Dataset loading and class distribution."""
import pickle
import numpy as np
import matplotlib.pyplot as plt

def load_radioml_2016(path='radioml/RML2016.10a_dict.pkl'):
    """
    Load RadioML 2016.10A dataset.
    Returns X (N, 2, 128) and parallel lists of modulations and SNRs.
    """
    with open(path, 'rb') as f:
        data = pickle.load(f, encoding='latin1')
    
    X_list, mods, snrs = [], [], []
    for (mod, snr), frames in data.items():
        for frame in frames:
            X_list.append(frame.astype(np.float32))
            mods.append(mod)
            snrs.append(snr)
    
    return np.array(X_list), mods, snrs

X, mods, snrs = load_radioml_2016()
modulation_list = sorted(set(mods))
snr_values = sorted(set(snrs))
mod2idx = {m: i for i, m in enumerate(modulation_list)}

print(f"Dataset loaded: {X.shape} (N, channels, time)")
print(f"Modulation classes: {modulation_list}")

# Class distribution
counts = {m: mods.count(m) for m in modulation_list}

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

ax = axes[0]
ax.bar(range(len(modulation_list)), [counts[m] for m in modulation_list])
ax.set_xticks(range(len(modulation_list)))
ax.set_xticklabels(modulation_list, rotation=45, ha='right')
ax.set_ylabel("Example count")
ax.set_title("Examples per modulation class")
ax.grid(True, axis='y')

ax = axes[1]
snr_counts = {s: snrs.count(s) for s in snr_values}
ax.bar(snr_values, [snr_counts[s] for s in snr_values])
ax.set_xlabel("SNR (dB)")
ax.set_ylabel("Example count")
ax.set_title("Examples per SNR level")
ax.grid(True, axis='y')

plt.tight_layout()
plt.savefig('lab11/plots/dataset_distribution.png', dpi=150)
print("\nSaved: lab11/plots/dataset_distribution.png")

# Per-SNR, per-class example count
per_snr_class = {}
for mod, snr in zip(mods, snrs):
    per_snr_class[(mod, snr)] = per_snr_class.get((mod, snr), 0) + 1

example_counts = set(per_snr_class.values())
print(f"Examples per (modulation, SNR) pair: {example_counts}")

Part A deliverable: The distribution plot, and a statement of how many examples per (modulation, SNR) pair. Is the dataset balanced? What SNR range do you expect most real ISM-band 433 MHz captures to fall in?

Part B: Model Training (8 points)

#!/usr/bin/env python3
"""Lab 11 Part B: Train ModulationCNN on RadioML 2016.10A."""
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader, random_split
import torch.optim as optim
import matplotlib.pyplot as plt

class ModulationCNN(nn.Module):
    """
    CNN for automatic modulation classification.
    Input: (batch, 2, 128) IQ frame
    Output: (batch, n_classes) log-probabilities
    """
    def __init__(self, n_classes=11, dropout=0.5):
        super().__init__()
        self.conv_block = nn.Sequential(
            nn.Conv1d(2, 64, kernel_size=3, padding=1),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Dropout(dropout),
            nn.Conv1d(64, 64, kernel_size=3, padding=1),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Dropout(dropout),
            nn.Conv1d(64, 128, kernel_size=3, padding=1),
            nn.BatchNorm1d(128),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Dropout(dropout),
        )
        self.classifier = nn.Sequential(
            nn.Flatten(),
            nn.Linear(128 * 16, 256),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(256, n_classes),
        )
    
    def forward(self, x):
        return self.classifier(self.conv_block(x))

class RadioMLDataset(Dataset):
    def __init__(self, X, mods, snrs, mod2idx, snr_min=None):
        mask = np.ones(len(X), dtype=bool)
        if snr_min is not None:
            mask = np.array([s >= snr_min for s in snrs])
        
        self.X = torch.from_numpy(X[mask])
        self.Y = torch.tensor([mod2idx[mods[i]] for i in np.where(mask)[0]])
    
    def __len__(self):
        return len(self.X)
    
    def __getitem__(self, idx):
        return self.X[idx], self.Y[idx]

def run_epoch(model, loader, optimizer, criterion, device, train=True):
    model.train(train)
    total_loss, correct, total = 0, 0, 0
    context = torch.no_grad() if not train else torch.enable_grad()
    with context:
        for X_b, Y_b in loader:
            X_b, Y_b = X_b.to(device), Y_b.to(device)
            if train:
                optimizer.zero_grad()
            logits = model(X_b)
            loss = criterion(logits, Y_b)
            if train:
                loss.backward()
                optimizer.step()
            total_loss += loss.item() * len(Y_b)
            correct += (logits.argmax(1) == Y_b).sum().item()
            total += len(Y_b)
    return total_loss / total, correct / total

# Full training pipeline
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"Training on: {device}")

dataset = RadioMLDataset(X, mods, snrs, mod2idx)
n_train = int(0.8 * len(dataset))
train_ds, val_ds = random_split(dataset, [n_train, len(dataset) - n_train],
                                 generator=torch.Generator().manual_seed(42))

train_loader = DataLoader(train_ds, batch_size=256, shuffle=True, num_workers=2, pin_memory=True)
val_loader = DataLoader(val_ds, batch_size=256, shuffle=False, num_workers=2)

model = ModulationCNN(n_classes=len(modulation_list)).to(device)
optimizer = optim.Adam(model.parameters(), lr=1e-3)
criterion = nn.CrossEntropyLoss()

n_epochs = 10
history = {'train_loss': [], 'val_loss': [], 'train_acc': [], 'val_acc': []}

for epoch in range(n_epochs):
    tr_loss, tr_acc = run_epoch(model, train_loader, optimizer, criterion, device, train=True)
    va_loss, va_acc = run_epoch(model, val_loader, optimizer, criterion, device, train=False)
    history['train_loss'].append(tr_loss)
    history['val_loss'].append(va_loss)
    history['train_acc'].append(tr_acc)
    history['val_acc'].append(va_acc)
    print(f"Epoch {epoch+1:02d}: train={tr_loss:.4f}/{tr_acc:.3f} val={va_loss:.4f}/{va_acc:.3f}")

# Save the model
torch.save(model.state_dict(), 'lab11/models/modulation_cnn.pt')
print("\nModel saved: lab11/models/modulation_cnn.pt")

# Plot learning curves
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(history['train_loss'], label='Train')
axes[0].plot(history['val_loss'], label='Val')
axes[0].set_xlabel("Epoch")
axes[0].set_ylabel("Loss")
axes[0].set_title("Training Loss")
axes[0].legend()
axes[0].grid(True)

axes[1].plot(history['train_acc'], label='Train')
axes[1].plot(history['val_acc'], label='Val')
axes[1].set_xlabel("Epoch")
axes[1].set_ylabel("Accuracy")
axes[1].set_title("Training Accuracy")
axes[1].legend()
axes[1].grid(True)

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

print(f"\nFinal test accuracy at SNR=+10dB (all classes): ", end="")
# Filter to +10 dB
snr10_mask = np.array([s == 10 for s in snrs])
X_10 = torch.from_numpy(X[snr10_mask]).to(device)
Y_10 = torch.tensor([mod2idx[mods[i]] for i in np.where(snr10_mask)[0]]).to(device)
model.eval()
with torch.no_grad():
    logits = model(X_10)
    acc_10 = (logits.argmax(1) == Y_10).float().mean().item()
print(f"{acc_10:.3f}")

Part B deliverable: Learning curves plot, and the final validation accuracy plus the SNR=+10dB accuracy. The SNR=+10dB accuracy should exceed 80% after 10 epochs with this architecture.

Part C: Confusion Matrix Analysis (6 points)

#!/usr/bin/env python3
"""Lab 11 Part C: Confusion matrix at two SNR levels."""
import numpy as np
import torch
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay

# Load trained model
model = ModulationCNN(n_classes=len(modulation_list))
model.load_state_dict(torch.load('lab11/models/modulation_cnn.pt', map_location='cpu'))
model.eval()

def get_predictions_at_snr(target_snr):
    """Collect predictions for all examples at a specific SNR."""
    mask = np.array([s == target_snr for s in snrs])
    X_snr = torch.from_numpy(X[mask])
    Y_true = np.array([mod2idx[mods[i]] for i in np.where(mask)[0]])
    
    with torch.no_grad():
        logits = model(X_snr)
        Y_pred = logits.argmax(1).numpy()
    
    acc = (Y_pred == Y_true).mean()
    return Y_true, Y_pred, acc

# Evaluate at +10 dB and -5 dB
print("Computing confusion matrices...")
Y_true_10, Y_pred_10, acc_10 = get_predictions_at_snr(10)
Y_true_n5, Y_pred_n5, acc_n5 = get_predictions_at_snr(-5)

print(f"Accuracy at SNR = +10 dB: {acc_10:.3f}")
print(f"Accuracy at SNR =  -5 dB: {acc_n5:.3f}")

# Plot confusion matrices
fig, axes = plt.subplots(1, 2, figsize=(18, 7))

for (Y_t, Y_p, acc, snr_val, ax) in [
    (Y_true_10, Y_pred_10, acc_10, "+10", axes[0]),
    (Y_true_n5, Y_pred_n5, acc_n5, "-5",  axes[1]),
]:
    cm = confusion_matrix(Y_t, Y_p)
    # Normalize rows (per-class recall)
    cm_norm = cm.astype(float) / cm.sum(axis=1, keepdims=True)
    
    disp = ConfusionMatrixDisplay(confusion_matrix=cm_norm,
                                   display_labels=modulation_list)
    disp.plot(ax=ax, colorbar=True, cmap='Blues', xticks_rotation=45,
              values_format='.2f')
    ax.set_title(f'Confusion Matrix SNR = {snr_val} dB (acc={acc:.3f})')

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

# Find the two most confused class pairs at +10 dB
cm_10 = confusion_matrix(Y_true_10, Y_pred_10)
np.fill_diagonal(cm_10, 0)  # zero out diagonal (correct predictions)

for _ in range(2):
    row, col = np.unravel_index(np.argmax(cm_10), cm_10.shape)
    confused_count = cm_10[row, col]
    print(f"Most confused pair: {modulation_list[row]} → predicted as {modulation_list[col]}: {confused_count} errors")
    cm_10[row, col] = 0  # remove so next iteration finds the second pair

Part C deliverable: The confusion matrix plot (both SNRs), the accuracy figures at each SNR, and identification of the two most confused class pairs at +10 dB. For each pair, explain in 1-2 sentences why those two modulations would be confused (what property do they share?).

Part D: Real Signal Test (3 points)

Apply the trained classifier to one segment of the unknown_signal.cf32 from Lab 7:

#!/usr/bin/env python3
"""Lab 11 Part D: Apply classifier to Lab 7 unknown signal."""
import numpy as np
import torch

# Load the Lab 7 unknown signal (adjust path if needed)
LAB7_FILE = '../labs/unknown_signal.cf32'  # adjust to your path
iq_unknown = np.fromfile(LAB7_FILE, dtype=np.complex64)
fs_unknown = 2.4e6

# Extract one 128-sample frame from the most active part
# Find the peak-energy window
window = 128
n_windows = len(iq_unknown) // window
energies = np.array([np.mean(np.abs(iq_unknown[i*window:(i+1)*window])**2) for i in range(n_windows)])
best_window = np.argmax(energies)

frame = iq_unknown[best_window * window : (best_window + 1) * window]
# Normalize
frame = frame / (np.max(np.abs(frame)) + 1e-10)

# Convert to (2, 128) float32
x_input = np.stack([frame.real, frame.imag]).astype(np.float32)
x_tensor = torch.from_numpy(x_input).unsqueeze(0)  # (1, 2, 128)

# Run classifier
model = ModulationCNN(n_classes=len(modulation_list))
model.load_state_dict(torch.load('lab11/models/modulation_cnn.pt', map_location='cpu'))
model.eval()

with torch.no_grad():
    logits = model(x_tensor)
    probs = torch.softmax(logits, dim=1).squeeze().numpy()

# Top-3 predictions
top3_idx = np.argsort(probs)[::-1][:3]
print("Top-3 classification results for Lab 7 unknown_signal.cf32:")
for i, idx in enumerate(top3_idx):
    print(f"  {i+1}. {modulation_list[idx]:12s}: {probs[idx]:.3f} ({probs[idx]*100:.1f}%)")

print(f"\nLab 7 Stage 2 manual hypothesis: [paste your Stage 2 result here]")
print("Comparison: [agree / disagree -- explain why if they differ]")

Part D deliverable: The top-3 class probabilities from the CNN classifier, your Lab 7 Stage 2 manual hypothesis, and a brief comparison (100 words): do they agree? If they disagree, which do you trust more and why?

Lab Report

Create lab-11-report.md with:

1. Part A: dataset distribution plot + balanced/imbalance assessment + SNR discussion
2. Part B: learning curves + final validation accuracy + SNR=+10dB accuracy
3. Part C: confusion matrix plot + accuracy at both SNRs + most confused pairs with explanations
4. Part D: top-3 predictions + comparison with Lab 7 Stage 2

Grading