Skip to content
Snippets Groups Projects
Verified Commit 5b9eb24d authored by Christian Rohlfing's avatar Christian Rohlfing
Browse files

Bugfix in PAM

parent 53e6f647
Branches development
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
IUE Pulsamplitudenmodulation.ipynb
+ 12
− 11
View file @ 5b9eb24d
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from IPython.display import Markdown, Latex
import matplotlib.pyplot as plt
from scipy import signal # convolution,
from scipy.io import wavfile # wavfile
import rwth_nb.plots.mpl_decorations as rwth_plots
import rwth_nb.misc.transforms as rwth_transforms
import rwth_nb.misc.media as rwth_media
import rwth_nb.misc.filters as rwth_filters
from rwth_nb.misc.signals import *
def convolution(s, h):
# Convolve s and h numerically
g = signal.convolve(s, h, mode='same')*deltat; #g = g[0:len(s)];
return g
def ient_ideal_lowpass(s, fg):
# Convolve with impulse response of ideal lowpass
return convolution(s, 2*fg*si(2*np.pi*fg*t))
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Pulsamplitudenmodulation
Zum Starten: Im Menü: Run <span class="fa-chevron-right fa"></span> Run All Cells auswählen.
## Übersicht
![Blockdiagramm](figures/pam_block_diagram.png)
Hier Rauschprozess $n(t)$ (Kanal) vernachlässigt (ansonsten $g(t)\rightarrow y(t)=g(t)+n_e(t)$ mit $n_e(t)=n(t)\ast s(-t)$)
## Nutzsignal
%% Cell type:code id: tags:
``` python
fs = 400000 # interne Abtastrate
fs0, data = wavfile.read('data/krawehl.wav')
f = 0.99*data/np.max(np.abs(data)) # Normalisieren
f = signal.resample(f, int(len(f)/fs0*fs)) # Etwas hochtasten für schönere Plots
#data = np.hstack((data,data))
f0 = 0.99*data/np.max(np.abs(data)) # Normalisieren
f = signal.resample(f0, int(len(f0)/fs0*fs)) # Etwas hochtasten für schönere Plots
(t, deltat) = np.linspace(-len(f)/fs/2, len(f)/fs/2, len(f), retstep=True) # Zeitachse in Sekunden
F, f_ax = rwth_transforms.dft(f, fs) # Fourier-Transformation
# Plot
fig,axs=plt.subplots(2,1, figsize=(8,4));
ax = axs[0]; ax.plot(1000*t, f, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow f(t)$'); rwth_plots.axis(ax);
ax = axs[1]; ax.plot(f_ax, np.abs(F), 'rwth:blue');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow |F(f)|$');
ax.set_xlim([-10E3, 10E3]); rwth_plots.axis(ax);
rwth_media.audio_play(f, fs, r'$f(t)$')
rwth_media.audio_play(f0, fs0, r'$f(t)$')
```
%% Cell type:markdown id: tags:
## Sender
Trägersignal $s(t)=\frac{1}{\sqrt{t_0}}\mathrm{rect}\left(\frac{t}{t_0}\right)$ mit Breite $t_0 = 0{,}5 T$ (oder $t_0 = 1{,}25 T$ für Aufgabe 11) und $T=125\mu\mathrm{s}$
%% Cell type:code id: tags:
``` python
T = 125E-6
t0 = 0.5*T
t0 = 1.25*T
#t0 = 1.25*T # Aufgabe 11: '#' entfernen und alle Cells erneut laufen lassen
# Trägersignal
s = lambda t,t0: 1/np.sqrt(t0)*rect(t/t0)
phi_ss = convolution(s(-t,t0),s(t,t0))
# Plot
fig, axs = plt.subplots(1,2,figsize=(8,4));
ax = axs[0]; ax.plot(t*1E6, s(t,t0), 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [$\mu$s]'); ax.set_ylabel(r'$\uparrow s(t)$');
ax.set_xlim([-2E6*t0, 2E6*t0]); rwth_plots.axis(ax);
rwth_plots.annotate_distance(ax, r'$t_0$', (-t0/2*1E6,10), (t0/2*1E6,10));
ax = axs[1]; ax.plot(t*1E6, phi_ss, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [$\mu$s]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow s(t)\ast s(-t)$');
ax.set_xlim([-2E6*t0, 2.5E6*t0]); rwth_plots.axis(ax);
rwth_plots.annotate_distance(ax, r'$2t_0$', (-t0*1E6,.1), (t0*1E6,.1)); fig.tight_layout();
Markdown("$t_0 = {:.3f}\ \mu\\mathrm{{s}}$".format(t0*1E6))
```
%% Cell type:markdown id: tags:
Tiefpassfilterung (ideal, Grenzfrequenz $f_\mathrm{g} = \frac{1}{2T}$) und ideale Abtastung
%% Cell type:code id: tags:
``` python
# Tiefpassfilterung
fg = 1/(2*T)
fTP = ient_ideal_lowpass(f, fg)
# Abtastung
nT_idx = np.arange(0, len(t), int(fs*T)); nT = t[nT_idx];
fa_nT = fTP[nT_idx];
fa = np.zeros_like(t); fa[nT_idx] = fa_nT
# Plot
rwth_media.audio_play(f, fs, r'$f(t)$')
rwth_media.audio_play(fTP, fs, r'$f(t)\ast h_\mathrm{LP}(t)$')
rwth_media.audio_play(f0, fs0, r'$f(t)$')
rwth_media.audio_play(signal.resample(fTP, len(f0)), fs0, r'$f(t)\ast h_\mathrm{LP}(t)$')
display(Markdown("$f_g = {:.1f}\ \\mathrm{{Hz}}$".format(fg)))
fig,ax=plt.subplots();
ax.plot(1000*t, fTP, 'rwth:blue', label=r'$f(t)$');
rwth_plots.plot_dirac(ax, 1000*nT, fa_nT, 'rwth:red', label=r'$f_\mathrm{a}(t)$')
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.legend();
ax.set_xlim([0, 10]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
Sendesignal
%% Cell type:code id: tags:
``` python
def sender_carrier(s, fa, t0):
m = convolution(s(t, t0), fa/deltat)
return m
m = sender_carrier(s, fa, t0)
# Plot
fig,ax=plt.subplots();
ax.plot(1000*t, m, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow m(t)$');
ax.set_xlim([0, 10]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
## Kanal
Hier wird $n(t)=0$ angenommen.
%% Cell type:code id: tags:
``` python
n = np.zeros_like(t)
```
%% Cell type:markdown id: tags:
## Empfänger
### Korrelationsfilter
%% Cell type:code id: tags:
``` python
# Schicke m(t) und n(t) getrennt durch Korrelationsfilter
def receiver_filter(m, n, s, t0):
h = s(-t, t0)
g = convolution(h, m)
ne = convolution(h, n)
y = g+ne
return (y,g,ne,h)
y,g,ne,h = receiver_filter(m, n, s, t0)
# Plot
fig,ax = plt.subplots();
ax.plot(1000*t, y, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow y(t)$');
ax.set_xlim([0, 10]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
### Abtastung und Rekonstruktion
%% Cell type:code id: tags:
``` python
def receiver_sampling_lp(y):
# Abastung
ya_nT = y[nT_idx]
ya = np.zeros_like(t)
ya[nT_idx] = ya_nT
# Rekonstruktion: Tiefpass und Multiplikation mit T
fe = ient_ideal_lowpass(ya/deltat, fg)*T
return (fe, ya, ya_nT)
fe, ya, ya_nT = receiver_sampling_lp(y)
# Plot
fig,ax=plt.subplots();
rwth_plots.plot_dirac(ax, 1000*nT, fe[nT_idx], 'rwth:black-50')
ax.plot(1000*t, fe, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow f_\mathrm{e}(t)$');
ax.set_xlim([0, 10]); rwth_plots.axis(ax);
rwth_media.audio_play(fTP, fs, r'Sender: $f(t)\ast h_\mathrm{LP}(t)$')
rwth_media.audio_play(fe, fs, r'Empfänger: $f_\mathrm{e}(t)$')
rwth_media.audio_play(signal.resample(fTP, len(f0)), fs0, r'Sender: $f(t)\ast h_\mathrm{LP}(t)$')
rwth_media.audio_play(signal.resample(fe, len(f0)), fs0, r'Empfänger: $f_\mathrm{e}(t)$')
```
%% Cell type:markdown id: tags:
### Kompensation am Empfänger
In diesem Abschnitt geht es um den in Aufgabe 11 beschriebenen Fall, in dem $t_0 > T$ gilt.
%% Cell type:code id: tags:
``` python
def receiver_compensation_filter(t0, T):
if t0 > T:
# Kompensationsfilter
H2 = 1/(1+0.4*np.cos(2*np.pi*f_ax*T))*rect(f_ax/(2*fg))
print("Kompensiere")
else:
H2 = np.ones_like(F)
print("Kompensation nicht nötig")
return H2
H2 = receiver_compensation_filter(t0, T)
# Plot
fig,ax=plt.subplots(); ax.plot(f_ax, np.abs(H2), 'rwth:blue');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow |H_2(f)|$');
ax.set_xlim([-2*fg, 2*fg]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
Filterung und Vergleich von $F_\mathrm{e}(f)$ mit $F(f)$
%% Cell type:code id: tags:
``` python
def receiver_compensate(fe, H2):
# Fourier-Transformation
Fe,_ = rwth_transforms.dft(fe, fs)
# Kompensation (im Frequenzbereich)
Fe2 = Fe*H2
# Inverse Fourier-Transformation
fe2 = np.real(rwth_transforms.idft(Fe2, len(t)))
return (fe2, Fe2, Fe)
fe2, Fe2, Fe = receiver_compensate(fe, H2)
# Plot
fig, axs = plt.subplots(2,2, figsize=(8,4));
ax = axs[0,0]; ax.plot(f_ax, np.abs(F), 'rwth:blue', label=r'$|F(f)|$');
ax.plot(f_ax, np.abs(Fe), 'rwth:green', label=r'$|F_\mathrm{e}(f)|$');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=rwth_plots.wbbox);
ax.set_xlim([-fg, fg]); ax.legend(); rwth_plots.axis(ax);
ax = axs[0,1]; ax.plot(f_ax, np.abs(F), 'rwth:blue');
ax.plot(f_ax, np.abs(Fe), 'rwth:green');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=rwth_plots.wbbox);
ax.set_xlim([0, fg/10]); rwth_plots.axis(ax);
ax = axs[1,0]; ax.plot(f_ax, np.abs(F), 'rwth:blue', label=r'$|F(f)|$');
ax.plot(f_ax, np.abs(Fe2), 'rwth:green', label=r'$|F_\mathrm{e}(f) \cdot H_2(f)|$');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=rwth_plots.wbbox);
ax.set_xlim(axs[0,0].get_xlim()); ax.set_ylim(axs[0,0].get_ylim()); ax.legend(); rwth_plots.axis(ax);
ax = axs[1,1]; ax.plot(f_ax, np.abs(F), 'rwth:blue');
ax.plot(f_ax, np.abs(Fe2), 'rwth:green');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=rwth_plots.wbbox);
ax.set_xlim(axs[0,1].get_xlim()); ax.set_ylim(axs[0,1].get_ylim()); rwth_plots.axis(ax);
rwth_media.audio_play(fTP, fs, r'$f(t)\ast h_\mathrm{LP}(t)$')
rwth_media.audio_play(fe2, fs, r'$f_\mathrm{e}(t)\ast h_2(t)$')
rwth_media.audio_play(signal.resample(fTP, len(f0)), fs0, r'$f(t)\ast h_\mathrm{LP}(t)$')
rwth_media.audio_play(signal.resample(fe2, len(f0)), fs0, r'$f_\mathrm{e}(t)\ast h_2(t)$')
```
%% Cell type:markdown id: tags:
### Interaktive Demo
%% Cell type:code id: tags:
``` python
fix,axs = plt.subplots(2,2, figsize=(8,8));
@interact(t0byT = widgets.FloatSlider(min=0.5,max=1.25,step=0.25,value=0.5,description=r'$t_0/T$:', continuous_update=False))
def update_plot(t0byT):
t0 = t0byT*T
# Sender: Träger mit Breite t0
m = sender_carrier(s, fa, t0)
# Empfänger: Korrelationsfilter, Abtastung und Tiefpass
y,g,ne,h = receiver_filter(m, n, s, t0)
fe, ya, ya_nT = receiver_sampling_lp(y)
# Empfänger: Kompensation
H2 = receiver_compensation_filter(t0, T)
fe2, Fe2, Fe = receiver_compensate(fe, H2)
# Plot
if not axs[0,0].lines:
ax = axs[0,0]; ax.plot(t/T, m, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t/T$', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow m(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim([0,50]); ax.grid(); rwth_plots.axis(ax)
ax = axs[0,1]; ax.plot(t/T, y, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t/T$', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow y(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim(axs[0,0].get_xlim()); ax.grid(); rwth_plots.axis(ax)
ax = axs[1,0]; ax.plot(t/T, fe, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t/T$', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow f_\mathrm{e}(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim(axs[0,0].get_xlim()); ax.grid(); rwth_plots.axis(ax)
ax = axs[1,1]; ax.plot(t/T, fe, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t/T$', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow f_{\mathrm{e}}(t)\ast h_2(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim(axs[0,0].get_xlim()); ax.grid(); rwth_plots.axis(ax)
else:
axs[0,0].lines[0].set_ydata(m); axs[0,1].lines[0].set_ydata(y);
axs[1,0].lines[0].set_ydata(fe); axs[1,1].lines[0].set_ydata(fe2);
```
%% Cell type:markdown id: tags:
---
This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources) (OER). Feel free to use the notebook for your own purposes. The code is licensed under the [MIT license](https://opensource.org/licenses/MIT).
Please attribute the work as follows:
*Christian Rohlfing, Übungsbeispiele zur Vorlesung "Informationsübertragung"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment