feat: Ch03 Power Spectral Density

Change-Id: I9c58de4ce3fd33db44b2cfbd6f139fd44460dd05
diff --git a/dcs/frames/ch03_psd.py b/dcs/frames/ch03_psd.py
new file mode 100644
index 0000000..975680f
--- /dev/null
+++ b/dcs/frames/ch03_psd.py
@@ -0,0 +1,195 @@
+#  SPDX-License-Identifier: MPL-2.0
+#    Copyright (c) 2022 Philipp Le <philipp@philipple.de>.
+#  This Source Code Form is subject to the terms of the Mozilla Public
+#  License, v. 2.0. If a copy of the MPL was not distributed with this
+#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+from __future__ import annotations
+
+from tkinter import ttk, LEFT, BOTH
+
+import scipy.signal.windows
+from pydantic import confloat, conint
+from dcs.config import default_store, ConfigObject, ui_create, ConfigControlFrame
+import numpy as np
+from dcs.frames.base import BaseFrame, Window
+from dcs.frames.groups import Ch03Group
+from dcs.utils import swap_freq
+from typing import List
+from enum import Enum
+
+import matplotlib
+matplotlib.use('TkAgg')
+from matplotlib.figure import Figure
+from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
+
+
+class ConfigCh03PsdFrame(ConfigControlFrame):
+    wdg_signal_rms: ttk.Label
+    wdg_noise_rms: ttk.Label
+
+
+class DisplayMode(str, Enum):
+    Lin = 'linear'
+    Log = 'logarithmic'
+
+
+@ui_create
+class ConfigCh03Psd(ConfigObject):
+    _KEY = 'ch03_psd'
+
+    amplitude_dBV: confloat(ge=-160.0, lt=40.0, multiple_of=0.1) = 5.0
+    freq_mult: conint(ge=0, lt=32) = 16
+    noise_dBV: confloat(ge=-160.0, le=40.0, multiple_of=0.1) = -10.0
+    random_seed: conint(ge=0, lt=10000) = 1000
+
+    disp_mode: DisplayMode = DisplayMode.Log
+
+    def make_config_widget(self, parent: ttk.Widget) -> ConfigControlFrame:
+        frm = ConfigCh03PsdFrame(parent, borderwidth=1, relief='raised')
+
+        ttk.Label(frm, text='Signal Level:').grid(row=0, column=0)
+        w = self.ui_create_amplitude_dBV(frm)
+        w.grid(row=0, column=1)
+        frm.add_widget(w)
+        ttk.Label(frm, text='dBV').grid(row=0, column=2)
+
+        ttk.Label(frm, text='Signal RMS:').grid(row=1, column=0)
+        frm.wdg_signal_rms = ttk.Label(frm, text='')
+        frm.wdg_signal_rms.grid(row=1, column=1)
+        ttk.Label(frm, text='V').grid(row=1, column=2)
+
+        ttk.Label(frm, text='Frequency:').grid(row=2, column=0)
+        w = self.ui_create_freq_mult(frm)
+        w.grid(row=2, column=1)
+        frm.add_widget(w)
+
+        ttk.Label(frm, text='Noise Level:').grid(row=3, column=0)
+        w = self.ui_create_noise_dBV(frm)
+        w.grid(row=3, column=1)
+        frm.add_widget(w)
+        ttk.Label(frm, text='dBV').grid(row=3, column=2)
+
+        ttk.Label(frm, text='Noise RMS:').grid(row=4, column=0)
+        frm.wdg_noise_rms = ttk.Label(frm, text='')
+        frm.wdg_noise_rms.grid(row=4, column=1)
+        ttk.Label(frm, text='V').grid(row=4, column=2)
+
+        ttk.Label(frm, text='Random Seed:').grid(row=5, column=0)
+        w = self.ui_create_random_seed(frm)
+        w.grid(row=5, column=1)
+        frm.add_widget(w)
+
+        ttk.Label(frm, text='Display Mode:').grid(row=6, column=0)
+        w = self.ui_create_disp_mode_dropdown(frm)
+        w.grid(row=6, column=1)
+        frm.add_widget(w)
+
+        return frm
+
+    def calc_signal_rms(self):
+        return np.power(10, (self.amplitude_dBV / 20))
+
+    def calc_noise_rms(self):
+        return np.power(10, (self.noise_dBV / 20))
+
+    def _calc_noise_sigma(self):
+        return np.sqrt(2) * self.calc_noise_rms()
+
+    def calc_signal(self, t: np.ndarray) -> np.ndarray:
+        phi = np.sqrt(2) * self.calc_signal_rms() * np.exp(1j * 2 * np.pi * (self.freq_mult) * t)
+        return np.real(phi) + np.random.normal(0, self._calc_noise_sigma(), len(t))
+
+
+class Ch03PsdFrame(BaseFrame):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+        self._config: ConfigCh03Psd = default_store().get_config(ConfigCh03Psd)
+
+        self.ctrl_frm = self._create_control()
+        self.ctrl_frm.pack(side=LEFT)
+
+        signal_frm = self._create_signal_canvas()
+        signal_frm.pack(expand=True, fill=BOTH)
+
+    def _create_control(self) -> ConfigCh03PsdFrame:
+        frm: ConfigCh03PsdFrame = self._config.make_config_widget(self)
+        frm.widgets_on_change(self._on_change)
+        return frm
+
+    def _on_change(self, _, __, ___):
+        default_store().save()
+
+        self.draw_signal()
+
+    def _create_signal_canvas(self) -> ttk.Widget:
+        frm = ttk.Frame(self)
+
+        self._signal_fig = Figure(figsize=(12, 6), dpi=100)
+
+        self._signal_canvas = FigureCanvasTkAgg(self._signal_fig, frm)
+        self._signal_canvas.get_tk_widget().pack(expand=True, fill=BOTH)
+
+        self.draw_signal()
+
+        return frm
+
+    def draw_signal(self):
+        LENGTH = 512
+
+        self.ctrl_frm.wdg_signal_rms.configure(text=f'{self._config.calc_signal_rms()}')
+        self.ctrl_frm.wdg_noise_rms.configure(text=f'{self._config.calc_noise_rms()}')
+
+        np.random.seed(self._config.random_seed)
+
+        self._signal_fig.clear()
+
+        ax_signal = self._signal_fig.add_subplot(2, 1, 1)
+        ax_psd = self._signal_fig.add_subplot(2, 1, 2)
+
+        ax_signal.set_xlim(-1.2, 1.2)
+        ax_signal.set_xlabel('time')
+        ax_signal.set_ylabel('value')
+        ax_signal.set_title('Signal')
+
+        #ax_psd.set_xlim(0, 512)
+        ax_psd.set_xlabel('frequency')
+        ax_psd.set_ylabel('power density')
+        ax_psd.set_title('Power Spectral Density')
+
+        # t = np.linspace(-1.0, 1.0, LENGTH)
+        t = np.arange(-256, 256, 1) / 256.0
+        # f = np.arange(-int(LENGTH/2), int(LENGTH/2), 1)
+        f = swap_freq(np.fft.fftfreq(t.shape[-1], 2.0/LENGTH))
+
+        sig = self._config.calc_signal(t)
+        ax_signal.plot(t, sig, label='Signal', linestyle='solid', linewidth=1)
+        ax_signal.legend()
+
+        #psd = np.abs(np.fft.fft(np.correlate(sig, sig, mode='same')))
+        #psd = np.abs(self._swap(np.fft.fft(sig * scipy.signal.windows.hamming(len(sig)))))
+
+        windowed = sig  # * scipy.signal.windows.hamming(len(sig))
+        fftransformed = swap_freq(np.fft.fft(windowed)) / len(sig)
+        psd = np.abs(fftransformed)
+        if self._config.disp_mode == DisplayMode.Log:
+            psd = 20 * np.log10(psd)
+            ax_psd.set_ylabel('power density (dbV^2)')
+        else:
+            ax_psd.set_ylabel('power density (V^2)')
+        ax_psd.plot(f, psd, label='PSD', linestyle='solid', linewidth=1)
+        ax_psd.legend()
+
+        self._signal_canvas.draw()
+
+
+class Ch03PsdWindow(Window):
+    GROUP = Ch03Group
+    TITLE = 'Power Spectral Density'
+    FRAME = Ch03PsdFrame
+
+
+if __name__ == '__main__':
+    Ch03PsdWindow.main()
+