Salta al contenuto
Note per Studenti Calcoli di comunicazioni in Python - Q, dB, quantizzatore, Huffman e link budget

Calcoli di comunicazioni in Python - Q, dB, quantizzatore, Huffman e link budget

In questa pagina 6

Tutti gli esercizi svolti di questo corso sono stati verificati con le funzioni sotto: è il modo più rapido per controllare un calcolo prima di consegnarlo (le formule da usare a mano sono in Formulario - fondamenti di comunicazioniTutte le formule del corso in una pagina, in ordine di capitolo: decibel e probabilità ($Q$), quantizzazione ($\Delta$, $\Lambda_q$), entropia e codifica, rumore e figura di rumore, link budget (Friis), spazio dei segnali, probabilità di errore (binaria, $M$-aria, PAM, QAM, PSK), codifica di canale, capacità, accesso al mezzo, ARQ, code. Serve per costruire il foglio A4 consentito all'esame; ogni gruppo rimanda alla nota con le ipotesi di validità, che contano quanto la formula.Formulario - fondamenti di comunicazioni →). Il codice è stato eseguito e i numeri riportati nei commenti sono quelli ottenuti.

Funzioni di base

Q è la funzione coda della gaussiana standard (Variabili aleatorie e vettori aleatori per le comunicazioniUna variabile aleatoria è descritta da una PMD (discreta) o da una PDF (continua) e dalla funzione di ripartizione; per la gaussiana $P[x>a]=Q\left(\frac{a-m}{\sigma}\right)$ con la funzione coda $Q$. Media $m_x$, varianza $\sigma_x^2$ e potenza statistica $M_x=\sigma_x^2+m_x^2$. Per un vettore aleatorio contano l'indipendenza, le probabilità condizionate (totali e di Bayes) e la correlazione; combinazioni lineari di gaussiane sono gaussiane, e gaussiane scorrelate sono indipendenti.Variabili aleatorie e vettori aleatori per le comunicazioni →): è la survival functionprobabilità che la variabile superi un valore, cioè 1 meno la funzione di ripartizione norm.sf; la sua inversa è norm.isf.

python
import heapq
import numpy as np
from scipy.stats import norm

Q = norm.sf          # Q(x) = P[N(0,1) > x]
Qinv = norm.isf      # Q^{-1}(p)
db = lambda x: 10 * np.log10(x)
undb = lambda x: 10 ** (np.asarray(x) / 10)

print(Q(3), Qinv(1e-6))      # 0.00135, 4.753

Quantizzatore gaussiano a media nulla

Dato σ\sigma, i bit bb e la probabilità di saturazione, il range dinamico è Vsat=σQ−1(Psat/2)V_{sat}=\sigma Q^{-1}(P_{sat}/2) e l'SNR è σ2/(Δ2/12)\sigma^2/(\Delta^2/12) (SNR di quantizzazione e progetto del quantizzatoreL'SNR di quantizzazione è $\Lambda_q=\frac{M_a}{M_e}$. Con errore granulare uniforme e saturazione trascurabile vale $\Lambda_q=\frac{\sigma^2}{\Delta^2/12}=3\frac{\sigma^2}{V_{sat}^2},2^{2b}$, cioè $[\Lambda_q]{dB}=6{,}02,b+4{,}77+20\log{10}\frac\sigma{V_{sat}}$: ogni bit in più dà $+6$ dB. Per progettare: $V_{sat}$ dalla probabilità di saturazione ($V_{sat}=\sigma,Q^{-1}\left(\frac{P_{sat}}2\right)$ per un gaussiano), poi $b$ dall'SNR richiesto, arrotondando per eccesso.SNR di quantizzazione e progetto del quantizzatore →).

python
def quantizzatore_gauss(sigma, b, psat):
    vsat = sigma * Qinv(psat / 2)
    delta = 2 * vsat / 2 ** b
    snr = db(sigma**2 / (delta**2 / 12))
    return vsat, delta, snr

print([round(float(x), 3) for x in quantizzatore_gauss(np.sqrt(2), 3, 2e-4)])
# [5.259, 1.315, 11.424]   (mV, mV, dB) per sigma^2 = 2 mV^2, b = 3, Psat = 2e-4

Entropia, Huffman e Kraft

python
def entropia(p):
    p = np.asarray([x for x in p if x > 0])
    return float(-(p * np.log2(p)).sum())

def huffman_lunghezze(p):
    heap = [(q, i, [i]) for i, q in enumerate(p)]
    heapq.heapify(heap)
    lun, n = [0] * len(p), len(p)
    while len(heap) > 1:                       # si fondono sempre i due meno probabili
        q1, _, a = heapq.heappop(heap)
        q2, _, b = heapq.heappop(heap)
        for i in a + b:
            lun[i] += 1                        # ogni fusione aggiunge un bit ai simboli coinvolti
        heapq.heappush(heap, (q1 + q2, n, a + b))
        n += 1
    return lun

p = [.25, .25, .25, .10, .05, .05, .04, .01]
lun = huffman_lunghezze(p)
Ly = sum(l * q for l, q in zip(lun, p))
print(lun, round(Ly, 3), round(entropia(p), 3), round(entropia(p) / Ly, 3))
# [2, 2, 2, 4, 4, 4, 5, 5]  Ly = 2.55  H = 2.517  efficienza 0.987
print("Kraft:", sum(2.0 ** -l for l in (2, 2, 3, 3, 3, 3, 4, 4)))   # 1.125 > 1: codice impossibile

Il codice restituisce le lunghezze (che sono quelle che servono per LyL_y e per l'efficienza, Codici di Shannon-Fano e di HuffmanIn un codice ottimo le parole più probabili non sono più lunghe di quelle meno probabili e le due parole più lunghe differiscono solo per l'ultimo simbolo. Shannon-Fano costruisce l'albero dall'alto dividendo ripetutamente i simboli in due gruppi di probabilità quasi uguali; Huffman lo costruisce dal basso unendo ogni volta i due simboli meno probabili ed è sempre ottimo tra i codici a prefisso. La lunghezza media $L_y$ è la somma delle probabilità dei nodi uniti, l'efficienza è $\eta=\frac{H}{L_y}$.Codici di Shannon-Fano e di Huffman →); la scelta tra pari probabilità può dare lunghezze diverse da quelle di un'altra costruzione a mano, ma con lo stesso LyL_y.

python
def attenuazione_friis(d_km, f_mhz, gtx_db=0, grx_db=0):
    return 32.4 + 20 * np.log10(d_km) + 20 * np.log10(f_mhz) - gtx_db - grx_db

def snr_rif_db(ptx_dbm, a_db, F_db, B_hz):
    return ptx_dbm - a_db + 114 - F_db - 10 * np.log10(B_hz / 1e6)

def pbit_qam(M, gamma):
    return 4 / np.log2(M) * (1 - 1 / np.sqrt(M)) * Q(np.sqrt(3 * gamma / (M - 1)))

def pe_qam(M, gamma):
    return 4 * (1 - 1 / np.sqrt(M)) * Q(np.sqrt(3 * gamma / (M - 1)))

G = snr_rif_db(10, 90, 17, 60e3)      # radio locale: 10 dBm, 90 dB, F = 17 dB, B = 60 kHz
print(round(G, 2))                    # 29.22 dB
for M in (16, 64, 256, 1024):
    print(M, f"{pbit_qam(M, undb(G)):.2e}")
# 16: 1.22e-38   64: 8.30e-11   256: 4.03e-04   1024: 2.28e-02
gmin = db(Qinv(1e-3) ** 2)            # 4-QAM, Pbit = 1e-3
print(round(gmin, 2), round(gmin + 90 - 114 + 17 + 10 * np.log10(0.06), 2))
# 9.8 dB  e  Ptx minima = -9.42 dBm

(Le formule sono quelle di QAM - modulazione di ampiezza in quadraturaNella QAM si modulano ampiezza e fase di una portante con due ampiezze $\alpha_{m,I},\alpha_{m,Q}$ (simbolo complesso $\alpha_m=\alpha_{m,I}+j\alpha_{m,Q}$): $s_m(t)=\operatorname{Re}\left[\alpha_mh(t)e^{j2\pi f_0t}\right]$. È un segnale in banda passante con base a due dimensioni ($\cos$ e $\sin$ per l'impulso) e punti $\sqrt{\frac{E_h}2}\left[\alpha_{m,I},\alpha_{m,Q}\right]$. Per $M=L^2$: $d_{min}=\sqrt{2E_h}$, $E_s=E_h\frac{M-1}3$ e $P[E]=1-\left[1-2\left(1-\frac1{\sqrt M}\right)Q\left(\sqrt{\frac{E_h}{N_0}}\right)\right]^2\approx4\left(1-\frac1{\sqrt M}\right)Q\left(\sqrt{\frac3{M-1}\frac{E_s}{N_0}}\right)$.QAM - modulazione di ampiezza in quadratura → e Progetto di un collegamento - link budget e modulazioneEsercizio tipico d'esame: dati potenza, guadagni, frequenza e banda, si calcola (1) l'attenuazione (spazio libero o cavo) e quindi la distanza massima, (2) l'SNR di riferimento $\Gamma_{dB}=P_{tx,dBm}-a_{ch}+114-F_{dB}-10\log_{10}B_{MHz}$, (3) la cardinalità massima $M$ della QAM (o PSK) tale che $P\le P_{target}$, (4) la potenza minima per una modulazione fissata invertendo $P(\Gamma)$: $\Gamma_{min}=\left[Q^{-1}\right]^2$ per la 4-QAM e poi $P_{tx,min}=\Gamma_{min}+a_{ch}-114+F+10\log_{10}B_{MHz}$.Progetto di un collegamento - link budget e modulazione →.) Per la PSK basta sostituire pbit_psk(M, g) = 2/np.log2(M) * Q(np.sqrt(2*g)*np.sin(np.pi/M)) (per M>2M>2).

Simulazione Monte Carlo

Per essere sicuri che la formula sia giusta (e per vedere quanto l'approssimazione P[E]≈4(1−1M)QP[E]\approx4(1-\frac1{\sqrt M})Q sbaglia a SNR basso) si simula la 16-QAM con rumore AWGN e decisione a minima distanza. Con N0=1N_0=1 le componenti del rumore hanno σI2=N02\sigma_I^2=\frac{N_0}2 e i punti stanno a ±Eh/2⋅{1,3}\pm\sqrt{E_h/2}\cdot\{1,3\} (QAM - modulazione di ampiezza in quadraturaNella QAM si modulano ampiezza e fase di una portante con due ampiezze $\alpha_{m,I},\alpha_{m,Q}$ (simbolo complesso $\alpha_m=\alpha_{m,I}+j\alpha_{m,Q}$): $s_m(t)=\operatorname{Re}\left[\alpha_mh(t)e^{j2\pi f_0t}\right]$. È un segnale in banda passante con base a due dimensioni ($\cos$ e $\sin$ per l'impulso) e punti $\sqrt{\frac{E_h}2}\left[\alpha_{m,I},\alpha_{m,Q}\right]$. Per $M=L^2$: $d_{min}=\sqrt{2E_h}$, $E_s=E_h\frac{M-1}3$ e $P[E]=1-\left[1-2\left(1-\frac1{\sqrt M}\right)Q\left(\sqrt{\frac{E_h}{N_0}}\right)\right]^2\approx4\left(1-\frac1{\sqrt M}\right)Q\left(\sqrt{\frac3{M-1}\frac{E_s}{N_0}}\right)$.QAM - modulazione di ampiezza in quadratura →).

python
rng = np.random.default_rng(0)
gamma = undb(10)                       # Es/N0 = 10 dB
Eh_N0 = 3 * gamma / 15                 # Eh = 3 Es / (M - 1)
amp = np.sqrt(Eh_N0 / 2)               # N0 = 1, punti a amp * (+-1, +-3)
n = 2_000_000
i, q = rng.integers(0, 4, n), rng.integers(0, 4, n)
liv = np.array([-3, -1, 1, 3]) * amp
sig = np.sqrt(0.5)                     # sigma_I^2 = N0/2
ri, rq = liv[i] + rng.normal(0, sig, n), liv[q] + rng.normal(0, sig, n)
dec = lambda r: np.clip(np.floor(r / (2 * amp)) + 2, 0, 3).astype(int)   # soglie a meta' tra i livelli
print(round(1 - np.mean((dec(ri) == i) & (dec(rq) == q)), 3),
      round(pe_qam(16, gamma), 3))
# simulata 0.222, formula approssimata 0.236

Il valore simulato (0,2220{,}222) coincide con la formula esatta 1−[1−2(1−1M)Q(Eh/N0)]2=0,2221-\left[1-2\left(1-\frac1{\sqrt M}\right)Q\left(\sqrt{E_h/N_0}\right)\right]^2=0{,}222 e l'approssimazione di libro (0,2360{,}236) lo sovrastima: a SNR alto la differenza sparisce (Γ=20\Gamma=20 dB: 1,16⋅10−51{,}16\cdot10^{-5} in entrambi i casi).

Errori comuni

  • Usare norm.cdf al posto di norm.sf (si ottiene 1−Q1-Q): per i valori molto piccoli 1−cdf1-\text{cdf} perde precisione, sf no.
  • Dimenticare la conversione undb prima di passare Γ\Gamma alle formule (che lavorano in lineare).
  • Fidarsi di una simulazione con pochi campioni per probabilità minori di 1n\frac1{n}: servono almeno 102…10310^2\ldots10^3 errori.

Versione ripasso

  • Base: Q = norm.sf, Qinv = norm.isf, db = 10*log10, undb = 10**(x/10). Q(3)=1,35⋅10−3Q(3)=1{,}35\cdot10^{-3}, Q−1(10−6)=4,753Q^{-1}(10^{-6})=4{,}753.
  • Quantizzatore gaussiano: Vsat=σQ−1(Psat/2)V_{sat}=\sigma Q^{-1}(P_{sat}/2), Δ=2Vsat2b\Delta=\frac{2V_{sat}}{2^b}, SNR =σ2/(Δ2/12)=\sigma^2/(\Delta^2/12); σ2=2\sigma^2=2 mV², b=3b=3, Psat=2⋅10−4P_{sat}=2\cdot10^{-4}: 5,2595{,}259 mV, 1,3151{,}315 mV, 11,4211{,}42 dB.
  • Huffman con heapq (fonde i due minimi, +1 bit ai simboli coinvolti): 0,253,0,1,0,052,0,04,0,01→[2,2,2,4,4,4,5,5]0{,}25^3,0{,}1,0{,}05^2,0{,}04,0{,}01\to[2,2,2,4,4,4,5,5], Ly=2,55L_y=2{,}55, η=0,987\eta=0{,}987; Kraft 2,2,3,3,3,3,4,42,2,3,3,3,3,4,4: 1,1251{,}125.
  • Link: Friis 32,4+…32{,}4+\dots; ΓdB=Ptx−a+114−F−10log⁡10BMHz\Gamma_{dB}=P_{tx}-a+114-F-10\log_{10}B_{MHz}; pbit_qam, pe_qam; 10 dBm, 90 dB, F=17F=17, 60 kHz: 29,2229{,}22 dB; M=256M=256: 4,03⋅10−44{,}03\cdot10^{-4}; 4-QAM 10−310^{-3}: 9,89{,}8 dB, −9,42-9{,}42 dBm.
  • Monte Carlo 16-QAM a 10 dB: simulata 0,2220{,}222 = formula esatta; approssimata 0,2360{,}236.
  • Errori tipici: cdf al posto di sf; Γ\Gamma in dB nelle formule; troppi pochi campioni.