Various sims

dft.py

@@ -0,0 +1,24 @@
+import numpy as np
+
+def DFT(x):
+    """
+    Compute the Discrete Fourier Transform (DFT) of the input signal x.
+    
+    Parameters:
+        x (array-like): Input signal
+        
+    Returns:
+        X (array-like): DFT of the input signal
+    """
+    N = len(x)
+    n = np.arange(N)
+    k = n.reshape((N, 1))
+    e = np.exp(-2j * np.pi * k * n / N)
+    X = np.dot(e, x)
+    return X
+
+# Example usage:
+x = np.array([0, 1, 2, 3])
+X = DFT(x)
+print("DFT of", x, ":", X)
+

dft2.py

@@ -0,0 +1,74 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def DFT(x):
+    """
+    Compute the Discrete Fourier Transform (DFT) of the input signal x.
+    
+    Parameters:
+        x (array-like): Input signal
+        
+    Returns:
+        X (array-like): DFT of the input signal
+    """
+    N = len(x)
+    n = np.arange(N)
+    k = n.reshape((N, 1))
+    e = np.exp(-2j * np.pi * k * n / N)
+    X = np.dot(e, x)
+    return X
+
+def DFTI(X):
+    """
+    Compute the Inverse Discrete Fourier Transform (IDFT) of the input signal X.
+    
+    Parameters:
+        X (array-like): Input signal
+        
+    Returns:
+        x (array-like): IDFT of the input signal
+    """
+    N = len(X)
+    n = np.arange(N)
+    k = n.reshape((N, 1))
+    e = np.exp(2j * np.pi * k * n / N)
+    x = np.dot(e, X) / N
+    return x
+
+# Define a sine function as the input signal
+t = np.linspace(0, 1, 1000)  # Time array
+frequency = 5  # Frequency of the sine wave
+x = np.sin(2 * np.pi * frequency * t) + np.sin(np.pi * frequency * t)
+
+# Compute the DFT of the sine wave
+X = DFT(x)
+
+# Plot the input signal
+plt.figure(figsize=(12, 4))
+plt.subplot(1, 3, 1)
+plt.plot(t, x)
+plt.title('Input Signal')
+plt.xlabel('Time')
+plt.ylabel('Amplitude')
+
+# Plot the magnitude spectrum (DFT)
+freq = np.fft.fftfreq(len(x), d=t[1]-t[0])
+plt.subplot(1, 3, 2)
+plt.stem(freq, np.abs(X))
+plt.title('Magnitude Spectrum (DFT)')
+plt.xlabel('Frequency')
+plt.ylabel('Magnitude')
+plt.xlim(-10, 10)  # Limit the x-axis for better visualization
+
+# Compute the IDFT of the DFT signal
+x_reconstructed = DFTI(X)
+# Plot the reconstructed signal
+plt.subplot(1, 3, 3)
+plt.plot(t, x_reconstructed.real, label='Reconstructed Signal')
+plt.plot(t, x, label='Original Signal')
+plt.title('Reconstructed Signal vs Original Signal')
+plt.xlabel('Time')
+plt.ylabel('Amplitude')
+plt.legend()
+plt.tight_layout()
+plt.show()

pid_sim.py

