Add Chua circuit examples, alpha-transform dynamics, and core WDF extensions

README.md

@@ -1,10 +1,22 @@
 # pywdf
 
+
 [![General badge](https://img.shields.io/badge/PDF-Paper-<COLOR>.svg)](https://repositori.upf.edu/handle/10230/57903)
 
 
 <code>pywdf</code> is a Python framework for modeling and simulating wave digital filter circuits. It allows users to easily create and analyze WDF circuit models in a high-level, object-oriented manner. The library includes a variety of built-in components, such as voltage sources, capacitors, diodes etc., as well as the ability to create custom components and circuits. Additionally, pywdf includes a variety of analysis tools, such as frequency response and transient analysis, to aid in the design and optimization of WDF circuits. Also included are several example circuits as shown below. 
 
+
+## About *frantic0's* pywdf fork≈
+
+This fork adds new examples of the Chua circuit that you can find in the structure below and a few extensions to the core *wdf.py*, including the following: 
+* a new element, `ChuaDiode`, a non-linear negative resistor that implements a piecewise-linear I-V relationship,  
+* a new two-port series adaptor with an embedded resistive voltage source, `SeriesVoltage` for exciting circuits with voltage impulses injection, useful for the Chua circuit,
+* an implementation of the alpha transform for digitizing linear dynamic elements and testing different stability levels, with optional parameters for choosing between the bilinear and backward Euler transforms. 
+
+
+
+
 ## Installation
 ```
 pip install git+https://github.com/gusanthon/pywdf
@@ -23,6 +35,9 @@ The <code>core</code> directory contains the main source code of the repository.
 │       ├── baxandalleq.py
 │       ├── diodeclipper.py
 │       ├── lc_oscillator.py
+│       ├── chua.py					implemented by @frantic0
+│       ├── chua_minimal.py			implemented by @frantic0
+│       ├── chua_ODE.py				implemented by @frantic0
 │       ├── passivelpf.py
 │       ├── rca_mk2_sef.py
 │       ├── rclowpass.py
@@ -33,6 +48,18 @@ The <code>core</code> directory contains the main source code of the repository.
 ├── setup.py
 ```
 
+## Interactive Debugging
+
+
+```
+python3 -m debugpy --wait-for-client --listen 5678 pywdf/examples/chua_circuit.py
+```
+
+
+
+
+
+
 ## Usage
 
 ```python
@@ -96,6 +123,7 @@ For further reading, check out:
 - Giovanni De Sanctis and Augusto Sarti, “Virtual analog modeling in the wave-digital domain,” IEEE transactions on audio, speech, and language processing, vol. 18, no. 4, pp. 715–727, 2009. [[URL]](https://ieeexplore.ieee.org/abstract/document/5276845)
 - Kurt James Werner, Vaibhav Nangia, Julius O Smith, and Jonathan S Abel, “A general and explicit formulation for wave digital filters with multiple/multiport nonlinearities and complicated topologies,” IEEE, 2015, pp. 1–5. [[URL]](https://ieeexplore.ieee.org/document/7336908)
 - D. Franken, Jörg Ochs, and Karlheinz Ochs, “Generation of wave digital structures for networks containing multiport elements,” Circuits and Systems I: Regular Papers, IEEE Transactions on, vol. 52, pp. 586 – 596, 04 2005. [[URL]](https://www.researchgate.net/publication/4018571_Generation_of_wave_digital_structures_for_connection_networks_containing_ideal_transformers)
+- Klaus Meerkotter and Rheinhard Scholz, “Digital simulation of nonlinear circuits by wave digital filter principles”, 1989, IEEE International Symposium on Circuits and Systems (ISCAS)
 
 
 ## Contributions

pywdf/__init__.py

@@ -1,13 +1,14 @@
 from .core.circuit import *
 from .core.rtype import *
 from .core.wdf import *
+from .core.solver import *
 
 from .examples.tr_808_hatresonator import TR_808_HatResonator
 from .examples.bassmantonestack import BassmanToneStack
 from .examples.baxandalleq import BaxandallEQ, UnadaptedBaxandallEQ
 from .examples.diodeclipper import DiodeClipper
 from .examples.lc_oscillator import LCOscillator
-from .examples.passivelpf import PassiveLPF
+from .examples.passive_lpf import PassiveLPF
 from .examples.rca_mk2_sef import RCA_MK2_SEF
-from .examples.rclowpass import RCLowPass
-from .examples.voltagedivider import VoltageDivider
+from .examples.rc_lowpass import RCLowPass
+from .examples.voltage_divider import VoltageDivider

pywdf/core/circuit.py

@@ -22,6 +22,7 @@ def __init__(self, source: baseWDF, root: rootWDF, output: baseWDF) -> None:
         self.root = root
         self.output = output
 
+
     def process_sample(self, sample: float) -> float:
         """Process an individual sample with this circuit.
 
@@ -38,6 +39,40 @@ def process_sample(self, sample: float) -> float:
         self.root.next.accept_incident_wave(self.root.propagate_reflected_wave())
         return self.output.wave_to_voltage()
 
+
+
+    def process_sample_i_v(self, sample: float) -> float:
+        """Process an individual sample with this circuit.
+
+        Note: not every circuit will follow this general pattern, in such cases users may want to overwrite this function. See example circuits
+
+        Args:
+            sample (float): incoming sample to process
+
+        Returns:
+            (i, v) I-V tupple: processed sample
+        """
+        self.source.set_voltage(sample)
+        self.root.accept_incident_wave(self.root.next.propagate_reflected_wave())
+        self.root.next.accept_incident_wave(self.root.propagate_reflected_wave())
+
+        return ( self.output.wave_to_voltage(), self.source.wave_to_current(), self.output.wave_to_current() ) 
+
+
+    def process_i_v_signals(self, signal: np.array) -> np.array:
+        """Process an entire signal with this circuit.
+
+        Args:
+            signal (np.array): incoming signal to process
+
+        Returns:
+           (i,v) tuples list: processed signal
+        """
+        self.reset()
+        l = [ self.process_sample_i_v(sample) for sample in signal]
+        return l
+
+
     def process_signal(self, signal: np.array) -> np.array:
         """Process an entire signal with this circuit.
 
@@ -48,7 +83,10 @@ def process_signal(self, signal: np.array) -> np.array:
             np.array: processed signal
         """
         self.reset()
-        return np.array([self.process_sample(sample) for sample in signal])
+
+        return np.array([ self.process_sample(sample) for sample in signal ])
+
+
 
     def process_wav(self, filepath: str, output_filepath: str = None) -> np.array:
         fs, x = wavfile.read(filepath)
@@ -67,8 +105,14 @@ def __call__(self, *args: any, **kwds: any) -> any:
         elif hasattr(args[0], "__iter__"):
             return self.process_signal(args[0])
 
+
+
+
+
+
     def get_impulse_response(self, delta_dur: float = 1, amp: float = 1) -> np.array:
-        """Get circuit's impulse response
+        """
+        Get circuit's impulse response
 
         Args:
             delta_dur (float, optional): duration of Dirac delta function in seconds. Defaults to 1.
@@ -81,6 +125,34 @@ def get_impulse_response(self, delta_dur: float = 1, amp: float = 1) -> np.array
         d[0] = amp
         return self.process_signal(d)
 
+
+
+
+
+    def plot_signal(
+        self,
+        signal: np.array,
+        n_samples: int = 500,
+        outpath: str = None,
+        title: str = "Signal",
+    ) -> None:
+
+        plt.figure(figsize=(9, 5.85))
+        plt.plot(signal[:n_samples])
+        plt.xlabel("Sample [n]")
+        plt.ylabel("Amplitude [V]")
+        plt.title(title)
+
+        plt.grid()
+        if outpath:
+            plt.savefig(outpath)
+        plt.show()
+
+
+
+
+
+
     def plot_impulse_response(
         self,
         n_samples: int = 500,
@@ -120,14 +192,23 @@ def reset(self) -> None:
             if isinstance(self.__dict__[key], baseWDF):
                 self.__dict__[key].reset()
 
+
+
+
+
     def compute_spectrum(self, fft_size: int = None) -> np.ndarray:
         x = self.get_impulse_response()
         N2 = int(fft_size / 2 - 1)
         H = fft(x, fft_size)[:N2]
         return H
 
+
+
+
+
     def plot_freqz(self, outpath: str = None, fft_size: int = None):
-        """Plot the circuit's frequency response
+        """
+        Plot the circuit's frequency response
 
         Args:
             outpath (str, optional): filepath to save figure. Defaults to None.
@@ -171,6 +252,9 @@ def plot_freqz(self, outpath: str = None, fft_size: int = None):
             plt.savefig(outpath)
         plt.show()
 
+
+
+
     def plot_freqz_list(
         self,
         values: list,
@@ -241,6 +325,59 @@ def plot_freqz_list(
 
         plt.show()
 
+
+
+
+
+    def i_v_analysis(
+        self,
+        freq: float = 1000,
+        amplitude: float = 1,
+        t_ms: float = 5,
+        outpath: str = None,
+    ):
+        """Plot transient analysis of Circuit's response to sine wave
+
+        Args:
+            freq (float, optional): frequency of sine wave. Defaults to 1000.
+            amplitude (float, optional): amplitude of sine wave. Defaults to 1.
+            t_ms (float, optional): time in ms of sine wave. Defaults to 5.
+        """
+        _, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 6.5))
+
+        n_samples = int(t_ms * self.fs / 1000)
+        n = np.arange(0, 2, 1 / self.fs)
+        x = np.sin(2 * np.pi * freq * n) * amplitude
+
+        y = self.process_i_v_signals(x)
+
+        v, ii, io = zip(*y)
+
+        v, ii, io = np.array(v), np.array(ii) * 100, np.array(io) * 0.01 
+
+        ax.plot(x[:n_samples], label="input signal")
+        ax.plot(v[:n_samples], label="voltage out", alpha=0.75)
+        ax.plot(ii[:n_samples], label="current source through", alpha=0.75)
+        ax.plot(io[:n_samples], label="current output through", alpha=0.75)
+        ax.set_xlabel("sample")
+        ax.set_ylabel("amplitude")
+        ax.set_yscale('log')
+        ax.set_ylim(0, 1000000)
+
+        ax.set_title(loc="left", label="output signal vs input signal waveforms")
+        ax.grid(True)
+        ax.legend()
+
+        if outpath:
+            plt.savefig(outpath)
+
+        plt.show()
+
+
+
+
+
+
     def AC_transient_analysis(
         self,
         freq: float = 1000,

pywdf/core/solver/__init__.py

@@ -0,0 +1 @@
+__all__ = ['runge_kutta']

pywdf/core/solver/runge_kutta.py

@@ -0,0 +1,57 @@
+import numpy as np
+
+def rk4_step(f, t, y, h):
+    """
+    Single step of fourth-order Runge-Kutta method
+
+    Parameters:
+    f: function that defines dy/dt = f(t, y)
+    t: current time
+    y: current value of y
+    h: step size
+
+    Returns:
+    y_next: next value of y
+    """
+    k1 = h * f(t, y)
+    k2 = h * f(t + h/2, y + k1/2)
+    k3 = h * f(t + h/2, y + k2/2)
+    k4 = h * f(t + h, y + k3)
+
+    y_next = y + (k1 + 2*k2 + 2*k3 + k4) / 6
+    return y_next
+
+def rk4_solve(f, t_span, y0, h):
+    """
+    Solve ODE using RK4 method
+
+    Parameters:
+    f: function that defines dy/dt = f(t, y)
+    t_span: tuple (t_start, t_end)
+    y0: initial condition (scalar or array)
+    h: step size
+
+    Returns:
+    t_values: array of time points
+    y_values: array of solution values
+    """
+    t_start, t_end = t_span
+    t_values = np.arange(t_start, t_end + h, h)
+
+    # Handle both scalar and vector initial conditions
+    y0 = np.asarray(y0)
+    if y0.ndim == 0:  # Scalar case
+        y_values = np.zeros(len(t_values))
+    else:  # Vector case
+        y_values = np.zeros((len(t_values), len(y0)))
+
+    # Initial condition
+    y_values[0] = y0
+
+    # Apply RK4 method
+    for i in range(len(t_values) - 1):
+        y_values[i + 1] = rk4_step(f, t_values[i], y_values[i], h)
+
+    return t_values, y_values
+
+