- minor bugfix in real sampling

GDET3 Reale Abtastung.ipynb

@@ -424,7 +424,7 @@
     "            axs[2].yaxis.label.set_text(r'$\\uparrow G(f)$');\n",
     "        axs[1].legend(loc=2)\n",
     "        \n",
-    "    tmp = np.concatenate((np.abs(S0),np.abs(S0dirac))); ient_update_ylim(axs[1], tmp, 0.19, np.max(np.abs(tmp))); ient_update_ylim(axs[2], G, 0.19, np.max(G));\n",
+    "    tmp = np.concatenate((np.abs(S0),np.abs(S0dirac),np.abs(H))); ient_update_ylim(axs[1], tmp, 0.19, np.max(np.abs(tmp))); ient_update_ylim(axs[2], G, 0.19, np.max(G));\n",
     "    ient_update_ylim(axs[3], np.concatenate((s(t),g)), 0.19, np.max(np.concatenate((s(t),g))));"
    ]
   },
 %% Cell type:code id: tags:
 
 ``` python
 # Copyright 2019 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 ient_nb.ient_plots import *
 from ient_nb.ient_transforms import *
 from ient_nb.ient_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[findIndOfLeastDiff(f/F, [-1, 1])]/F) * 0.5,
              'sin-Funktion':   lambda f, F: np.isin(f/F, f[findIndOfLeastDiff(f/F, [1])]/F) / 2j - np.isin(f/F, f[findIndOfLeastDiff(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="ient_nb/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.
 
 ## Shape-top Abtastung
 
 Abtastung in Intervallen endlicher Dauer $T_0$ mit Abstand $T=\frac{1}{r}$
 
 $$
 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
 $$
 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]
 $$
 
 Jede spektrale Kopie, zentriert bei $\frac{k}{T}$, wird skaliert mit von $f$ unabhängigen Faktor $\mathrm{si} \left(\pi T_0 \frac{k}{T}\right)$.
 
 Grenzübergang liefert 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
 
 Abtastung in Intervallen endlicher Dauer und um $T=\frac{1}{r}$ gehaltene Signalwerte
 
 $$
 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)
 $$
 
 Im Frequenzbereich
 $$
 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}
 $$
 
 Dieses Verfahren wird häufig in Analog-Digital-Wandlern eingesetzt.
 
 %% Cell type:markdown id: tags:
 
 ### Abtastung
 
 Zunächst wird [ideale Abtastung](GDET3%20Ideale%20Abtastung.ipynb) wiederholt mit
 $s_\mathrm{a}(t) = \sum\limits_{n=-\infty}^\infty s(nT)\cdot\delta(t-nT)$ und Spektrum
 $S_\mathrm{a}(f) = \frac{1}{T} \sum\limits_{k=-\infty}^\infty S(f-kr)$.
 
 %% 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 = ient_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)
 ax = axs[0]; ax.set_title('Zeitbereich');
 ax.plot(t, s(t), color='rwth', linestyle='--', label=r'$s(t)$');
 ient_plot_dirac(ax, nT, snT, 'rot', label=r'$s_\mathrm{a}(t)$')
 ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ient_grid(ax); ient_axis(ax);
 ## Frequency domain
 ax = axs[1];  ax.set_title('Frequenzbereich');
 ient_plot_dirac(ax, fSadirac, Sadirac, 'rot', label=r'$S_\mathrm{a}(f)$');
 ient_plot_dirac(ax, f[np.where(S(f))], S(f)[np.where(S(f))], 'rwth', label=r'$S(f)$');
 ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_xlabel(r'$\rightarrow f$'); ient_grid(ax); ient_axis(ax);
 txt,_=ient_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)
 $$
 
 %% 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(); ax.set_title('Zeitbereich')
 ax.plot(t, s(t), 'rwth', linestyle='--', label=r'$s(t)$'); ax.plot(t, s0, 'rot', label=r'$s_0(t)$')
 ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-2.9, 2.9]); ax.legend(loc=2); ient_grid(ax); ient_axis(ax);
 ```
 
 %% Cell type:markdown id: tags:
 
 Nun $S_0(f)$ im Frequenzbereich
 $$
 S_0(f)
 = S_\mathrm{a}(f) \cdot T \mathrm{si}(\pi f T) \cdot \mathrm{e}^{-j\pi f T}
 $$
 
 %% 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);
 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
 ient_plot_dirac(ax, fS0dirac, np.abs(S0dirac), 'rot', label=r'$|S_0(f)$|');
 ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_xlabel(r'$\rightarrow f$'); ient_grid(ax); ient_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}$')
 ient_stem(ax, fS0dirac, np.angle(S0dirac), 'rot', label=r'$\angle S_0(f)$')
 ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ax.set_xlabel(r'$\rightarrow f$'); ient_grid(ax); ient_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 zum Ausgleich nötig
 $$
 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)
 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$'); ient_grid(ax); ient_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$'); ient_grid(ax); ient_axis(ax);
 ```
 
 %% Cell type:markdown id: tags:
 
 Dieses Filter wird nun zur Rekonstruktion angewendet:
 
 %% Cell type:code id: tags:
 
 ``` python
 G = S0 * H*T;
 g = ient_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)
 ax = axs[0]; fGdirac  = f[np.where(G)]; Gdirac = G[np.where(G)]
 ient_plot_dirac(ax, fGdirac, Gdirac, 'grun');
 ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow G(f)$', bbox=ient_wbbox);
 ax.set_xlim([-7.5,7.5]); ient_grid(ax); ient_axis(ax);
 
 ax = axs[1]; ax.plot(t, s(t), color='rwth', linestyle='--', label=r'$s(t)$');
 ax.plot(t/deltat/(f[-1]*2), g, color='grun', linestyle='-', label=r'$g(t)$');
 ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ient_grid(ax); ient_axis(ax);
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Demo
 
 %% Cell type:code id: tags:
 
 ``` python
 T0 = T/4 # width of rects for shape-top sampling
 
 plt.close(); fig, axs = plt.subplots(4, 1); fig.canvas.layout.height = '800px'; plt.tight_layout();
 @widgets.interact(sampling_type=widgets.Dropdown(options=['Shape-top', 'Flat-top'], description=r'Art der Abtastung:', style=ient_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=ient_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 = ient_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 = ient_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', linestyle='--', label=r'$s(t)$');
         ax.plot(t, s0, 'rot', label=r'$s_0(t)$')
         ax.set_xlabel(r'$\rightarrow t$');
         ax.set_xlim([-2.9, 2.9]); ax.legend(loc=2); ient_grid(ax); ient_axis(ax);
 
         ax = axs[1];  ax.set_title('Frequenzbereich');
         ax.plot(f, np.abs(H), '-', color='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='rot', label=r'$|S_0(f)|$');
             ient_plot_dirac(ax, fS0dirac, np.abs(S0dirac), 'rot');
         else:
             ax.plot(f, np.abs(S0), '-', color='rot', label=r'$|S_0(f)|$');
             ient_plot_dirac(ax, [], [], 'rot');
         ax.plot(f, S(f), color='rwth', 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$'); ient_grid(ax); ient_axis(ax);
         txt,_=ient_annotate_xtick(ax, r'$r=2$', r, -.15, 'black'); txt.get_bbox_patch().set_alpha(1);
         txt,_=ient_annotate_xtick(ax, r'$f_\mathrm{g}$', r/2, -.15, '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='grun');
             ient_plot_dirac(ax, fGdirac, Gdirac, 'grun');
         else:
             ax.plot(f, G, '-', color='grun');
             ient_plot_dirac(ax, [], [], 'rot');
         ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow G(f)$', bbox=ient_wbbox);
         ax.set_xlim([-7.5,7.5]); ient_grid(ax); ient_axis(ax);
 
         ax = axs[3]; ax.set_title('Zeitbereich (nach Rekonstruktion)');
         ax.plot(t, s(t), color='rwth', linestyle='--', label=r'$s(t)$');
         ax.plot(t/deltat/(f[-1]*2), g, 'grun', label=r'$g(t)$');
         ax.set_xlabel(r'$\rightarrow t$'); ax.set_xlim([-7.5,7.5]); ax.legend(loc=2); ient_grid(ax); ient_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
             ient_dirac_set_data(axs[1].containers, fS0dirac, np.abs(S0dirac));
             axs[1].lines[1].set_ydata(np.ones_like(f)*np.nan);
             ient_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)); ient_dirac_set_data(axs[1].containers, [], []);
             axs[2].lines[0].set_ydata(G); ient_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))); ient_update_ylim(axs[1], tmp, 0.19, np.max(np.abs(tmp))); ient_update_ylim(axs[2], G, 0.19, np.max(G));
+    tmp = np.concatenate((np.abs(S0),np.abs(S0dirac),np.abs(H))); ient_update_ylim(axs[1], tmp, 0.19, np.max(np.abs(tmp))); ient_update_ylim(axs[2], G, 0.19, np.max(G));
     ient_update_ylim(axs[3], np.concatenate((s(t),g)), 0.19, np.max(np.concatenate((s(t),g))));
 ```
 
 %% 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, 2019, Institut für Nachrichtentechnik, RWTH Aachen University.