Exercises
ex-sp-ch09-01
EasyUse scipy.stats.norm to compute for
. Verify by drawing samples and
computing the empirical probability.
Use norm.sf(2) + norm.cdf(-2) or 2 * norm.sf(2).
Implementation
from scipy.stats import norm
import numpy as np
p_theory = 2 * norm.sf(2)
rng = np.random.default_rng(42)
samples = rng.standard_normal(1_000_000)
p_empirical = np.mean(np.abs(samples) > 2)
print(f"Theory: {p_theory:.6f}")
print(f"Empirical: {p_empirical:.6f}")
ex-sp-ch09-02
EasyGenerate 10000 pairs of correlated Gaussian random variables with using the Cholesky decomposition. Compute the sample correlation and verify it matches the target.
Build and use np.linalg.cholesky.
Implementation
import numpy as np
rho = 0.6
R = np.array([[1, rho], [rho, 1]])
L = np.linalg.cholesky(R)
rng = np.random.default_rng(42)
z = rng.standard_normal((2, 10000))
x = L @ z
print(f"Sample correlation: {np.corrcoef(x[0], x[1])[0,1]:.4f}")
print(f"Target: {rho}")
ex-sp-ch09-03
EasyWrite a BPSK BER simulation for a single dB point with bits. Compare the simulated BER to the theoretical .
BPSK: bits , add AWGN with variance .
Implementation
import numpy as np
from scipy.special import erfc
rng = np.random.default_rng(42)
ebno_db = 8.0
ebno = 10 ** (ebno_db / 10)
N = 1_000_000
bits = rng.integers(0, 2, size=N)
x = 2.0 * bits - 1.0
noise = rng.standard_normal(N) / np.sqrt(2 * ebno)
y = x + noise
errors = np.sum(bits != (y > 0).astype(int))
ber = errors / N
ber_theory = 0.5 * erfc(np.sqrt(ebno))
print(f"Simulated: {ber:.4e}, Theory: {ber_theory:.4e}")
ex-sp-ch09-04
EasyGenerate 50000 Rayleigh fading channel coefficients . Plot a histogram of and overlay the theoretical Rayleigh PDF. Verify with a KS test.
Use scipy.stats.rayleigh for the theoretical distribution.
rayleigh.pdf(x) with default scale=1 gives the Rayleigh PDF.
Implementation
import numpy as np
from scipy.stats import rayleigh, kstest
rng = np.random.default_rng(42)
h = (rng.standard_normal(50000)
+ 1j * rng.standard_normal(50000)) / np.sqrt(2)
envelope = np.abs(h)
D, p = kstest(envelope, 'rayleigh', args=(0, 1/np.sqrt(2)))
print(f"KS stat: {D:.4f}, p-value: {p:.4f}")
ex-sp-ch09-05
EasyCompute the 95% confidence interval for with bits using both the normal approximation and the Clopper-Pearson exact method.
For Clopper-Pearson, use scipy.stats.beta.ppf.
Implementation
import numpy as np
from scipy.stats import beta, norm
p_hat = 0.0023
N = 100000
k = int(p_hat * N) # 230 errors
# Normal approximation
z = norm.ppf(0.975)
ci_normal = (p_hat - z*np.sqrt(p_hat*(1-p_hat)/N),
p_hat + z*np.sqrt(p_hat*(1-p_hat)/N))
# Clopper-Pearson
ci_cp = (beta.ppf(0.025, k, N-k+1),
beta.ppf(0.975, k+1, N-k))
print(f"Normal: [{ci_normal[0]:.6f}, {ci_normal[1]:.6f}]")
print(f"CP: [{ci_cp[0]:.6f}, {ci_cp[1]:.6f}]")
ex-sp-ch09-06
MediumFit Rayleigh, Rice, and Nakagami distributions to synthetic Rice-distributed data ( dB). Use the KS test to rank the fits. Plot all three fitted PDFs against the histogram.
Use scipy.stats.rice, rayleigh, nakagami with .fit() and kstest().
Implementation
from scipy.stats import rice, rayleigh, nakagami, kstest
import numpy as np
K_db = 3.0
K = 10 ** (K_db / 10)
nu = np.sqrt(K / (K+1))
data = rice.rvs(nu, scale=np.sqrt(1/(K+1)), size=10000,
random_state=42)
for name, dist in [("Rayleigh", rayleigh),
("Rice", rice), ("Nakagami", nakagami)]:
params = dist.fit(data)
D, p = kstest(data, dist.cdf, args=params)
print(f"{name:10s}: D={D:.4f}, p={p:.4f}")
ex-sp-ch09-07
MediumImplement a two-sample t-test to compare the mean BER of two MIMO detectors (ZF and MMSE). Run 30 independent trials at dB for each detector and test for a significant difference at .
Implementation
from scipy.stats import ttest_ind
import numpy as np
rng = np.random.default_rng(42)
ber_zf = 2e-3 + 5e-4 * rng.standard_normal(30)
ber_mmse = 1.5e-3 + 5e-4 * rng.standard_normal(30)
stat, p = ttest_ind(ber_zf, ber_mmse)
print(f"t={stat:.3f}, p={p:.4f}")
print(f"Significant: {p < 0.05}")
ex-sp-ch09-08
MediumWrite a BER simulation with error-count stopping (minimum 100 errors) for BPSK over AWGN. Run for dB. Report BER, error count, total bits, and 95% CI for each point.
Implementation
import numpy as np
def simulate_until_errors(ebno_db, min_errors=100,
block_size=100_000):
rng = np.random.default_rng(42)
ebno = 10 ** (ebno_db / 10)
total_err, total_bits = 0, 0
while total_err < min_errors and total_bits < 1e9:
bits = rng.integers(0, 2, size=block_size)
x = 2.0*bits - 1.0
y = x + rng.standard_normal(block_size)/np.sqrt(2*ebno)
total_err += np.sum(bits != (y>0).astype(int))
total_bits += block_size
ber = total_err / total_bits
ci = 1.96 * np.sqrt(ber*(1-ber)/total_bits)
return ber, total_err, total_bits, ci
for ebn0 in [0, 2, 4, 6, 8, 10]:
ber, ne, nb, ci = simulate_until_errors(ebn0)
print(f"Eb/N0={ebn0:2d} dB: BER={ber:.3e} "
f"[+/-{ci:.1e}], errors={ne}, bits={nb}")
ex-sp-ch09-09
MediumGenerate a Kronecker MIMO channel with , . Generate 5000 realizations and verify the vectorized channel covariance matches .
Implementation
import numpy as np
rho_t, rho_r = 0.3, 0.8
Rt = np.array([[1, rho_t], [rho_t, 1]])
Rr = np.array([[1, rho_r], [rho_r, 1]])
Lt = np.linalg.cholesky(Rt)
Lr = np.linalg.cholesky(Rr)
rng = np.random.default_rng(42)
N = 5000
h_vec = np.zeros((4, N), dtype=complex)
for n in range(N):
Hw = (rng.standard_normal((2,2))
+ 1j*rng.standard_normal((2,2)))/np.sqrt(2)
H = Lr @ Hw @ Lt.T
h_vec[:, n] = H.ravel(order='F')
R_hat = (h_vec @ h_vec.conj().T) / N
R_theory = np.kron(Rt.conj(), Rr)
print(f"Max error: {np.max(np.abs(R_hat - R_theory)):.4f}")
ex-sp-ch09-10
MediumUse antithetic variates to estimate the BER of BPSK at dB with samples. Compare the variance of the estimator with crude Monte Carlo over 100 independent runs.
Implementation
import numpy as np
from scipy.stats import norm
ebno = 10 ** (6/10)
d = np.sqrt(2 * ebno)
ss = np.random.SeedSequence(42)
seeds = ss.spawn(100)
N = 100_000
mc_bers, av_bers = [], []
for seed in seeds:
rng = np.random.default_rng(seed)
u = rng.uniform(size=N)
n1 = norm.ppf(u)
mc_bers.append(np.mean(n1 > d))
n2 = norm.ppf(1 - u)
av_bers.append(np.mean(0.5*((n1>d)+(n2>d))))
print(f"MC var: {np.var(mc_bers):.2e}")
print(f"AV var: {np.var(av_bers):.2e}")
print(f"Ratio: {np.var(mc_bers)/np.var(av_bers):.2f}x")
ex-sp-ch09-11
MediumParallelize a BPSK BER simulation across 4 SNR points using
concurrent.futures.ProcessPoolExecutor. Use SeedSequence.spawn()
for independent random streams.
Implementation
from concurrent.futures import ProcessPoolExecutor
import numpy as np
def worker(args):
ebno_db, seed = args
rng = np.random.default_rng(seed)
ebno = 10**(ebno_db/10)
N = 1_000_000
bits = rng.integers(0, 2, size=N)
x = 2.0*bits - 1.0
y = x + rng.standard_normal(N)/np.sqrt(2*ebno)
errors = np.sum(bits != (y>0).astype(int))
return ebno_db, errors/N, errors
ss = np.random.SeedSequence(42)
snrs = [4, 6, 8, 10]
seeds = ss.spawn(len(snrs))
with ProcessPoolExecutor(max_workers=4) as pool:
results = list(pool.map(worker, zip(snrs, seeds)))
for snr, ber, ne in results:
print(f"Eb/N0={snr} dB: BER={ber:.4e} ({ne} errors)")
ex-sp-ch09-12
HardImplement importance sampling for BPSK BER estimation at dB. Use exponential tilting with shift (decision boundary). With only IS samples, compare the accuracy to crude MC with samples.
Implementation
import numpy as np
from scipy.special import erfc
ebno_db = 14.0
ebno = 10 ** (ebno_db / 10)
d = np.sqrt(2 * ebno)
ber_theory = 0.5 * erfc(np.sqrt(ebno))
rng = np.random.default_rng(42)
# Crude MC
N_mc = 10_000_000
noise_mc = rng.standard_normal(N_mc)
ber_mc = np.mean(noise_mc > d)
# IS
N_is = 1000
noise_is = d + rng.standard_normal(N_is)
w = np.exp(-d * noise_is + d**2 / 2)
ber_is = np.mean((noise_is > d) * w)
print(f"Theory: {ber_theory:.6e}")
print(f"Crude MC: {ber_mc:.6e} (N={N_mc})")
print(f"IS: {ber_is:.6e} (N={N_is})")
ex-sp-ch09-13
HardGenerate a time-varying Rayleigh channel using Jakes' model for a mobile at 120 km/h, carrier frequency 3.5 GHz, sampled at 10 kHz. Plot the channel magnitude (dB) vs time and verify the level crossing rate matches theory: where is the normalized threshold.
ex-sp-ch09-14
HardSimulate BPSK over a Rayleigh fading channel for to dB. Verify the simulated BER matches the theoretical where . Use error-count stopping with errors.
ex-sp-ch09-15
HardImplement control variates for estimating the average capacity of a Rayleigh fading channel. Use as the control variate (known mean = 1) and quantify the variance reduction over 1000 independent runs.
ex-sp-ch09-16
HardBuild a 3GPP TDL-A channel model with 7 taps. Generate 10000 channel impulse response realizations. Compute the power delay profile and verify the per-tap powers match the standard values. Compute the RMS delay spread.
ex-sp-ch09-17
ChallengeDesign an adaptive importance sampling scheme for estimating the BER of 16-QAM over AWGN. The challenge: 16-QAM has multiple decision boundaries, so a single shift is suboptimal. Implement a mixture proposal that shifts toward each nearest boundary. Compare convergence with crude MC at dB.
ex-sp-ch09-18
ChallengeBuild a complete MIMO link-level simulator with: (a) Kronecker channel model with , (b) ZF and MMSE detection, (c) BPSK modulation, (d) error-count stopping with . Plot BER vs SNR for both detectors on the same graph. Include 95% confidence intervals as error bars.