@@ -0,0 +1,75 @@
+# Simulation of a PID-controller
+import matplotlib.pyplot as plt
+
+class PIDController:
+    def __init__(self, Kp, Ki, Kd, setpoint):
+        self.Kp = Kp
+        self.Ki = Ki
+        self.Kd = Kd
+        self.setpoint = setpoint
+        self.prev_error = 0
+        self.integral = 0
+
+    def update(self, measured_value):
+        error = self.setpoint - measured_value
+        self.integral += error
+        derivative = error - self.prev_error
+        output = self.Kp * error + self.Ki * self.integral + self.Kd * derivative
+        self.prev_error = error
+        return output
+    
+    def set_setpoint(self, setpoint):
+        self.setpoint = setpoint
+
+# Simulation parameters
+Kp = 4 # Proportional gain
+Ki = 1/4 # Integral gain
+Kd = 1/2 # Derivative gain
+setpoint = 100.0
+initial_value = 0.0
+simulation_time = 100
+time_step = 0.1
+
+# Create PID controller
+pid_controller = PIDController(Kp, Ki, Kd, setpoint)
+
+# Initialize lists to store data for plotting
+time_points = []
+control_signals = []
+current_values = []
+
+# Simulate PID controller
+current_value = initial_value
+for t in range(simulation_time):
+    if t == simulation_time // 3:
+        pid_controller.set_setpoint(50.0)
+    if t == 2 * simulation_time // 3:
+        pid_controller.set_setpoint(150.0)
+    time_points.append(t)
+    control_signal = pid_controller.update(current_value)
+    control_signals.append(control_signal)
+    # Apply control signal (simulate process)
+    current_value += control_signal * time_step
+    current_values.append(current_value)
+
+# Plotting
+plt.figure(figsize=(10, 6))
+plt.subplot(2, 1, 1)
+plt.plot(time_points, current_values, label='Current Value')
+plt.plot([0, simulation_time // 3], [100, 100], 'r--', label='Setpoint 100')
+plt.plot([simulation_time // 3, 2 * simulation_time // 3], [50, 50], 'r--', label='Setpoint 50')
+plt.plot([2 * simulation_time // 3, simulation_time], [150, 150], 'r--', label='Setpoint 150')
+plt.xlabel('Time')
+plt.ylabel('Value')
+plt.title('PID Controller Simulation')
+plt.legend()
+
+plt.subplot(2, 1, 2)
+plt.plot(time_points, control_signals, 'g', label='Control Signal')
+plt.xlabel('Time')
+plt.ylabel('Control Signal')
+plt.legend()
+
+plt.tight_layout()
+# plt.savefig('pid_simulation.png') # Save to file
+plt.show() # Show immediately

statpy.py

@@ -0,0 +1,107 @@
+#!/bin/python3
+import numpy as np
+
+# import pandas as pd
+import time
+
+# Print the current time
+print(time.strftime("%Y-%m-%d %H:%M:%S"))
+
+a = {1, 2, 3}
+assert 1 in a, "1 is not in a"
+assert 4 not in a, "4 is in a"
+
+b = {"a": 1, "b": 2, "c": 3}
+assert "a" in b, "a is not in b"
+assert "d" not in b, "d is in b"
+
+c = a.union(b.values())
+assert 1 in c, "1 is not in c"
+
+# Xor
+assert 1 ^ 1 == 0
+assert 1 ^ 0 == 1
+assert 0 ^ 1 == 1
+assert 0 ^ 0 == 0
+
+a = 0b1010  # 10
+a = a ^ 0b1010  # 10
+assert a == 0, "a should be 0"
+
+# Poisson distribution
+# print(np.random.poisson(5, 1))
+# print(np.random.binomial(10, 0.5, 10))
+
+
+# Probability mass function of poisson distribution
+def pmf_poisson(k, l):
+    return (l**k) * np.exp(-l) / np.math.factorial(k)
+
+
+def mean(l):
+    return sum(l) / len(l)
+
+
+def median(l):
+    l.sort()
+    n = len(l)
+    if n % 2 == 0:
+        return (l[n // 2 - 1] + l[n // 2]) / 2
+    else:
+        return l[n // 2]
+
+def mode(l):
+    return max(set(l), key=l.count)
+
+
+assert median([1, 2, 3, 4]) == 2.5
+assert median([1, 2, 3, 4, 5]) == 3
+
+quantile = np.quantile
+# def quantile(l, q):
+
+assert quantile([1, 2, 3, 4, 5], 0.5) == 3
+assert quantile([1, 2, 3, 4, 5], 0.25) == 2
+assert quantile([1, 2, 3, 4, 5], 0.75) == 4
+
+
+def whisker_plot(l):
+    q1 = quantile(l, 0.25)  # Lower quartile
+    q3 = quantile(l, 0.75)  # Upper quartile
+    iqr = q3 - q1  # Interquartile range
+    lw = q1 - 1.5 * iqr  # Lower whisker
+    uw = q3 + 1.5 * iqr  # Upper whisker
+    return lw, q1, q3, uw
+
+
+def outliers(l):
+    lw, q1, q3, uw = whisker_plot(l)
+    return [x for x in l if x < lw or x > uw]
+
+
+values = [16, 18, 28, 13, 50, 31, 25, 22, 18, 23, 29, 38]
+
+assert median(values) == 24
+# assert quantile(values, 0.25) == 18
+# assert quantile(values, 0.75) == 30
+
+print(quantile(values, 0.25))
+print(quantile(values, 0.75))
+
+
+print(whisker_plot(values))
+print(outliers(values))
+
+# Plot the box diagram
+import matplotlib.pyplot as plt
+
+# plt.boxplot(values, vert=False)
+# plt.show()
+
+
+# Plot the poisson distribution when lambda = 5
+# import matplotlib.pyplot as plt
+
+# s = np.random.poisson(3, 10000)
+# count, bins, ignored = plt.hist(s, 10, density=True)
+# plt.show()