[Animation_ChatGPT] updated

This commit is contained in:
mbergmann 2024-10-28 15:11:01 +01:00
parent 851857c19d
commit 2515472d32

View File

@ -9,12 +9,14 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy.signal import welch from scipy.signal import welch
from matplotlib.animation import FuncAnimation, PillowWriter from matplotlib.animation import FuncAnimation, PillowWriter
import useful_functions as f
# Generate distances and corresponding frequencies # Generate distances and corresponding frequencies
distances = np.arange(-400, 451, 1) distances = np.arange(-400, 2000, 1)
f1 = 800 f1 = 800
f2 = f1 + distances f2 = f1 + distances
# Time parameters # Time parameters
dt = 0.00001 dt = 0.00001
t = np.arange(0, 2, dt) t = np.arange(0, 2, dt)
@ -34,7 +36,7 @@ def plot_powerspectrum(i):
axs[1].cla() axs[1].cla()
# Generate the signal # Generate the signal
x = np.sin(2*np.pi*f1*t) + 0.2 * np.sin(2*np.pi*f2[i]*t) x = np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2[i]*t)
x[x < 0] = 0 # Apply half-wave rectification x[x < 0] = 0 # Apply half-wave rectification
# Plot the signal (first 20 ms for clarity) # Plot the signal (first 20 ms for clarity)
@ -42,7 +44,7 @@ def plot_powerspectrum(i):
axs[0].set_title(f"Signal (f2={f2[i]} Hz)") axs[0].set_title(f"Signal (f2={f2[i]} Hz)")
axs[0].set_xlabel('Time [s]') axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Amplitude') axs[0].set_ylabel('Amplitude')
axs[0].set_ylim(0, 1.2) axs[0].set_ylim(-2, 2)
# Compute power spectrum # Compute power spectrum
freq, power = welch(x, fs=1/dt, nperseg=2**16) freq, power = welch(x, fs=1/dt, nperseg=2**16)
@ -55,9 +57,50 @@ def plot_powerspectrum(i):
axs[1].set_xlabel('Frequency [Hz]') axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power [1/Hz]') axs[1].set_ylabel('Power [1/Hz]')
# # Create the animation
# ani = FuncAnimation(fig, plot_powerspectrum, frames=len(distances), interval=500)
# # Display the animation
# ani.save("signal_animation.gif", writer=PillowWriter(fps=30))
# plt.show()
# Generate the signal as a sum of two sine waves
def plot_powerspectrum_2(i):
# Clear the previous plots
axs[0].cla()
axs[1].cla()
# Generate the signal as a sum of two sine waves
x = np.sin(2 * np.pi * f1 * t) + 0.8 * np.sin(2 * np.pi * f2[i] * t) # Second wave is 20% as strong
# Plot the signal (first 20 ms for clarity)
axs[0].plot(t[t < 0.02], x[t < 0.02])
axs[0].set_title(f"Signal (f2={f2[i]} Hz)")
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Amplitude')
axs[0].set_ylim(-2, 2)
x[x < 0] = 0 # Apply half-wave rectification (optional)
# Compute power spectrum
freq, power = welch(x, fs=1/dt, nperseg=2**16)
pref = np.max(power)
decibel_power = 10 * np.log10(power/pref)
AM = f.find_AM(f1, 0.5 * f1, f2[i])
# Plot the power spectrum
axs[1].plot(freq, power)
axs[1].set_xlim(0, 3000)
axs[1].set_title(f'Power Spectrum (f2={f2[i]} Hz)')
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power [1/Hz]')
axs[1].set_ylim(0, 0.00007)
axs[1].plot(f1, power[np.argmin(np.abs(freq-f1))], 'o')
axs[1].plot(f2[i], power[np.argmin(np.abs(freq-f2[i]))], 'd')
axs[1].plot(AM, power[np.argmin(np.abs(freq-AM))], '*')
axs[1].axvline(AM, alpha = 0.5, color = 'r')
# Create the animation # Create the animation
ani = FuncAnimation(fig, plot_powerspectrum, frames=len(distances), interval=500) ani = FuncAnimation(fig, plot_powerspectrum_2, frames=len(distances), interval=500)
# Display the animation # Save the animation as a GIF file (optional)
ani.save("signal_animation.gif", writer=PillowWriter(fps=30)) ani.save("sum_of_sinewaves.gif", writer=PillowWriter(fps=30))
plt.show()