Verified Commit 2b31941e authored by Christian Rohlfing's avatar Christian Rohlfing
Browse files

minor changes in discrete FFT demo

parent bd7998f8
%% 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
from rwth_nb.plots import colors
import rwth_nb.plots.mpl_decorations as rwth_plt
from rwth_nb.misc.signals import *
from rwth_nb.misc.transforms import *
t,deltat = np.linspace(-10,10,50001, retstep=True) # t-axis
f,deltaf = np.linspace(-10,10,len(t), retstep=True) # f-axis
kMax = 16 # number of k values in sum for Sa(f)
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Demonstrator Diskrete Fourier-Transformation
Zum Starten: Im Menü: Run <span class="fa-chevron-right fa"></span> Run All Cells auswählen.
## Einleitung
$s(t)$ wird ideal mit Abtastrate $r$ abgetastet.
Diskrete Fourier-Transformation mit $F = \frac{1}{M}$
Diskrete Fourier-Transformation mit $F = \frac{1}{M}$ resultiert in diskreter Fourier-Transformierten
$$
S_\mathrm{d}(k)
=
\sum_{n=0}^{M-1} s_\mathrm{d}(n) \mathrm{e}^{-\mathrm{j}2\pi kFn} \quad k=0,\ldots,M-1 \text{ ,}
$$
wobei $s_\mathrm{d}(n)$ um $M$ periodisches, diskretes Zeitsignal mit $t = n \cdot \frac{1}{r}$.
Für die Rücktransformation in den Zeitbereich gilt dann
$$
s_\mathrm{d}(n)
=
\frac{1}{M} \sum_{k=0}^{M-1}S_\mathrm{d}(k) \mathrm{e}^{\mathrm{j}2\pi kFn} \quad n=0,\ldots,M-1 \text{ ,}
$$
wobei
$f = k \cdot \frac{r}{M}$.
## Beispiel
Zeige eine Periode mit $M=24$ Werten
%% Cell type:code id: tags:
``` python
M = 24; # period of discrete signal
r = 4; # sampling rate
n = np.arange(M) # discrete time axis
sd = tri(n/r) + tri((n-M)/r) # discrete signal sd(n)
F = 1/M
Sd = np.zeros_like(sd, dtype=np.complex)
for k in np.arange(M): # very slow computation for teaching purpose! Fasten your seatbelts with FFT :)
# Very slow DFT computation for teaching purposes!
# Fasten your seatbelts with FFF by uncommenting the line with "fft.fft" below
Sd = np.zeros_like(sd, dtype=complex)
for k in np.arange(M):
Sd[k] = np.sum(sd * np.exp(-2j*np.pi*k*F*n))
#Sd = np.fft.fft(sd)
Sd = np.real(Sd) # discard super small imaginary part (1e-16)
#Sd = np.fft.fft(sd) # this would be so much faster!
Sd = np.real(Sd) # discard super small imaginary part (1e-16), works "only" for examples like this
# Plot
fig,axs = plt.subplots(2,1)
## Plot discrete time domain
ax = axs[0]; rwth_plt.stem(ax, n, sd);
ax.set_xlabel(r'$\rightarrow n$'); ax.set_ylabel(r'$\uparrow s_\mathrm{d}(n)$')
rwth_plt.axis(ax);
# Plot discrete Fourier domain
ax = axs[1]; rwth_plt.stem(ax, np.arange(M), Sd/r);
ax.set_xlabel(r'$\rightarrow k$'); ax.set_ylabel(r'$\uparrow S_\mathrm{d}(k)/r$')
rwth_plt.axis(ax); fig.tight_layout()
```
%% Cell type:markdown id: tags:
## Demonstrator
%% Cell type:code id: tags:
``` python
def decorate_second_axis(ax):
ax.spines['top'].set_position(('axes',.9)); ax.spines['right'].set_color('none'); ax.spines['left'].set_color('none'); ax.spines['bottom'].set_color('none');
ax.xaxis.set_label_coords(.98,1); ax.xaxis.label.set_horizontalalignment('right')
fig, axs = plt.subplots(2, 1, **rwth_plt.landscape); #fig.tight_layout()
highlight_rect_style = {'facecolor': "none", 'edgecolor': 'rwth:orange', 'linestyle':'--', 'linewidth':2.0}
@widgets.interact( r = widgets.IntSlider(min=0, max=10, step=2, value=0, description='$r$ [Hz]'),
Mbr = widgets.IntSlider(min=0, max=10, step=2, value=0, description='$M/r$'))
def update_plots(r, Mbr):
M = r*Mbr
global n, snT, k, SkF
# Continuous functions s(t) and S(f)
s = tri
S = lambda f: si(np.pi*f) ** 2
# Ideal sampling
Sa_f = np.zeros_like(S(f))
nT = []; snT = nT; n = nT;
if r > 0:
T = 1/r
# Sampled function s(n)
nT, snT = sample(t, s(t), T); n = nT/T;
# Sampled continuous spectrum S_a(f)
for k in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
Sa_f += S(f-k/T)
#Sa_f = Sa_f/T # normalize here
# Periodic copies of s(n)
kF = []; SkF = kF; k = kF;
sd = np.zeros_like(snT)
if M > 0:
F = 1/M # in Hz
kF, SkF = sample(f, Sa_f, r*F); k = kF/F/r
for d in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
sd += s((n-d*M)*T)
display('M={}'.format(M))
# 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), 'k', label='$s(t)$');
ax.set_xlabel(r'$\rightarrow t $ [s]')
ax.set_xlim([-5,10]); ax.set_ylim([-.09, 1.19]); ax.legend(loc='upper right', bbox_to_anchor=(1.0, 0.8)); rwth_plt.axis(ax);
axn = ax.twiny(); rwth_plt.stem(axn, [1], [1], color='rwth:blue', label=r'$s(n)$');
axn.set_xlabel(r'$\rightarrow n$'); decorate_second_axis(axn);
axn.set_xlabel(r'$\rightarrow n$ [Samples]'); decorate_second_axis(axn);
axn.legend(loc='upper right', bbox_to_anchor=(1.0, 0.6));
axn.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
ax = axs[1]; ax.set_title('Frequenzbereich');
ax.plot(f, S(f), 'k', label='$S(f)$')
ax.plot(f, Sa_f, 'rwth:red', label='$S_\mathrm{a}(f)/r$')
ax.set_xlabel(r'$\rightarrow f $ [Hz]')
ax.set_xlim([-5,10]); ax.set_ylim([-.09, 1.19]); ax.legend(loc='upper right', bbox_to_anchor=(1.0, 0.8)); rwth_plt.axis(ax);
axk = ax.twiny(); rwth_plt.stem(axk, [1], [1], label=r'$S_\mathrm{d}(k)/r$');
axk.set_xlabel(r'$\rightarrow k$'); decorate_second_axis(axk);
axk.legend(loc='upper right', bbox_to_anchor=(1.0, 0.5));
axk.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
else: # Update only the plot lines themselves. Should not take too much processing time
ax = axs[0] # time domain
axn = fig.axes[2];
if r > 0:
if len(axn.collections) > 1: axn.collections[1].remove()
if M == 0:
rwth_plt.stem_set_data(axn.containers[0], n, snT); #axn.containers[0].set_label(r'$s(n)$')
else:
rwth_plt.stem_set_data(axn.containers[0], n, sd); #axn.containers[0].set_label(r'$s_\mathrm{d}(n)$')
axn.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
ax = axs[1] # frequency domain
ax.lines[1].set_ydata(Sa_f)
axk = fig.axes[3];
if M > 0 and r > 0:
if len(axk.collections) > 1: axk.collections[1].remove()
axk.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
rwth_plt.stem_set_data(axk.containers[0], k, SkF)
axn = fig.axes[2]; axn.set_visible(r > 0)
axk = fig.axes[3]; axk.set_visible(M > 0 and r > 0)
if r > 0:
axn.set_xlim((ax.get_xlim()[0]*r,ax.get_xlim()[1]*r))
if M > 0:
axk.set_xlim((ax.get_xlim()[0]*M/r,ax.get_xlim()[1]*M/r))
```
%% 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, Emin Kosar, Übungsbeispiele zur Vorlesung "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
*Christian Rohlfing, Emin Kosar, Übungsbeispiele zur Vorlesung "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2021, Institut für Nachrichtentechnik, RWTH Aachen University.
......
......@@ -95,7 +95,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
......@@ -109,7 +109,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
"version": "3.9.7"
}
},
"nbformat": 4,
......
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