Skip to content
GitLab
Commit 30f273ce authored by Hafiz Emin Kosar's avatar Hafiz Emin Kosar
Browse files

- minor bug fixes

parent 559cb6c6
Pipeline #290198 failed with stage
in 3 minutes and 47 seconds
Hide whitespace changes
Inline Side-by-side
GDET3 Laplace-Transformation.ipynb
View file @ 30f273ce
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ipywidgets import interact, interactive, fixed, HBox, VBox
import ipywidgets as widgets
from IPython.display import clear_output, display, HTML
from IPython.display import Markdown as md
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import rwth_nb.plots.mpl_decorations as rwth_plots
from rwth_nb.misc.signals import *
import rwth_nb.misc.transforms as rwth_transforms
from src.laplace.laplace_plot import pzPoint, pzPlot
t = np.linspace(-6,6,1024)
f = np.linspace(-6,6,1024)
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Laplace-Transformation
%% Cell type:markdown id: tags:
Beispiele
$s_1(t) = \mathrm{e}^{-bt}\cdot\epsilon(t)$
transformiert sich zu
$S_1(p) = \frac{1}{b+p}$ mit Konvergenzbereich $\mathrm{Re}\{p\}>\mathrm{Re}\{-b\}$
%% Cell type:code id: tags:
``` python
b = 2
pp = np.array([-b]); pz = np.array([])
roc = np.array([-b, np.inf])
fig,axs = plt.subplots(1,2, **rwth_plots.landscape)
# Laplace-Bereich
ax = axs[0]; rwth_plots.plot_lroc(ax, roc)
ax.plot(np.real(pp), np.imag(pp), **rwth_plots.style_poles); ax.plot(np.real(pp), -np.imag(pp), **rwth_plots.style_poles)
ax.set_xlabel(r'$\rightarrow \mathrm{Re}$'); ax.set_ylabel(r'$\uparrow \mathrm{Im}$');
ax.set_xlim(-5, 5); ax.set_ylim(-5, 5); rwth_plots.grid(ax); rwth_plots.axis(ax)
# Zeitbereich
s1, _, _ = rwth_transforms.ilaplace_ht(t, 1, pp, pz, [1], [], roc)
ax = axs[1]; ax.plot(t, np.real(s1), **rwth_plots.style_graph);
ax.set_xlabel(r'$\rightarrow t$'); ax.set_ylabel(r'$\uparrow s_1(t)$'); rwth_plots.grid(ax); rwth_plots.axis(ax)
```
%% Cell type:markdown id: tags:
Da der Konvergenzbereich die imaginäre Achse umschließt, existiert auch die Fouriertransformierte $S_1(f)$:
%% Cell type:code id: tags:
``` python
# Frequenzbereich
S1f = rwth_transforms.ilaplace_Hf(f, 1, pp, pz, [1], [], dB=False)
fig, ax = plt.subplots(1,1)
ax.plot(f, S1f, **rwth_plots.style_graph);
ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow S_1(f)$');
ax.set_xlim([-5.5,5.5]); ax.set_ylim([-0.1,0.55]); rwth_plots.grid(ax); rwth_plots.axis(ax)
```
%% Cell type:markdown id: tags:
$s_2(t) = -\mathrm{e}^{-bt}\cdot\epsilon(-t)$
transformiert sich zu
$S_2(p) = \frac{1}{b+p}$ mit Konvergenzbereich $\mathrm{Re}\{p\}<\mathrm{Re}\{-b\}$:
%% Cell type:code id: tags:
``` python
roc = np.array([-np.inf, -b])
fig,axs = plt.subplots(1,2, **rwth_plots.landscape)
# Laplace-Bereich
ax = axs[0]; rwth_plots.plot_lroc(ax, roc); rwth_plots.annotate_order(ax, pp, [1]);
ax.plot(np.real(pp), np.imag(pp), **rwth_plots.style_poles); ax.plot(np.real(pp), -np.imag(pp), **rwth_plots.style_poles)
ax.set_xlabel(r'$\rightarrow \mathrm{Re}$'); ax.set_ylabel(r'$\uparrow \mathrm{Im}$');
ax.set_xlim(-5, 5); ax.set_ylim(-5, 5); rwth_plots.axis(ax); rwth_plots.grid(ax)
# Zeitbereich
s2, _, _ = rwth_transforms.ilaplace_ht(t, 1, pp, pz, [1], [], roc)
ax = axs[1]; ax.plot(t, np.real(s2), **rwth_plots.style_graph);
ax.set_xlabel(r'$\rightarrow t$'); ax.set_ylabel(r'$\uparrow s_2(t)$');
ax.set_ylim([-45,5]); ax.grid(); rwth_plots.axis(ax)
ax.set_ylim([-45,5]); rwth_plots.grid(ax); rwth_plots.axis(ax)
```
%% Cell type:markdown id: tags:
Da in diesem Fall der Konvergenzbereich links vom linkesten Pol liegt und somit nicht die imaginäre Achse beinhaltet, existiert die Fouriertransformierte nicht.
Eine ausführlichere Demo zur Laplacetransformation findet man [hier](GDET3%20Laplace-Transformation%20GUI.ipynb).
%% 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:
*Emin Kosar, Christian Rohlfing, Übungsbeispiele zur Vorlesung "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
......
GDET3 RLC-System.ipynb
View file @ 30f273ce
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ipywidgets import interact, interactive
import ipywidgets as widgets
from IPython.display import clear_output, display, HTML
import matplotlib.pyplot as plt
import numpy as np
import cmath # for sqrt(-1)
import rwth_nb.plots.mpl_decorations as rwth_plots
from rwth_nb.misc.signals import *
plt.close('all')
```
%% Cell type:markdown id: tags:
<div class="inline-block">
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# RLC-System
%% Cell type:markdown id: tags:
Es wird folgendes RLC-System betrachtet:
![Blockdiagramm](figures/rlc_system_block_diagram.png)
%% Cell type:code id: tags:
``` python
# Exemplary values
R = 16 # Ohm
L = 1.5E-3 # Henry, mH
C = 1E-6 # Farad, myF
```
%% Cell type:markdown id: tags:
## Berechnung Laplace-Übertragungsfunktion
Die Übertragungsfunktion des Systems kann im Laplace-Bereich mittels Spannungsteiler berechnet werden:
$$
H(p)
=
\frac{U_2(p)}{U_1(p)}
=
\frac{1/(Cp)}{R+Lp+1/(Cp)}
=
\frac{1}{LCp^2+RCp+1}
=
\frac{1/(LC)}{p^2+(R/L)p + 1/(LC)}
$$
mit $p = \sigma + \mathrm{j}\omega = \sigma + \mathrm{j}2 \pi f$.
%% Cell type:code id: tags:
``` python
# Laplace transfer function
H_laplace = lambda p: 1/(L*C*p**2+R*C*p+1);
```
%% Cell type:markdown id: tags:
### Polstellen
Die zugehörigen Polstellen berechnen sich zu
$$p_{\mathrm{P},1,2} =
- \underbrace{\frac{R}{2L}}_{=a} \pm \underbrace{\sqrt{\frac{R^2}{4L^2} - \frac{1}{LC}}}_{=b}
$$
mit $a=\frac{R}{2L}$ und $b=\sqrt{\frac{R^2}{4L^2}-\frac{1}{LC}}$.
%% Cell type:code id: tags:
``` python
# Poles
a = R/(2*L)
b = cmath.sqrt(R**2/(4*L**2)-1/(L*C))
p_p1 = -a+b
p_p2 = -a-b
# Print out the numbers
print("a={0:.3f}".format(a))
print("b={0:.3f}".format(b))
if R**2/(4*L**2) >= 1/(L*C): print("b reell")
else: print("b imaginär")
print("\nPolstellen:")
print("p_p1={0:.3f}, p_p2={0:.3f}".format(p_p1, p_p2))
```
%% Cell type:markdown id: tags:
### Pol-Nulstellendiagramm
Nachfolgend wird das Pol-Nulstellendiagramm geplottet. Es enthält die beiden konjugiert komplexen Polstellen, den Konvergenzbereich und das zugehörige $H_0$.
Da der Konvergenzbereich die imaginäre Achse beinhaltet, ist das System stabil.
%% Cell type:code id: tags:
``` python
beta = np.imag(b) # Imaginary part of the poles
pp = np.array([p_p1, p_p2]); pz = np.array([]) # Zeros # Poles and Zeros
ord_p = np.array([1, 1]); ord_z = np.array([]) # Poles' and Zeros' orders
roc = np.array([np.max(np.real(pp)), np.inf]) # region of convergence
# Plot
fig, ax = plt.subplots()
ax.plot(np.real(pp), np.imag(pp), **rwth_plots.style_poles); ax.plot(np.real(pp), -np.imag(pp), **rwth_plots.style_poles); rwth_plots.annotate_order(ax, pp, ord_p);
ax.plot(np.real(pz), np.imag(pz), **rwth_plots.style_zeros); ax.plot(np.real(pz), -np.imag(pz), **rwth_plots.style_zeros); rwth_plots.annotate_order(ax, pz, ord_z);
rwth_plots.plot_lroc(ax, roc, 500, np.imag(p_p1)); ax.text(-1000,ax.get_ylim()[1]*0.8,'$H_0 = 1/(LC)$',bbox = rwth_plots.wbbox);
ax.set_xlim(ax.get_xlim()[0], 300); rwth_plots.annotate_xtick(ax, r'$-a$', -a,0,'k');
ax.set_xlabel(r'$\rightarrow \mathrm{Re}\{p\}$'); ax.set_ylabel(r'$\uparrow \mathrm{Im}\{p\}$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
## Fourier-Übertragungsfunktion
Aus der Laplaceübertragungsfunktion kann die Fourierübertragungsfunktion berechnet werden, indem $p = \mathrm{j}2\pi f$ gesetzt wird:
$$H(p = \mathrm{j}2\pi f) \quad\text{mit}\quad \omega_0 = \frac{1}{\sqrt{LC}}$$
%% Cell type:code id: tags:
``` python
# Fourier transfer function
H_fourier = lambda f: H_laplace(1j*2*np.pi*f)
f = np.linspace(0, 10000, 10000) # frequency axis
# Resonance frequency
omega0 = 1/np.sqrt(L*C)
f0 = omega0/(2*np.pi)
# Print f0
print("f0 = {0:.2f} Hz".format (f0))
# Plot
fig,ax = plt.subplots()
ax.plot(f/1000, np.abs(H_fourier(f))); rwth_plots.axis(ax)
ax.set_xlabel(r'$\rightarrow f$ [kHz]'); ax.set_ylabel(r'$\uparrow |H(f)|$');
rwth_plots.annotate_xtick(ax, r'$f_0 = \omega_0/(2\pi)$', omega0/(2*np.pi)/1000,-0.25,'k');
rwth_plots.annotate_xtick(ax, r'$f_0 = \omega_0/(2\pi)$', omega0/(2*np.pi)/1000,-0.25,'rwth:black');
ax.axvline(f0/1000,0,1, color='k',linestyle='--');
```
%% Cell type:markdown id: tags:
## Impulsantwort
Die Impulsantwort kann mittels Partialbruchzerlegung (siehe Vorlesung) bestimmt werden zu
$$
h(t) = \frac{\mathrm{e}^{-at}}{bLC}
\frac{\mathrm{e}^{bt}-\mathrm{e}^{-bt}}{2}
\cdot \varepsilon(t)
$$
mit $a=\frac{R}{2L}$ und $b=\sqrt{\frac{R^2}{4L^2}-\frac{1}{LC}}$.
Der nachfolgende Plot zeigt diese Impulsantwort.
%% Cell type:code id: tags:
``` python
# Impulse response (time-domain)
h = lambda t: np.real(np.exp(-a*t)/(b*L*C)*(np.exp(b*t)-np.exp(-b*t))/2*unitstep(t))
t = np.linspace(-0.001, 0.01, 10000) # time axis
# Plot
fig,ax = plt.subplots()
ax.plot(t*1000, h(t), 'rwth:blue'); rwth_plots.axis(ax); ax.set_xlim([-0.1, 5])
ax.set_xlabel(r'$\rightarrow t$ [ms]'); ax.set_ylabel(r'$\uparrow h(t)$');
#np.mean(np.abs(h(t) - 1/(beta*L*C)*np.exp(-a*t)*np.sin(beta*t)*unitstep(t))**2)
```
%% Cell type:markdown id: tags:
## Interaktive Demo
In dieser interaktiven Demo kann das Verhalten des Systems für variable Werte von $R$ betrachtet werden. Über den Schieberegler kann der Wert für $R$ geändert werden, entsprechend sieht man die Änderungen für die Fourier-Übertragungsfunktion und die Impulsantwort.
%% Cell type:code id: tags:
``` python
fig0,axs = plt.subplots(2, 1, figsize=(8,16)); fig0.canvas.layout.height= '600px'
fig0,axs = plt.subplots(2, 1, **rwth_plots.landscape)
@widgets.interact(Rsel=widgets.FloatSlider(min=0, max=200, step=1, value=R, description='$R$ [$\Omega$]'))
def update_plot(Rsel):
H_laplace = lambda p: 1/(L*C*p**2+Rsel*C*p+1);
H_fourier = lambda f: H_laplace(1j*2*np.pi*f);
a = Rsel/(2*L)
b = cmath.sqrt(Rsel**2/(4*L**2)-1/(L*C))
h = lambda t: np.real(np.exp(-a*t)/(b*L*C)*(np.exp(b*t)-np.exp(-b*t))/2*unitstep(t))
if not axs[0].lines: # Call plot() and decorate axes. Usually, these functions take some processing time
ax = axs[0]; ax.plot(f/1000, np.abs(H_fourier(f))); rwth_plots.axis(ax)
ax = axs[0]; ax.plot(f/1000, np.abs(H_fourier(f)), 'rwth:blue'); rwth_plots.axis(ax)
ax.set_xlabel(r'$\rightarrow f$ [kHz]'); ax.set_ylabel(r'$\uparrow |H(f)|$');
ax.axvline(f0/1000,0,1, color='k',linestyle='--');
ax.axvline(f0/1000,0,1, color='rwth:black',linestyle='--');
ax = axs[1]; ax.plot(t*1000, h(t), 'rwth:blue'); rwth_plots.axis(ax); ax.set_xlim([-.225, 5])
ax.set_xlabel(r'$\rightarrow t$ [ms]'); ax.set_ylabel(r'$\uparrow h(t)$');
else: # If lines exist, replace only xdata and ydata since plt.plot() takes longer time
axs[0].lines[0].set_ydata(np.abs(H_fourier(f)))
axs[1].lines[0].set_ydata(h(t))
tmp = np.max(np.abs(H_fourier(f))); rwth_plots.update_ylim(axs[0], np.abs(H_fourier(f)), 0.1*tmp, tmp+10 )
tmp = np.max(np.abs(h(t))); rwth_plots.update_ylim(axs[1], (h(t)), 0.1*tmp, tmp+10 )
display(HTML('{}<br />{}'.format('a={0:.3f}'.format(a), 'b={0:.3f}'.format(b))))
```
%% 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 "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
......
GDET3 Reale Abtastung.ipynb
View file @ 30f273ce
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ipywidgets import interact, interactive, fixed, HBox, VBox
import ipywidgets as widgets
from IPython.display import clear_output, display, HTML
import matplotlib.pyplot as plt
import numpy as np
import rwth_nb.plots.mpl_decorations as rwth_plots
import rwth_nb.misc.transforms as rwth_transforms
from rwth_nb.misc.signals import *
signals_t = {'cos-Funktion': lambda t: np.cos(2 * np.pi * t),
'sin-Funktion': lambda t: np.sin(2 * np.pi * t),
'si-Funktion': lambda t: si(2 * np.pi * t),
'Rechteckimpuls': lambda t: rect(t / 1.05),
'Dreieckimpuls': tri}
signals_f = {'cos-Funktion': lambda f, F: np.isin(f/F, f[find_ind_least_diff(f/F, [-1, 1])]/F) * 0.5,
'sin-Funktion': lambda f, F: np.isin(f/F, f[find_ind_least_diff(f/F, [1])]/F) / 2j - np.isin(f/F, f[find_ind_least_diff(f/F, [-1])]/F) / 2j,
'si-Funktion': lambda f, F: 1/(2*np.abs(F))*rect(f / (2*F)),
'Rechteckimpuls': lambda f, F: 1/np.abs(F)*si(np.pi * f / F),
'Dreieckimpuls': lambda f, F: 1/np.abs(F)*si(np.pi * f/F) ** 2}
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Reale Abtastung
Im Gegensatz zur [idealen Abtastung](GDET3%20Ideale%20Abtastung.ipynb) werden hier zwei Verfahren zur realen Abtastung betrachtet. Tatsächlich kann nicht mit einem idealen Dirac an einem definierten Zeitpunkt abgetastet werden. Im Folgenden werden zwei Verfahren der Abtastung, die *Shape-top Abtastung* und die *Flat-top Abtastung* beschrieben.
## Shape-top Abtastung
Bei der Shape-top Abtastung wird das kontinuierliche Signal $s(t)$ mit Abstand $T=\frac{1}{r}$ abgetastet.
Anstatt einer Diracfolge wird eine Folge schmaler Rechteckimpulse mit endlicher Dauer $T_0$ verwendet. Das abgetastete Signal $s_0(t)$ ergibt sich zu
$$
s_0(t)
= s(t) \cdot \sum\limits_{n=-\infty}^\infty \mathrm{rect}\left(\frac{t-nT}{T_0}\right)
= s(t) \cdot \left[\mathrm{rect}\left(\frac{t}{T_0}\right) \ast \sum\limits_{n=-\infty}^\infty \delta(t-nT) \right].
$$
Im Frequenzbereich folgt durch Transformation
$$
S_0(f)
= S(f) \ast \left[T_0 \mathrm{si}\left(\pi T_0 f\right) \cdot \frac{1}{T}\sum_{k=-\infty}^\infty \delta(f-kr)\right]
=
S(f) \ast \left[\frac{T_0}{T} \sum_{k=-\infty}^\infty \delta\left(f-\frac{k}{T}\right) \mathrm{si} \left(\pi T_0 \frac{k}{T}\right) \right]\text{.}
$$
Auch hier entstehen wie bei der idealen Abtastung spektrale Kopien, welche um $\frac{k}{T}$ zentriert sind. Jede spektrale Kopie wird mit einem von $f$ unabhängigen Faktor $\mathrm{si} \left(\pi T_0 \frac{k}{T}\right)$ skaliert.
Der Grenzübergang $T_0\rightarrow 0$ liefert hier die ideale Abtastung: $\lim\limits_{T_0\rightarrow 0}\left(\frac{1}{T_0}s_0(t)\right) = s_\mathrm{a}(t)$.
%% Cell type:markdown id: tags:
## Flat-top Abtastung
Bei der Flat-top Abtastung wird in Intervallen endlicher Dauer abgetastet und der Signalwert dann um $T=\frac{1}{r}$ gehalten. Dieses Verfahren wird häufig in Analog-Digital-Wandlern eingesetzt. Das abgetastete Signal $s_0(t)$ ergibt sich somit zu
$$
s_0(t)
= \sum\limits_{n=-\infty}^\infty s(nT) \mathrm{rect}\left(\frac{t}{T}-\frac{1}{2}-nT\right)
= s_\mathrm{a}(t) \ast \mathrm{rect}\left(\frac{t}{T}-\frac{1}{2}\right)
= \left[s(t) \cdot \sum\limits_{n=-\infty}^\infty \delta(t-nT)\right] \ast \mathrm{rect}\left(\frac{t}{T}\right)\ast \delta\left(t-\frac{T}{2}\right) \text{.}
$$
Im Frequenzbereich folgt dann
$$
S_0(f)
= S_\mathrm{a}(f) \cdot T \cdot \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-\mathrm{j}\pi f T}
= \left[S(f) \ast \sum\limits_{k=-\infty}^\infty \delta(f-kr)\right] \cdot \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-\mathrm{j}\pi f T}\text{.}
$$
In der Demonstration kann der Unterschied zwischen der Shape-top und der Flat-top Abtastung nocheinmal anschaulich betrachtet werden.
%% Cell type:markdown id: tags:
## Beispiel
### Abtastung
Zunächst wird nocheinmal die [ideale Abtastung](GDET3%20Ideale%20Abtastung.ipynb) wiederholt. In der Abbildung ist als gestrichelte Linie das Signal $s(t)$ dargestellt, das zugehörige Spektrum in blau. Das Signal wird mit Abtastrate $r=2$ abgetastet.
Das abgetastete Signal
$s_\mathrm{a}(t) = \sum\limits_{n=-\infty}^\infty s(nT)\cdot\delta(t-nT)$ und das Spektrum des abgetasteten Signals
$S_\mathrm{a}(f) = \frac{1}{T} \sum\limits_{k=-\infty}^\infty S(f-kr)$ sind ebenfalls abgebildet.
Das Spektrum des abgetasteten Signals zeigt die für die Abtastung typischen spektralen Wiederholungen.
%% Cell type:code id: tags:
``` python
# Construction of s(t) and corresponding spectrum S(f)
t,deltat = np.linspace(-10,10,50001, retstep=True) # t-axis
f,deltaf = np.linspace(-50,50,len(t), retstep=True) # f-axis
F = 0.9 # frequency of the signal
s = lambda t: signals_t['cos-Funktion'](t*F);
S = lambda f: signals_f['cos-Funktion'](f, F);
# Ideal sampling: Construction of sa(t) and Sa(f)
r = 2; T = 1/r; # sampling rate
## Time domain
nT, snT = rwth_transforms.sample(t, s(t), T)
## Frequency domain
Sa = np.zeros_like(S(f))
kMax = 16 # number of k values in sum for Sa(f), should be infinity :)
for k in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
Sa += S(f-k/T)
Sa = Sa/T
fSadirac = f[np.where(Sa)]; Sadirac = Sa[np.where(Sa)]
# Plot
## Time domain
fig, axs = plt.subplots(2,1, **rwth_plots.landscape)
ax = axs[0]; ax.set_title('Zeitbereich');
ax.plot(t, s(t), color='rwth:blue', linestyle='--', label=r'$s(t)$');
rwth_plots.plot_dirac(ax, nT, snT, 'rwth:red', label=r'$s_\mathrm{a}(t)$')
ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); rwth_plots.grid(ax); rwth_plots.axis(ax);
## Frequency domain
ax = axs[1]; ax.set_title('Frequenzbereich');
rwth_plots.plot_dirac(ax, fSadirac, Sadirac, 'rwth:red', label=r'$S_\mathrm{a}(f)$');
rwth_plots.plot_dirac(ax, f[np.where(S(f))], S(f)[np.where(S(f))], 'rwth:blue', label=r'$S(f)$');
ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_xlabel(r'$\rightarrow f$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
txt,_=rwth_plots.annotate_xtick(ax, r'$r=2$', r, -.15, 'black'); txt.get_bbox_patch().set_alpha(1);
```
%% Cell type:markdown id: tags:
Nun wird das abgetastete Signal $s_0(t)$ im Zeitbereich betrachtet, welches mittels Flat-top Abtastung erzeugt wurde:
$$
s_0(t)
= s_\mathrm{a}(t) \ast \mathrm{rect}\left(\frac{t}{T}-\frac{1}{2}\right)
= \sum\limits_{n=-\infty}^\infty s(nT) \mathrm{rect}\left(\frac{t}{T}-\frac{1}{2}-nT\right).
$$
Die Abtastrate ist ebenfalls $r=2$.
Für die Flat-top Abtastung charakteristisch wird der Signalwert bis zum nächsten Abtastwert gehalten.
%% Cell type:code id: tags:
``` python
# Time-Domain
s0 = np.zeros_like(t) # sum of rects
Nmax = np.ceil(t[-1]/T) # Parts of the infinite rect sum, should be infinity :)
for n in np.arange(-Nmax,Nmax+1):
s0 = s0 + rect((t-n*T)/T-0.5) * s(n*T)
# Plot
fig, ax = plt.subplots(**rwth_plots.landscape); ax.set_title('Zeitbereich')
ax.plot(t, s(t), 'rwth:blue', linestyle='--', label=r'$s(t)$'); ax.plot(t, s0, 'rwth:red', label=r'$s_0(t)$')
ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-2.9, 2.9]); ax.legend(loc=2); rwth_plots.grid(ax); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
Dies hat im Frequenzbereich zur Konsequenz, dass die spektralen Kopien mit $T \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-j\pi f T}$ gewichtet werden.
$S_0(f)$ ergibt sich also zu
$$
S_0(f)
= S_\mathrm{a}(f) \cdot T \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-j\pi f T}.
$$
Die nachfolgende Abbildung verdeutlicht diesen Effekt.
%% Cell type:code id: tags:
``` python
# Frequency domain
S_si = T*si(np.pi*T*f)
S_exp = np.exp(-1j*np.pi*f*T)
S0 = Sa * S_si * S_exp
# Plot
## Magnitude
fig, axs = plt.subplots(2,1, **rwth_plots.landscape);
ax = axs[0]; ax.set_title('Betrag')
ax.plot(f, np.abs(S_si*S_exp), 'k--', label=r'$|T\mathrm{si}(\pi f T)e^{-\mathrm{j}\pi f T}|$')
fS0dirac = f[np.where(S0)]; S0dirac = S0[np.where(S0)] # sample S0 and f
rwth_plots.plot_dirac(ax, fS0dirac, np.abs(S0dirac), 'rwth:red', label=r'$|S_0(f)$|');
ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_xlabel(r'$\rightarrow f$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
## Phase
ax = axs[1]; ax.set_title('Phase')
ax.plot(f, np.angle(S_si*S_exp), 'k--', label=r'$\angle T\mathrm{si}(\pi f T)e^{-\mathrm{j}\pi f T}$')
rwth_plots.stem(ax, fS0dirac, np.angle(S0dirac), 'rwth:red', label=r'$\angle S_0(f)$')
ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_xlabel(r'$\rightarrow f$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
### Rekonstruktion
Durch die Multiplikation von $S_\mathrm{a}(f)$ mit $T \cdot \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-j\pi f T}$ wird das Spektrum $S_0(f)$ im Basisband verzerrt. So ist es nicht möglich, $S(f)$ mittels eines einfachen idealen Tiefpasses zu rekonstruieren. Zusätzlich ist ein Filter nötig, welches diesen Faktor ausgleicht, dieses ist in der folgenden Abbildung dargestellt.
$$
H_\mathrm{eq}(f) = \frac{1}{T \cdot \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-j\pi f T}}
$$
%% Cell type:code id: tags:
``` python
# Reconstruction filters #plt.close('all')
## Ideal Low pass to crop the base band
H_lp = rect(f/(r+0.001)) # ideal low pass between -r/2 and r/2
## Equalizing filter to compensate the influence of si(...) and exp(...) terms in S_0(f)
H_eq = 1/(T*si(np.pi*T*f) * np.exp(-1j*np.pi*f*T))
## Overall reconstruction filter
H = H_lp * H_eq
# Plot
fig,axs = plt.subplots(2,1, **rwth_plots.landscape)
ax = axs[0]
ax.plot(f, np.abs(H_eq), 'k--', linewidth=1, label=r'$|H_\mathrm{eq}(f)|$')
ax.plot(f, np.abs(H), 'k', label=r'$|H(f)|$')
ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_ylim([-.1,21]);
ax.set_xlabel(r'$\rightarrow f$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
ax = axs[1]
ax.plot(f, np.angle(H_eq), 'k--', linewidth=1, label=r'$\angle H_\mathrm{eq}(f)$')
ax.plot(f, np.angle(H)*H_lp, 'k', label=r'$\angle H(f)$')
ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); #ax.set_ylim([-.1,11]);
ax.set_xlabel(r'$\rightarrow f$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
Dieses Filter wird nun zur Rekonstruktion angewendet.
Das rekonstruierte Signal unter Verwendung des ausgleichenden Filters ist unten abgebildet. Das Originalsignal konnte erfolgreich rekonstruiert werden.
%% Cell type:code id: tags:
``` python
G = S0 * H*T;
g = rwth_transforms.idft(G); g = np.fft.ifftshift(np.real(g)); # IDFT
G = np.real(G); # discards small imaginary numbers close to zero
# Plot
fig, axs = plt.subplots(2, 1, **rwth_plots.landscape)
ax = axs[0]; fGdirac = f[np.where(G)]; Gdirac = G[np.where(G)]
rwth_plots.plot_dirac(ax, fGdirac, Gdirac, 'rwth:green');
ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow G(f)$', bbox=rwth_plots.wbbox);
ax.set_xlim([-7.5,7.5]); rwth_plots.grid(ax); rwth_plots.axis(ax);
ax = axs[1]; ax.plot(t, s(t), color='rwth:blue', linestyle='--', label=r'$s(t)$');
ax.plot(t/deltat/(f[-1]*2), g, color='rwth:green', linestyle='-', label=r'$g(t)$');
ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); rwth_plots.grid(ax); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
## Demo
In der folgenden interaktiven Demonstration kann nun betrachtet werden, was der Unterschied zwischen der Flat-top und der Shape-top Abtastung ist.
Dargestellt sind untereinander:
* der Zeitbereich mit dem Originalsignal und dem abgestasteten Signal
* der zugehörige Frequenzbereich mit dem Originalspektrum, dem Spektrum des abgetasteten Signal und dem Rekonstruktionsfilter
* das rekonstruierte Spektrum
* das rekonstruierte Signal im Zeitbereich.
Im Drop-Down-Menü kann die Art der Abtastung (Shape-top oder Flat-top) sowie die abzutastende Funktion (Cosinus-, Sinus-und si-Funktion, sowie Rechteck- oder Dreieckimpuls) ausgewählt werden.
Über den Schieberegler kann $F$ geändert werden, was bei Cosinus-, Sinus- und si-Funktion die Frequenz und bei Rechteck- und Dreieckimpuls die Breite beeinflusst.
%% Cell type:code id: tags:
``` python
T0 = T/4 # width of rects for shape-top sampling
plt.close(); fig, axs = plt.subplots(4, 1, **rwth_plots.landscape); plt.tight_layout();
fig, axs = plt.subplots(4, 1, **rwth_plots.landscape); plt.tight_layout();
@widgets.interact(sampling_type=widgets.Dropdown(options=['Shape-top', 'Flat-top'], description=r'Art der Abtastung:', style=rwth_plots.wdgtl_style),
s_type=widgets.Dropdown(options=list(signals_t.keys()), description=r'Wähle $s(t)$:'),
F=widgets.FloatSlider(min=0.1, max=2, value=0.9, step=.1, description=r'$F$', style=rwth_plots.wdgtl_style, continuous_update=False))
def update_plots(sampling_type, s_type, F):
s = lambda t: signals_t[s_type](t*F);
S = lambda f: signals_f[s_type](f, F);
nT, snT = rwth_transforms.sample(t, s(t), T)
# Construct sampled signal
if sampling_type == 'Shape-top':
u = np.zeros_like(t) # sum of rects
for n in np.arange(-Nmax,Nmax+1):
u = u + rect((t-n*T)/T0)
s0 = s(t) * u
elif sampling_type == 'Flat-top':
s0 = np.zeros_like(t) # sum of rects
for n in np.arange(-Nmax,Nmax+1):
s0 = s0 + rect((t-n*T)/T-0.5) * s(n*T)
# Construct sampled spectrum
if sampling_type == 'Shape-top':
S0 = np.zeros_like(S(f))
for k in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
S0 += S(f-k/T) * si(np.pi*T0*k/T)
S0 = S0*T0/T
elif sampling_type == 'Flat-top':
Sa = np.zeros_like(S(f))
for k in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
Sa += S(f-k/T)
S0 = Sa * si(np.pi*T*f) * np.exp(-1j*np.pi*f*T)
# Reconstruct g(t)
if sampling_type == 'Shape-top':
H = H_lp
G = S0 * H * T / T0;
elif sampling_type == 'Flat-top':
H = H_lp * H_eq * T
G = S0 * H
g = rwth_transforms.idft(G);
g = np.fft.ifftshift(np.real(g)); # IDFT
# Sample for plot
if s_type == 'cos-Funktion' or s_type == 'sin-Funktion':
fS0dirac = f[np.where(S0)]; S0dirac = S0[np.where(S0)]
fSdirac = f[np.where(S(f))]; Sdirac = S(f)[np.where(S(f))]
fGdirac = f[np.where(G)]; Gdirac = G[np.where(G)]
if s_type == 'sin-Funktion':
Sdirac = np.imag(Sdirac); S = lambda f: np.imag(signals_f[s_type](f, F));
G = np.imag(G); Gdirac = np.imag(Gdirac)
else:
Sdirac = np.real(Sdirac); S0 = np.real(S0); G = np.real(G); Gdirac = np.real(Gdirac)
else:
g /= (len(f)/(2*f[-1])) # Parseval :)
S0dirac = np.zeros_like(f)
G = np.real(G)
if sampling_type == 'Shape-top': # normalize to T0 for plotting reasons
S0 = S0/T0
S0dirac = S0dirac/T0
# Plot
if not axs[0].lines: # Call plot() and decorate axes. Usually, these functions take some processing time
ax = axs[0]; ax.set_title('Zeitbereich');
ax.plot(t, s(t), 'rwth:blue', linestyle='--', label=r'$s(t)$');
ax.plot(t, s0, 'rwth:red', label=r'$s_0(t)$')
ax.set_xlabel(r'$\rightarrow t$');
ax.set_xlim([-2.9, 2.9]); ax.legend(loc=2); rwth_plots.grid(ax); rwth_plots.axis(ax);
ax = axs[1]; ax.set_title('Frequenzbereich');
ax.plot(f, np.abs(H), '-', color='black', label=r'$|H(f)|$')
ax.plot(f, np.abs(H), '-', color='rwth:black', label=r'$|H(f)|$')
if s_type == 'cos-Funktion' or s_type == 'sin-Funktion':
ax.plot(f, np.ones_like(f)*np.nan, '-', color='rwth:red', label=r'$|S_0(f)|$');
rwth_plots.plot_dirac(ax, fS0dirac, np.abs(S0dirac), 'rwth:red');
else:
ax.plot(f, np.abs(S0), '-', color='rwth:red', label=r'$|S_0(f)|$');
rwth_plots.plot_dirac(ax, [], [], 'rwth:red');
ax.plot(f, S(f), color='rwth:blue', linestyle='--', linewidth=1, label=r'$S(f)$')
ax.set_xlim([-7.5,7.5]); ax.set_ylim([-1,2]); ax.legend(loc=2);
ax.set_xlabel(r'$\rightarrow f$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
txt,_=rwth_plots.annotate_xtick(ax, r'$r=2$', r, -.15, 'black'); txt.get_bbox_patch().set_alpha(1);
txt,_=rwth_plots.annotate_xtick(ax, r'$f_\mathrm{g}$', r/2, -.15, 'black'); txt.get_bbox_patch().set_alpha(1);
txt,_=rwth_plots.annotate_xtick(ax, r'$r=2$', r, -.15, 'rwth:black'); txt.get_bbox_patch().set_alpha(1);
txt,_=rwth_plots.annotate_xtick(ax, r'$f_\mathrm{g}$', r/2, -.15, 'rwth:black'); txt.get_bbox_patch().set_alpha(1);
ax = axs[2];
if s_type == 'cos-Funktion' or s_type == 'sin-Funktion':
ax.plot(f, np.ones_like(f)*np.nan, '-', color='rwth:green');
rwth_plots.plot_dirac(ax, fGdirac, Gdirac, 'rwth:green');
else:
ax.plot(f, G, '-', color='rwth:green');
rwth_plots.plot_dirac(ax, [], [], 'rwth:red');
ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow G(f)$', bbox=rwth_plots.wbbox);
ax.set_xlim([-7.5,7.5]); rwth_plots.grid(ax); rwth_plots.axis(ax);
ax = axs[3]; ax.set_title('Zeitbereich (nach Rekonstruktion)');
ax.plot(t, s(t), color='rwth:blue', linestyle='--', label=r'$s(t)$');
ax.plot(t/deltat/(f[-1]*2), g, 'rwth:green', label=r'$g(t)$');
ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); rwth_plots.grid(ax); rwth_plots.axis(ax);
plt.tight_layout()
else:
axs[0].lines[0].set_ydata(s(t)); axs[0].lines[1].set_ydata(s0); axs[1].lines[-3].set_ydata(S(f));
axs[1].lines[0].set_ydata(np.abs(H))
if s_type == 'cos-Funktion' or s_type == 'sin-Funktion': # dirac plot
rwth_plots.dirac_set_data(axs[1].containers, fS0dirac, np.abs(S0dirac));
axs[1].lines[1].set_ydata(np.ones_like(f)*np.nan);
rwth_plots.dirac_set_data(axs[2].containers, fGdirac, Gdirac); axs[2].lines[0].set_ydata(np.ones_like(f)*np.nan);
else:
axs[1].lines[1].set_ydata(np.abs(S0)); rwth_plots.dirac_set_data(axs[1].containers, [], []);
axs[2].lines[0].set_ydata(G); rwth_plots.dirac_set_data(axs[2].containers, [], []);
axs[3].lines[0].set_ydata(s(t)); axs[3].lines[1].set_ydata(g);
if s_type == 'sin-Funktion': # Adapt labels
axs[1].lines[-3].set_label(r'$\mathrm{Im}\{S(f)\}$');
axs[2].yaxis.label.set_text(r'$\uparrow \mathrm{Im}\{G(f)\}$');
else:
axs[1].lines[-3].set_label(r'$S(f)$')
axs[2].yaxis.label.set_text(r'$\uparrow G(f)$');
axs[1].legend(loc=2)
tmp = np.concatenate((np.abs(S0),np.abs(S0dirac),np.abs(H))); rwth_plots.update_ylim(axs[1], tmp, 0.19, np.max(np.abs(tmp))); rwth_plots.update_ylim(axs[2], G, 0.19, np.max(G));
rwth_plots.update_ylim(axs[3], np.concatenate((s(t),g)), 0.19, np.max(np.concatenate((s(t),g))));
```
%% Cell type:markdown id: tags:
### Aufgaben
Wähle eine beliebige Funktion aus.
* Wie sieht das abgetastete Signal mit Shape-top Abtastung aus und wie mit Flat-top Abtastung?
* Welche Auswirkungen hat die Abtastungsart auf das Spektrum und die Rekonstruktion?
* Variiere $F$. Was passiert?
* Führe diese Beobachtung für unterschiedliche Funktionen aus.
%% 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 "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
......
GDET3 Zufallssignale in LTI-Systemen.ipynb
View file @ 30f273ce
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ipywidgets import interact, interactive, fixed
import ipywidgets as widgets
from IPython.display import clear_output, display, HTML, Audio
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
import rwth_nb.plots.mpl_decorations as rwth_plots
import rwth_nb.misc.media as rwth_media
import rwth_nb.misc.filters as rwth_filters
from rwth_nb.misc.signals import *
fs = 22050 # Samplingrate
(t,deltat) = np.linspace(0,10, 10*fs, retstep=True) # Zeitachse
(x,deltax) = np.linspace(-5,5, 1024, retstep=True) # VDF
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Zufallssignale in LTI-Systemen
%% Cell type:markdown id: tags:
Für jede Musterfunktion ${}^k s(t)$ eines stochastischen Prozesses gilt bei Übertragung über ein LTI-System der Impulsantwort $h(t)$ das Faltungsprodukt ${}^k s(t)∗h(t)={}^{k}g(t)$. Wird ein Zufallssignal ${}^k n(t)$ über ein LTI-System übertragen, kann dies durch folgendes Blockschaltbild dargestellt werden:
![Blockdiagramm](figures/white_noise_lti_block_diagram.png)
Der Ausgangsprozess $g(t)$ kann wie jeder Zufallsprozess durch Mittelwerte und Verbundmittelwerte beschrieben werden.
Der hier verwendete Zufallsprozess $n(t)$ ist gleichverteiltes, weißes Rauschen mit
$$p_n(x) = \frac{1}{a} \mathrm{rect}\left(\frac{x-m_n}{a}\right)$$
mit Mittelwert $m_n$ und Varianz $\sigma_n^2=\frac{a^2}{12}$.
Beides ist in der nachfolgenden Abbildung dargestellt. Ebenso kann man das Rauschsignal anhören.
%% Cell type:code id: tags:
``` python
# Sample ^k n(t)
m_n = 0 # Mittelwert
sigma_n = 1 # Varianz
a = np.sqrt(12)*sigma_n
n = np.random.uniform(-0.5*a+m_n, 0.5*a+m_n, len(t))
# Gemessene Verteilungsdichtefunktion
pn_meas,bins = np.histogram(n,bins=200,range=(-5,5),density=True)
x_meas = (bins[:-1] + bins[1:]) / 2; x_meas0 = x_meas
# Berechnete Verteilungsdichtefunktion
pn = 1/a*rect((x-m_n)/a)
# Plot
fig,axs = plt.subplots(2,1,figsize=(8,4));
ax = axs[0]; ax.plot(1000*t, n, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow {}^k n(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim([0,11]); rwth_plots.axis(ax);
ax = axs[1]
ax.plot(x_meas, pn_meas, 'rwth:blue'); ax.plot(x,pn,'k--');
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_n(x)$');
rwth_plots.update_ylim(ax,pn_meas,0.1); rwth_plots.axis(ax);
rwth_media.audio_play(n,fs)
```
%% Cell type:markdown id: tags:
## LTI System
### Tiefpassfilterung
Zunächst wird der gleichverteilte Rauschprozess über einen Tiefpass $h(t)$ mit Grenzfrequenz $f_\mathrm{g}$ übertragen. Der Betrag der Übertragungsfunktion ist in der nachfolgenden Grafik dargestellt.
%% Cell type:code id: tags:
``` python
# Filter requirements.
order = 4
fg = 500 # desired cutoff frequency of the filter, Hz
b, a = rwth_filters.butter_lowpass(fg, fs, order) # generate filter coefficients
w, H = signal.freqz(b, a, worN=len(t)) # compute H(f)=H(z=e^(j 2 pi f)) out of b, a
f = 0.5*fs*w/np.pi
fig,ax = plt.subplots(1,1); ax.plot(f, np.abs(H), 'rwth:blue');
ax.axvline(fg, color='k',linestyle='--',lw=0.5); # cutoff frequency
ax.axvline(fg, color='rwth:black',linestyle='--',lw=0.5); # cutoff frequency
ax.set_xlabel(r'$\rightarrow f$ [Hz]'); ax.set_ylabel(r'$\uparrow |H(f)|$', bbox=rwth_plots.wbbox);
ax.set_xlim([0,4000]); ax.set_ylim([-.25,1.19]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
${}^k n(t)$, die $k$-te Realisierung von $n(t)$, wird nun über $h(t)$ übertragen:
![Blockdiagramm](figures/white_noise_lti_block_diagram.png)
Das Eingangssignal ${}^k n(t)$ ist in der folgenden Abbildung dargestellt. Darunter befindet sich der Plot des tiefpassgefilterten Ausgangssignals. Beide Signale stehen auch als Hörbeispiel zur Verfügung. Der Effekt der Tiefpassfilterung ist sowohl optisch als auch akustisch deutlich zu erkennen.
%% Cell type:code id: tags:
``` python
g = rwth_filters.filter(n, b, a)
fig,axs = plt.subplots(2,1);
ax = axs[0]; ax.plot(1000*t, n, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow {}^k n(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim([0,22]); rwth_plots.axis(ax);
ax = axs[1]; ax.plot(1000*t, g, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=rwth_plots.wbbox); ax.set_ylabel(r'$\uparrow {}^k g(t)$', bbox=rwth_plots.wbbox);
ax.set_xlim([0,22]); rwth_plots.axis(ax);
rwth_media.audio_play(n, fs, r'${}^k n(t)$')
rwth_media.audio_play(g, fs, r'${}^k g(t)$')
```
%% Cell type:markdown id: tags:
Nun werden die Verteilungsdichtefunktionen des Eingangs- und Ausgangssignals gemessen.
Für die Verteilungsdichtefunktion $p_n(x)$ des Rauschsignals wird näherungsweise eine Gleichverteilung festgestellt.
Die Verteilungsdichtefunktion des Ausgangssignals $p_g(x)$ ist deutlich schmaler und nicht mehr gleichverteilt, sondern eher gauß-verteilt.
%% Cell type:code id: tags:
``` python
# Gemessene Verteilungsdichtefunktionen
pg_meas,bins = np.histogram(g, bins=200, range=(-5,5), density=True)
x_meas = (bins[:-1] + bins[1:]) / 2
fig,axs = plt.subplots(2,1);
ax = axs[0]; ax.plot(x_meas0, pn_meas, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_n(x)$');
ax.set_ylim([0,0.39]); rwth_plots.axis(ax);
ax = axs[1]; ax.plot(x_meas, pg_meas, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_g(x)$'); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
### Hochpassfilterung
Nun sei $h(t)$ ein Hochpass mit Grenzfrequenz $f_\mathrm{g}$, wie in der folgenden Abbildung dargestellt. Ein- und Ausgangssignal können auch hier angehört werden.
%% Cell type:code id: tags:
``` python
fg = 1000 # Hz
b, a = rwth_filters.butter_highpass(fg, fs, order) # generate filter coefficients
w, H = signal.freqz(b, a, worN=len(t)) # compute H(f)=H(z=e^(j 2 pi f)) out of b, a
f = 0.5*fs*w/np.pi
fig,ax = plt.subplots(1,1); ax.plot(f, np.abs(H), 'rwth:blue');
ax.axvline(fg, color='rwth:black',linestyle='--',lw=0.5); # cutoff frequency
ax.set_xlabel(r'$\rightarrow f$ [Hz]'); ax.set_ylabel(r'$\uparrow |H(f)|$');
ax.set_xlim([0,4000]); ax.set_ylim([-.25,1.19]); rwth_plots.axis(ax);
g = rwth_filters.filter(n, b, a)
rwth_media.audio_play(n, fs, r'${}^k n(t)$')
rwth_media.audio_play(g, fs, r'${}^k g(t)$')
```
%% Cell type:markdown id: tags:
Nachfolgend wird noch die gemessene Verteilungsdichtefunktion des Rauschprozesses nach der Hochpassfilterung geplottet. Es ist zu erkennen, dass hier zwar eine eher gauß-verteilte Verteilungsdichte erzeugt wird, diese aber von der Breite und Höhe ähnlich zu der Gleichverteilung des Eingangsprozesses ist. Im Hörbeispiel ist auch zu hören, dass diese Hochpassfilterung auf das Rauschen einen sehr viel geringeren akustischen Einfluss hat, als die vorherige Tiefpassfilterung.
%% Cell type:code id: tags:
``` python
# Gemessene Verteilungsdichtefunktion
pg_meas,bins = np.histogram(g, bins=200, range=(-5,5), density=True)
x_meas = (bins[:-1] + bins[1:]) / 2
fig,axs = plt.subplots(2,1);
ax = axs[0]; ax.plot(x_meas0, pn_meas, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_n(x)$');
ax.set_ylim([0,0.39]); rwth_plots.axis(ax);
ax = axs[1]; ax.plot(x_meas, pg_meas, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_g(x)$'); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
## Interaktive GUI
In dieser interaktiven GUI kann ausprobiert werden, welchen Effekt eine Hoch- oder Tiefpassfilterung auf einen gleichverteilten Rauschprozess haben. Über das Drop-Down-Menü kann der gewünschte Filtertyp ausgewählt werden. Die Grenzfrequenz des Filters ist über den Schieberegler anpassbar.
Geplottet wird die erzeugte Verteilungsdichtefunktion des Ausgangssignal $p_g(x)$, ebenso kann der Effekt der Filterung angehört werden.
%% Cell type:code id: tags:
``` python
fig,axs = plt.subplots(2,1);
@widgets.interact(fg=widgets.FloatSlider(min=300, max=10000, step=100, description='$f_g$ [Hz]', continuous_update=False),
btype= widgets.Dropdown(options=['Tiefpass','Hochpass'], description='Filtertyp'))
def update_plot(fg, btype):
# LTI
b,a = rwth_filters.butter(fg, fs, order, btype)
g = rwth_filters.filter(n, b,a)
rwth_media.audio_play(g, fs)
# Histogram
pg_meas,bins = np.histogram(g, bins=200, range=(-5,5), density=True)
x_meas = (bins[:-1] + bins[1:]) / 2
# Update plot
if not axs[0].lines:
ax = axs[0]; ax.plot(x_meas0, pn_meas, 'rwth:blue');
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_n(x)$');
rwth_plots.update_ylim(ax,pn_meas,0.1); rwth_plots.axis(ax);
ax = axs[1]; ax.plot(x_meas, pg_meas, 'rwth:blue');
rwth_plots.update_ylim(ax,pg_meas,0.1);
ax.set_xlabel(r'$\rightarrow x$'); ax.set_ylabel(r'$\uparrow p_g(x)$'); rwth_plots.axis(ax);
else:
axs[1].lines[0].set_ydata(pg_meas)
```
%% Cell type:markdown id: tags:
## Aufgaben
* Höre dir zur Erinnerung nocheinmal das Eingangssignal an.
* Wähle eine kleine Grenzfrequenz aus. Wie hört sich das tiefpassgefilterte Signal im Vergleich zum Eingangssignal an, wie das hochpassgefilterte?
* Erhöhe die Grenzfrequenz. Welchen Effekt hat dies auf das tiefpass- bzw hochpassgefilterte Signal?
* Betrachte ebenfalls die Verteilungsdichtefunktion $p_g(x)$ des gefilterten Signals. Kannst du das Verhalten der Verteilungsdichten erklären?
%% 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 "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
......
GDET3 Übertragungsfunktion.ipynb
View file @ 30f273ce
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
import ipywidgets as widgets
from ipywidgets import interact, interactive, fixed, Layout
from IPython.display import clear_output, display, HTML
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal # convolution
import rwth_nb.plots.mpl_decorations as rwth_plots
from rwth_nb.misc.signals import *
import rwth_nb.misc.transforms as rwth_transforms
def convolution(s, h):
# Convolve s and h numerically
return signal.convolve(s(t), h(t), mode='same')*deltat
(t,deltat) = np.linspace(-20, 20, 50001, retstep=True) # t Achse
(f,deltaf) = np.linspace(-20, 20, 50001, retstep=True) # f Achse
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Übertragungsfunktion
Zum Starten: Im Menü: Run <span class="fa-chevron-right fa"></span> Run All Cells auswählen.
In diesem Beispiel wird die Übertragung eines Eingangssignal $s(t)$ über ein System mit der Impulsantwort $h(t)$ und der zugehörigen Übertragungsfunktion $H(f)$ gezeigt.
## Eingangssignal
Das verwendete Eingangssignal ist ein Rechteck mit der Breite $T_0=4$:
$\displaystyle s(t) = \frac{1}{T_0}\mathrm{rect}\left(\frac{t}{T_0}\right)$. Dieses ist in der folgenden Abbildung dargestellt.
%% Cell type:code id: tags:
``` python
s = lambda t: 1/T0*rect(t/T0)
# PLot
T0 = 4
plt.close(); fig,ax=plt.subplots();
fig,ax=plt.subplots();
ax.plot(t, s(t), 'rwth:blue');
ax.set_xlabel(r'$\rightarrow t$'); ax.set_ylabel(r'$\uparrow s(t)$')
ax.set_xlim([-5,5]); ax.set_ylim([0,0.5]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
Das System hat die Impulsantwort $\displaystyle h(t) = \frac{1}{T}\varepsilon(t)\mathrm{e}^{-t/T}$ und die Übertragungsfunktion $\displaystyle H(f) = \frac{1}{1+\mathrm{j}2 \pi f T}$ mit $T=RC$. Der Betrag der Übertragungsfunktion ist nachfolgend geplottet.
%% Cell type:code id: tags:
``` python
T=2 # RC system with T=R*C
h = lambda t: 1/T*unitstep(t)*np.exp(-t/T) # impulse response
H = lambda f: 1/(1+2j*np.pi*f*T) # transfer function
# Plot
fig,ax=plt.subplots();
ax.plot(f, np.abs(H(f)), 'rwth:blue');
ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow |H(f)|$')
ax.set_xlim([-7.5,7.5]); rwth_plots.axis(ax);
#Hf, f0 = rwth_transforms.dft(h(t), 1/(deltat)); Hf = Hf *2*f[-1]; ax.plot(f0, np.abs(Hf), 'k--');
```
%% Cell type:markdown id: tags:
Das Ausgangssignal $g(t)$ kann nun als Faltung des Eingangssignals $s(t)$ mit der Impulsantwort $h(t)$ des Systems beschrieben werden.
$$g(t) = s(t) \ast h(t)$$
Die Impulsantwort sowie das aus der Faltung resultierende Ausgangssignal sind in der folgenden Abbildung dargestellt.
%% Cell type:code id: tags:
``` python
gt = convolution(s, h) # numerical convolution
# Plot
fig,ax = plt.subplots();
ax.plot(t, h(t), 'k--', label=r'$h(t)$');
ax.plot(t,gt, 'rwth:blue', label=r'$g(t)$');
ax.set_xlabel(r'$\rightarrow t$'); ax.legend(loc=2)
ax.set_xlim([-11.5,11.5]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
Anstatt im Zeitbereich eine Faltung durchzuführen, kann das Ergebnis auch über eine Multiplikation im Frequenzbereich berechnet werden. Hierfür muss die Fouriertransformierte $S(f)$ des Eingangssignals mit der Übertragungsfunktion $H(f)$ des Systems multipliziert werden:
$G(f) = S(f) H(f)$
Aus $G(f)$ kann dann mittels Rücktransformation $G(t)$ berechnet werden.
%% Cell type:code id: tags:
``` python
Gf, f0 = rwth_transforms.dft(gt, 1/(deltat)); Gf = Gf *2*f[-1];
# Plot
fig,ax = plt.subplots();
ax.plot(f, np.abs(H(f)), 'k--', label=r'$|H(f)|$');
ax.plot(f0, np.abs(Gf), 'rwth:blue', label=r'$|G(f)|$');
ax.set_xlabel(r'$\rightarrow f$'); ax.legend(loc=2)
ax.set_xlim([-7.5,7.5]); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
## Demonstration
In dieser Demo kann das Verhalten der Übertragung eines Rechteckimpuls über ein System mit der Impulsantwort $h(t)$ bei Variation der Breite $T_0$ des Rechteckimpulses betrachtet werden. Auf der linken Seite sind die Eingangsfunktion $s(t)$, sowie der Betrag der Übertragungsfunktion $|S(f)|$ dargestellt. Auf der rechten Seite ist oben die Impulsantwort $h(t)$ des Systems sowie das resultierende Ausgangssignal $g(g)$ im Zeitbereich zu sehen. Darunter sind die zugehörigen Beträge der Funktionen im Frequenzbereich dargestellt.
%% Cell type:code id: tags:
``` python
plt.close(); fig,axs = plt.subplots(2,2, gridspec_kw = {'width_ratios':[1, 2]});
@widgets.interact(T0=widgets.FloatSlider(min=0.1, max=2, step=0.1, value=1, description='$T_0$', continuous_update=True))
def update_plot(T0):
s = lambda t: 1/T0*rect(t/T0)
h = lambda t: 1/T*unitstep(t)*np.exp(-t/T) # RC system with T=R*C
H = lambda f: 1/(1+2j*np.pi*f*T)
gt = convolution(s, h) # numerical convolution
Sf, f0 = rwth_transforms.dft(s(t), 1/(deltat)); Sf = Sf *2*f[-1];
Gf, _ = rwth_transforms.dft(gt, 1/(deltat)); Gf = Gf *2*f[-1];
# Update plot
if not axs[0,0].lines:
ax = axs[0,0];
ax.plot(t, s(t), 'rwth:blue')
ax.set_xlabel(r'$\rightarrow t$'); ax.set_ylabel(r'$\uparrow s(t)$');
ax.set_xlim([-1.1,1.1]); rwth_plots.axis(ax);
ax = axs[1,0];
ax.plot(f0, np.abs(Sf), 'rwth:blue')
ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow |S(f)|$');
ax.set_xlim([-11.5,11.5]); rwth_plots.axis(ax);
ax = axs[0,1]; ax.set_title('Zeitbereich')
ax.plot(t, h(t), 'k--', label=r'$h(t)$');
ax.plot(t, gt, 'rwth:blue', label=r'$g(t)$');
ax.set_xlabel(r'$\rightarrow t$'); ax.legend(loc=2)
ax.set_xlim([-11.5,11.5]); rwth_plots.axis(ax);
ax = axs[1,1]; ax.set_title('Frequenzbereich')
ax.plot(f, np.abs(H(f)), 'k--', label=r'$|H(f)|$');
ax.plot(f0, np.abs(Gf), 'rwth:blue', label=r'$|G(f)|$');
ax.set_xlabel(r'$\rightarrow f$'); ax.legend(loc=2)
ax.set_xlim([-11.5,11.5]); rwth_plots.axis(ax);
plt.tight_layout();
else:
axs[0,0].lines[0].set_ydata(s(t))
axs[1,0].lines[0].set_ydata(np.abs(Sf))
axs[0,1].lines[1].set_ydata(gt)
axs[1,1].lines[1].set_ydata(np.abs(Gf))
rwth_plots.update_ylim(axs[0,0], s(t), 0.19, np.max(s(t)));
```
%% Cell type:markdown id: tags:
## Aufgaben
* Variiere $T_0$. Wie sieht der Rechteckimpuls bei kleinen Werten, wie bei großen Werten aus?
* Was passiert bei Änderung von $T_0$ mit dem Betrag der Übertragungsfunktion $|S(f)|$?
* Wie ändert sich das Ausgangssignal $g(t)$ mit kleiner werdendem $T_0$? Was passiert im Frequenzbereich?
* Das kleinstmögliche $T_0$ ist hier 0.1. Wieso überlagern sich die Kurven für $g(t)$ und $h(t)$ bzw in diesem Fall fast?
* Was würde passieren, wenn $T_0$ noch kleiner wird? Wie sähe die zugehörige Übertragungsfunktion des Eingangssignals $|S(f)|$ aus?
%% 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 "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment