fft_de.ipynb

{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.7.8","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Visualisieren einer FFT an einer periodischen Schwingung\n\n**allgemeine harmonische Sinusschwingung:**<br><br>\n$\\quad {y(t)= \\hat{y} ⋅ sin(\\omega⋅t + \\phi_0)}$<br> oder <br> $\\quad {y(t)=\\hat{y}⋅sin(2{\\pi}f⋅t+φ_0)}$\n\n$\\hat{y} : Amplitude$ <br>\n${\\omega: Kreisfrequenz}$ in  ${1 \\over s}$<br>\n${f: Frequenz}$ in ${1 \\over s}$;  Periodendauer Ts = ${1 \\over f}$ in s<br>\n${\\phi}_0: Phasenwinkel$<br><br>\n\n\n**Rechteckfunktion:**\n$ x(t) = rect(at) = \\begin{cases}\n    1       & \\quad \\text{fürr |t|} \\leq {1 \\over 2a} \\\\\n    0       & \\quad \\text{sonst }\n  \\end{cases}$\n<br><br>\n[wikipedia, Rechteckfunktion](https://de.wikipedia.org/wiki/Rechteckfunktion)\n\n\n**Dreiecksfunktion:**\n$ x(t) = tri(t) = \\begin{cases}\n    1-\\text{|t|}       & \\quad \\text{für |t| < 1} \\\\\n    0       & \\quad \\text{sonst}\n  \\end{cases}$\n<br><br>\n[wikipedia, Dreiecksfunktion](https://de.wikipedia.org/wiki/Dreiecksfunktion)","metadata":{}},{"cell_type":"code","source":"# resourcen\nimport ipywidgets as widgets\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy import signal\n\n#set backend for interactive toolbar\n%matplotlib widget","metadata":{"trusted":true},"execution_count":1,"outputs":[]},{"cell_type":"code","source":"#presets\nMAX_Amplitude = 5\nMAX_Frequenz = 1000\nMAX_Phase = 180\n\nfNoise = 780\ntimeBase = np.arange(0, 3, 1/10000)     # zeitvektor für  0...3 s, sample rate:  1/fsample (default = 10kHz)","metadata":{"trusted":true},"execution_count":2,"outputs":[]},{"cell_type":"code","source":"# funktionen definieren\n#\ndef getSinus(amp,f,phase,time):\n    return(amp*np.sin(2*np.pi*f*time + phase))\n\n#https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.square.html\ndef getRect(amp,f,time):\n    return(amp*signal.square(2*np.pi*f*time))\n\n#https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.sawtooth.html#scipy.signal.sawtooth\ndef getTri(f,time):\n    return(signal.sawtooth(2*np.pi*f*time, width=0.5))\n\n\ndef calc_FFT(fSample, signal):\n    #https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft\n    # fft mit numpy.fft - routinens\n    #\n    # DFT für signal berechnen\n    fftSig = np.fft.fft(signal)\n    #\n    # anzahl der punkte des komplexen fft-vektors\n    n = len(fftSig)\n    # frequenzen den einzelnen Punkten zuordnen\n    freq = np.fft.fftfreq(n, 1/fSample)\n    #\n    # positiven Bereich des komplexen Spektrums auswaehlen\n    # die Nullfrequenzkomponente in die Mitte des Spektrums verschieben.\n    # - amplituden und frequenzen\n    Y1_shift = np.fft.fftshift(fftSig)\n    F1_shift = np.fft.fftshift(freq)\n    #\n    # indes des Nullpunktes bestimmen\n    #\n    iZero = int(np.ceil(n/2.0))\n    #\n    # amplituden und frequenzen von 0 ... n/2 auswaehlen\n    Y1_pos = Y1_shift[iZero:-1]\n    F1_pos = F1_shift[iZero:-1]\n    #\n    # amplitude normalisieren mit (2* 1/n) und real-Anteil bestimmen\n    #\n    reSpectrum = 2 * 1/n * np.abs(Y1_pos)\n    return(reSpectrum , F1_pos)\n#","metadata":{"trusted":true},"execution_count":3,"outputs":[]},{"cell_type":"code","source":"#interaktive Schalter\n# \n# signal menu\n#\nsignalForms = ['Sinus', 'Rechteck', 'Dreieck']\nsignal_dropdown = widgets.Dropdown(description='Signalform auswählen', options=signalForms, value='Sinus',\\\n                                   style={'description_width': 'initial'})\n#Sliders\namp_1 = widgets.IntSlider(min=0, max=MAX_Amplitude, value=1, description=\"$\\hat{y}_1$:\")\nfrq_1 = widgets.IntSlider(min=0, max=MAX_Frequenz, value=10, step=5, description=\"$f_1$: in Hz\")\nphase_1 = widgets.IntSlider(min=0, max=MAX_Phase, value=0, description='$\\phi$_1: in °')\n#\namp_2 = widgets.IntSlider(min=0, max=MAX_Amplitude, value=1, description=\"$\\hat{y}_2$:\")\nfrq_2 = widgets.IntSlider(min=0, max=MAX_Frequenz, value=10, step=5, description=\"$f_2$: in Hz\")\nphase_2 = widgets.IntSlider(min=0, max=MAX_Phase, value=0, description='$\\phi$_2: in °')\n#\nfSample = widgets.IntSlider(min=100, max=100000, value=10000, step=10, description='Abtastfrequenz',\\\n                            style={'description_width': 'initial'})\n#\n#checkbox for noise overlay\nchkNoise = widgets.Checkbox(description='Störsignal hinzufügen')\n#\n# reset button um auf default-Werte zurückzu setzen\nresetBtn = widgets.Button(description=\"Reset\")","metadata":{"trusted":true},"execution_count":4,"outputs":[]},{"cell_type":"code","source":"# Visualisierung\n# Diagramme zeichnen \nfig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))\n#\n# 1. Callback-Funktion definieren mit der die Figuren aktualisiert werden, wenn sich ein Stueerelement geändert hat\ndef update_view(*args):\n    \n    # Zeitvektor aktualisieren\n    timeBase = np.arange(0, 3, 1/fSample.value)     # zeitvektor für  0...3 s, sample rate:  1/fsample (default = 10kHz)\n    \n    #\n    #Abtastfrequenz muss mindestens doppelt so hoch wie die höchste vorkommende Frequenz sein\n    #(Nyquist - Theorem)\n    # andernfalls: IndexError bei der FFT-Berechnung\n    #\n    if ((frq_1.value > 50) | (frq_2.value > 50)):\n        fSample.min = 2* frq_1.value\n        if (frq_2.value > frq_1.value):\n            fSample.min = 2* frq_2.value\n    \n    if (chkNoise.value == True):        \n        fSample.min = 2* np.max([frq_1.value, frq_2.value, fNoise])\n\n    # dropdowns\n    #\n    signaltype = signal_dropdown.value\n    \n    #-----------------------------------------------------------\n    # Ausgangssignal nach der gewählten Signalform rechnen\n    #\n    if (signaltype == \"Sinus\"):\n        # evtl. ausgeschaltene Steuerlemente wieder aktivieren\n        amp_2.disabled = False\n        frq_2.disabled = False\n        phase_2.disabled = False\n        chkNoise.disabled = False\n        phase_1.disabled = False\n\n        ySin1 = getSinus(amp_1.value, frq_1.value, phase_1.value, timeBase)\n        ySin2 = getSinus(amp_2.value, frq_2.value, phase_2.value, timeBase)\n        ySum = ySin1 + ySin2\n        if (chkNoise.value == True): # das Störsignal hat ein feste Frequenz von 780Hz und eine Amplitude von 5\n            ySum = ySin1 + ySin2 + (5*np.sin(2*np.pi*fNoise*timeBase))\n    \n    if (signaltype == \"Rechteck\"):\n        yRect1 = getRect(amp_1.value, frq_1.value, timeBase)\n        # Steuerelemente des 2. (Sinus)-Signals deaktivieren\n        #\n        amp_2.disabled = True\n        frq_2.disabled = True\n        phase_2.disabled = True\n        chkNoise.disabled = True\n        phase_1.disabled = True\n        #\n        ySum = yRect1\n        #\n    if (signaltype == \"Dreieck\"):\n        yTri = getTri(frq_1.value, timeBase)\n        # Steuerelemente des 2. (Sinus)-Signals deaktivieren        \n        #\n        amp_2.disabled = True\n        frq_2.disabled = True\n        phase_2.disabled = True\n        chkNoise.disabled = True        \n        phase_1.disabled = True\n        #\n        ySum = yTri\n        \n    #-----------------------------------------------------------\n    # FFT berechnen\n    #\n    # Amplitudenspektrum des original-Signals:\n    #\n    fftSig, frqSig = calc_FFT(fSample.value, ySum)\n    #\n    # Rücktransformation\n    # (hier wieder die komplexe DFT verwenden)\n    #\n#     reOrig = np.fft.ifft(np.fft.fft(ySum))\n    \n    #-----------------------------------------------------------\n    # Figuren aktualisiern\n    #\n    # angezeigten Zeitbereich auf 3 signalperioden begrenzen:\n    if (frq_1.value < 2):\n        xEnd = 3\n    if (frq_1.value > 1):\n        xEnd = 3* 1/frq_1.value\n    if (signaltype == 'Sinus'): \n        if ((frq_1.value > 0) & (frq_2.value > 0) & (frq_1.value > frq_2.value)):\n            xEnd = 3* 1/frq_2.value\n    #\n    # Frequenzbereich des FFT-diagrams begrenzen:\n    fftError = True\n    xFFTEnd = len(frqSig)\n    add_f = 0\n    peaks,_ = signal.find_peaks(fftSig, height=(0.01, 10))\n    if (np.count_nonzero(peaks) > 0):\n        fftError = False\n        xFFTEnd = np.max(peaks)\n        add_f = (len(frqSig) - xFFTEnd) -1\n        if add_f > 100:\n            add_f = 100\n    \n    # fig_1: original-signal im zeitbereich\n    #\n    axes[0].clear()\n    axes[0].set_title(\"Ausgangssignal im Zeitbereich\")\n    axes[0].plot(timeBase, ySum, linestyle='-', color='r', label=r\"$y(t)_{{{}}}$\".format(signaltype))\n    if (signaltype == \"Sinus\"):\n        axes[0].plot(timeBase,ySin1, linewidth = 1, linestyle='--', color = 'b', label='$y_{1(t)}$')\n        axes[0].plot(timeBase,ySin2, linewidth = 1, linestyle='dotted', color = '0.5', label='$y_{2(t)}$')\n    axes[0].set_xlabel(\"Zeit /s\")\n    axes[0].set_ylabel(\"Amplitude\")\n    axes[0].set_xlim(timeBase[0], xEnd)\n    axes[0].legend(loc=\"best\")\n    axes[0].grid(True)\n    \n    # fig_2: fft des original-signals\n    #\n    axes[1].clear()\n    if (signaltype == \"Sinus\"):\n        axes[1].set_title(\"FFT des Summensignals\")\n    else:\n        axes[1].set_title(\"FFT des linken Signals\")\n    axes[1].plot(frqSig, fftSig, color='C1')\n    axes[1].set_xlabel(\"f /Hz\")\n    axes[1].set_ylabel(\"Amplitude\")\n    if fftError != True:\n        axes[1].set_xlim(frqSig[0], frqSig[xFFTEnd+add_f])\n    else:\n        yLim = plt.ylim()\n        yPos = yLim[0] + yLim[-1] / 2\n        bbox_props = dict(boxstyle=\"round\", fc=(1, 0.4, 0.4), ec=\"0.5\", alpha=0.9)\n        axes[1].text(frqSig[int(np.ceil(len(frqSig)/2))], yPos, \"Abtastfrequenz beachten\",\\\n                     ha=\"center\", va=\"center\", size=10, bbox=bbox_props)\n\n    axes[1].grid(True)\n    \n    fig.tight_layout()\n    \ndef reset_controls(btn):\n    # Reset Button:\n    fSample.value = 10000\n    signal_dropdown.value = 'Sinus'\n    frq_1.value = 10\n    frq_2.value = 10\n    amp_1.value = 1\n    amp_2.value = 1\n    phase_1.value = 1\n    phase_2.value = 1\n    chkNoise.value = False\n    update_view()\n    \n\n#--------------------------------------------------------\n# 2. die Callback-Funktion mit der Funktion 'observe' den Steuerelementen zuweisen\nsignal_dropdown.observe(update_view, 'value')\nfSample.observe(update_view,'value')\namp_1.observe(update_view,'value')\namp_2.observe(update_view,'value')\nfrq_1.observe(update_view,'value')\nfrq_2.observe(update_view,'value')\nphase_1.observe(update_view,'value')\nphase_2.observe(update_view,'value')\nchkNoise.observe(update_view,'value')\nresetBtn.on_click(reset_controls)\n\n#--------------------------------------------------------\n# Anwendung starten\n#\n# Diagramme einmal zeichnen\nupdate_view()\n#\n#\n# mit 'widgets.VBox / .HBox' die Steuerelemente arangieren\n#\n# box-layout definieren\nbox_layout = widgets.Layout(display='inline-flex',flex_flow='column',align_items='flex-start',width='90%')\n\nwidgets.VBox([widgets.HBox([signal_dropdown, chkNoise]), fSample, widgets.HBox([amp_1, amp_2]),\\\n                            widgets.HBox([frq_1, frq_2]),widgets.HBox([phase_1, phase_2]),\\\n             resetBtn],\\\n            layout = box_layout)","metadata":{"trusted":true},"execution_count":5,"outputs":[{"output_type":"display_data","data":{"text/plain":"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"3dfc6675725a47e6a9f9bb606f874ba5"}},"metadata":{}},{"output_type":"display_data","data":{"text/plain":"VBox(children=(HBox(children=(Dropdown(description='Signalform auswählen', options=('Sinus', 'Rechteck', 'Drei…","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"f0b373b92fe44780a73a20e493bf1974"}},"metadata":{}}]},{"cell_type":"markdown","source":"Copyright © 2020 IUBH Internationale Hochschule","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}