In [ ]:
# Copyright 2019 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib notebook
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from IPython.display import clear_output, display, HTML
import scipy.special as special # erfc
from scipy import signal # convolution
from ient_nb.ient_plots import *
from ient_nb.ient_signals import *
Hilfsfunktionen¶
Faltung sowie Konvertierung String zu Bits und umgekehrt
In [ ]:
def convolution(s, h):
# Convolve s and h numerically
g = signal.convolve(s, h, mode='full')*deltat; g = g[0:len(s)];
return g
def str2bits(s):
return np.array([int(x) for x in list(''.join(format(ord(i),'b').zfill(8) for i in s))])
def bits2str(b):
try: # decode string
s=''.join(np.array([chr(i) for i in np.packbits(b)]))
except: # decoding failed
s = '-'
return s
Sender¶
Streng geheime Nachricht wird in binäres Signal $a_n$ codiert:
In [ ]:
a = str2bits('Institut für Nachrichtentechnik')
La = len(a) # Anzahl an Bits
a[0:10] # Ausgabe der ersten 10 Werte
Parameter und Achsen
In [ ]:
# Parameter
fs = 1000 # Samplingrate
T=1 # Breite des Trägersignals in Sekunden
Tmax = La*T # Gesamtlänge des Sendesignals
# Achsen
(t,deltat) = np.linspace(0, La*T, La*T*fs, retstep=True) # Zeitachse
nT = np.linspace(1,La*T,La*T) # nT-Achse
nT_idxs = (nT*fs).astype(int); nT_idxs[-1] = nT_idxs[-1]-1 # nT in Samples
Trägersignal $s(t)=\mathrm{rect}\left(\frac{t}{T}-\frac{1}{2}\right)$
In [ ]:
# Trägersignal
s = lambda t: rect(t/T-0.5)
#s = lambda t: rect(2*(t-0.25))-rect(2*(t-0.75))
# Sendeenergie wenn a_n=1
Es = T # s(t) Rechteck der Breite T!
# Plot Trägersignal s(t)
fig,ax = plt.subplots(1,1); plt.plot(t/T, s(t), 'rwth');
ax.set_xlabel(r'$\rightarrow t/T$'); ax.set_ylabel(r'$\uparrow s(t)$');
ax.set_xlim([0,2]); ax.set_ylim([-0.1,1.1]); ient_axis(ax);
Sendesignal $m(t)=\sum_{n=-\infty}^\infty a_n s(t-nT)$
In [ ]:
m = np.zeros(len(t)) # Sendesignal
for l in range(len(a)):
m = m + a[l]*s(t-l*T)
# Plot
fig,ax = plt.subplots(1,1); plt.plot(t/T, m, 'rwth');
ax.set_xlabel(r'$\rightarrow t/T$'); ax.set_ylabel(r'$\uparrow m(t)$');
ax.set_xlim([-0.1,11.9]); ax.set_ylim(0,1.09); ient_axis(ax);
Kanal¶
$n(t)$ Gauß-verteiltes additives weißes Rauschen mit Leistungsdichte $N_0$
In [ ]:
def simulate_channel(N0, l=len(m)):
# Konvertierung Leistungsdichte in Leistung (im Diskreten). Einheit Leistungsdichte: Ws=W/Hz, Einheit Leistung: W
deltaf = 1/fs
sigma_n_squared = N0/deltaf # Varianz von n(t)
# n(t) Werte sind Gauß-verteilt
n = np.random.normal(0, np.sqrt(sigma_n_squared), l)
return n
Augenblicksleistung des gefilterten Signals $n_\mathrm{e}(t)=n(t)\ast g(t)$: $$N=\mathcal{E}\{n_\mathrm{e}^2(T)\} = N_0 \int\limits_{-\infty}^\infty |h(t)|^2 \mathrm{d}t$$
In [ ]:
def calculate_N(N0,h):
return N0*sum(h**2)*deltat
def calculate_N0(N,h):
return N/(sum(h**2)*deltat)
Empfänger¶
Filter $h(t) = s(T-t)$ mit $s(t)=\mathrm{rect}\left(\frac{t}{T}-\frac{1}{2}\right)$ und $$y(t)=g(t)+n_\mathrm{e}(t)=m(t)\ast h(t) + n(t) \ast h(t)$$
Zunächst wird $g(t) = m(t) \ast h(t)$ statisch berechnet:
In [ ]:
# Korrelationsfilter
h = s(T-t)
g = convolution(m,h); # Faltung: g(t) = m(t) * h(t)
# Plot
fig,ax = plt.subplots(1,1)
plt.plot(t/T, m, 'rwth'); plt.plot(t/T, g, 'grun')
ax.set_ylabel(r'$\uparrow m(t),g(t)$'); ax.set_xlabel(r'$\rightarrow t/T$');
ax.set_xlim([-0.1,11.9]); ax.set_ylim(0,1.09); plt.legend({r'$m(t)$',r'$g(t)$'},loc='upper right'); ient_axis(ax);
Da später die Störleistungsdichte $N_0$ des Kanals variiert werden soll, wird die Faltung $n_\mathrm{e}(t) = n(t) \ast h(t)$ in einer Funktion berechnet:
In [ ]:
def simulate_receiver(n,h,g,C):
# Faltung mit Filter h(t)
# eigentlich y(t) = [m(t)+n(t)]*h(t); allerdings ist m(t) determiniert und g(t) vorher berechnet worden
ne = convolution(n, h);
y = g + ne
# Abtaster
yT = y[nT_idxs]; gT = np.round(g[nT_idxs]); neT = ne[nT_idxs];
# Entscheidungsstufe
a_dec = np.array([0 if yTn < C else 1 for yTn in yT]) # Entscheidung: y(nT) > C: an=1, sonst: an=0
str_dec = bits2str(a_dec) # Decodiere String aus Bits
return (str_dec, a_dec, y, ne, yT,gT,neT)
Sendeleistung und -energie:¶
- Augenblicksleistung $S_\mathrm{a} = \mathcal{E}\{g^2(T)\}=g^2(T)$
- Energie (wenn $a_n=1$, ansonsten $E_s=0$!): $E_s = \int\limits_{-\infty}^{\infty}|s(t)|^2\mathrm{d}t$
- Energie pro gesendetem Bit: $E_b = \frac{E_s}{2}$
In [ ]:
Sa = 1 # Augenblickssendeleistung: Sa = g(T)^2=1 wenn a_n=1
Eb = Es/2 # Energie pro gesendetem Bit
C = np.sqrt(Sa)/2 # Entscheidungsschwelle
# Die gleichen Werte, diesmal nur mit numerischer Ungenauigkeit
gT = g[nT_idxs];
Sa_measured = np.mean(gT[a==1]**2)
Es_measured = sum(np.abs(s(t))**2)*deltat # Integral -> Summe im Diskreten; dt -> deltat (Abstand zwischen zwei Samplen)
Eb_measured = Es_measured/2
C_measured = np.sqrt(Sa_measured)/2
Statistiken¶
Bitfehlerwahrscheinlichkeit $$ P_b = \frac{1}{2}\mathrm{erfc}\left(\sqrt{\frac{E_b}{4N_0}}\right) $$
In [ ]:
def calculate_prob(Eb,N0):
Pb = 0.5*special.erfc(np.sqrt(Eb/(4*N0)))
return Pb
Interaktive Demo: Unipolare Binärübertragung¶
In [ ]:
(fig,axs) = plt.subplots(2, 2, figsize=(8,6))
interact_manual.opts['manual_name'] = 'Update'
@interact_manual(N0 = widgets.FloatSlider(min=0.001, max=.1, step=0.001, value=0.01, description='$N_0$'))
def update_signals(N0=0.01):
# Simulation
n = simulate_channel(N0) # Simuliere Kanal
str_dec, a_dec, y, ne, yT,gT,neT = simulate_receiver(n,h,g,C) # Simuliere Empfänger
# Plot
if not axs[0,1].lines: # Plot
ax=axs[0,0]; ax.plot(t/T, m, 'rwth');
ax.set_ylabel(r'$\uparrow m(t)$'); ax.set_xlabel(r'$\rightarrow t/T$'); ax.set_ylim(0,1.09);
ax=axs[0,1]; ax.plot(t/T, m+n, 'rwth');
ax.set_ylabel(r'$\uparrow m(t)+n(t)$'); ax.set_xlabel(r'$\rightarrow t/T$');
ax=axs[1,0]; ax.plot(t/T, g, 'rwth');
ax.plot(nT/T, gT, color='rwth', marker='o', lw=0)
ax.set_ylabel(r'$\uparrow g(t)$'); ax.set_xlabel(r'$\rightarrow t/T$'); ax.set_ylim(-1,2.09);
ax=axs[1,1]; ax.plot(t/T, y, 'rwth');
ax.plot(nT/T, yT, color='rwth', marker='o', lw=0)
ax.set_ylabel(r'$\uparrow y(t)$'); ax.set_xlabel(r'$\rightarrow t/T$'); ax.set_ylim(-1,2.09);
ax.plot(ax.get_xlim(), [C, C], 'k--',lw=1); ax.text(0, C+0.1, r'$C$', fontsize=12, horizontalalignment='left',verticalalignment='bottom',backgroundcolor='w');
[ax.set_xlim([-0.1,11.9]) and ient_axis(ax) for ax in axs.flatten()];
else: # update lines
ax = axs[0,1]; ax.lines[0].set_ydata(m+n)
ax = axs[1,1]; ax.lines[0].set_ydata(y); ax.lines[1].set_ydata(yT);
# Leistungen, S/N Verhältnis
N = calculate_N(N0,h); SN = Sa/N;
Pb = calculate_prob(Eb, N0);
N_measured = np.mean(neT**2); SN_measured = Sa_measured/N_measured;
Pb_measured = np.mean(a_dec != a);
# Output
display(HTML(
"""<table><tr><td>TX a</td>{}</tr><tr><td>RX a</td>{}</tr></table><br />Empfangen: {}<br />
Sa/N={:.2f} dB (berechnet: {:.2f} dB), log10(Pb)={:.4f} (berechnet: {:.4f})""".format(
'<td>{}</td>'.format('</td><td>'.join(str(tmp) for tmp in a)),
'<td>{}</td>'.format('</td><td>'.join(str(tmp) for tmp in a_dec)), str_dec,
10*np.log10(SN_measured), 10*np.log10(SN), np.log10(Pb_measured+np.finfo(float).eps), np.log10(Pb),
)
))
Wahrscheinlichkeitsdichtefunktion¶
Mit Verteilungsdichtefunktion von $n_\mathrm{e}(nT)$ $$p_{n_\mathrm{e}}(x) = \frac{1}{\sqrt{2\pi N}} \exp\left[ - \frac{x^2}{2N} \right]$$ und VDF von $g(nT)$ $$p_g(t) = \mathrm{Prob}\{a_n=0\} \cdot \delta(x) + \mathrm{Prob}\{a_n=1\} \cdot \delta\left(x-\sqrt{S_a}\right)$$ folgt VDF von $y(nT) = g(nT) + n_\mathrm{e}(nT)$ $$p_y(x) = p_g(x) \ast p_{n_\mathrm{e}}(x) = \frac{1}{\sqrt{2 \pi N}} \left( \mathrm{Prob}\{a_n=0\} \cdot \exp\left[ -\frac{x^2 }{2N} \right] + \mathrm{Prob}\{a_n=1\} \cdot \exp\left[ -\frac{\left(x-\sqrt{S_\mathrm{a}}\right)^2}{2N} \right] \right)$$
In [ ]:
def calculate_pdf(x,a,N,mean0,mean1):
# a_n = 0
pa0 = np.sum(a==0)/len(a);
py0 = 1/np.sqrt(2*np.pi*N)*np.exp(-(x-mean0)**2/(2*N))
# a_n = 1
pa1 = np.sum(a==1)/len(a);
py1 = 1/np.sqrt(2*np.pi*N)*np.exp(-(x-mean1)**2/(2*N))
# Gesamtwahrscheinlichkeit
py = pa0*py0 + pa1*py1
return py
Demo¶
In [ ]:
(logSN,deltaSN) = np.linspace(-5,25,1000, retstep=True) # logSN axis
(fig_pdf,axs_pdf) = plt.subplots(2, 1, figsize=(8,8))
np.seterr(divide = 'ignore')
@interact(N0 = widgets.FloatSlider(min=0.003, max=.5, step=0.001, value=0.005,
description=r'$N_0$', readout_format='.3f', continuous_update=False))
def update_signals_pdf(N0=0.01):
# Simulation
n = simulate_channel(N0) # Simuliere Kanal
str_dec, a_dec, y, ne, yT,gT,neT = simulate_receiver(n,h,g,C) # Simuliere Empfänger
# Leistungen, S/N Verhältnis
N = calculate_N(N0,h); SN = Sa/N;
Pb = calculate_prob(Eb, N0)
N_measured = np.mean(neT**2); SN_measured = Sa_measured/N_measured;
Pb_measured = np.mean(a_dec != a)
# Verteilungsdichte p_y(x)
py_measured, bins = np.histogram(yT,bins=100,range=(-0.5,1.5),density=True) # gemessen mit Histogramm
x = (bins[:-1] + bins[1:]) / 2
py = calculate_pdf(x,a,N,0,np.sqrt(Sa)) # berechnet
# Plot
if not axs_pdf[0].lines:
ax = axs_pdf[0]; #ax.plot(x, py_measured, 'rwth');
ient_stem(ax, x, py_measured, 'rwth', markerfmt=" ")
ax.plot(x, py, 'grun')
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_y(x)$'); ient_axis(ax);
ax = axs_pdf[1];
ax.semilogy(logSN, 0.5*special.erfc(np.sqrt(10**(logSN/10)/8)), color="rwth")
ax.semilogy(10*np.log10(SN_measured), Pb_measured , marker='o', markersize=5, color="rot");
ax.set_xlabel(r'$\rightarrow E_b/N_0$ [dB]'); ax.set_ylabel(r'$\rightarrow P_{b}$');
ax.text(ax.get_xlim()[1]-1, .24, r'$0{,}5$', fontsize=12,
horizontalalignment='right',verticalalignment='top',backgroundcolor='w');
ax.plot(ax.get_xlim(), [0.5, 0.5], 'k--',lw=1); ax.grid(); ax.set_ylim([1e-10,2]);
else: # update lines
ax = axs_pdf[0];
ient_stem_set_data(ax.containers[0], x, py_measured)
ax.lines[-1].set_data(x, py);
ax = axs_pdf[1]; ax.lines[1].set_xdata(10*np.log10(SN_measured+eps)); ax.lines[1].set_ydata(Pb_measured);
# Output
display(HTML('gemessen: log10(Pb)={:.4f} (berechnet: {:.4f})'.format(
np.log10(Pb_measured),np.log10(Pb))))
This notebook is provided as Open Educational Resource (OER). Feel free to use the notebook for your own purposes. The code is licensed under the MIT license.
Please attribute the work as follows: Christian Rohlfing, Übungsbeispiele zur Vorlesung "Informationsübertragung", gehalten von Jens-Rainer Ohm, 2019, Institut für Nachrichtentechnik, RWTH Aachen University.