This commit is contained in:
Diana 2024-10-28 15:21:24 +01:00
commit 91f27157dc

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)
@ -27,37 +29,45 @@ axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power [1/Hz]') axs[1].set_ylabel('Power [1/Hz]')
axs[1].set_xlim(0, 1500) axs[1].set_xlim(0, 1500)
# Function to compute and plot the power spectrum
def plot_powerspectrum(i):
# Generate the signal as a sum of two sine waves
def plot_powerspectrum_2(i):
# Clear the previous plots # Clear the previous plots
axs[0].cla() axs[0].cla()
axs[1].cla() axs[1].cla()
# Generate the signal # Generate the signal as a sum of two sine waves
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) + 0.8 * np.sin(2 * np.pi * f2[i] * t) # Second wave is 20% as strong
x[x < 0] = 0 # Apply half-wave rectification
# Plot the signal (first 20 ms for clarity) # Plot the signal (first 20 ms for clarity)
axs[0].plot(t[t < 0.02], x[t < 0.02]) axs[0].plot(t[t < 0.02], x[t < 0.02])
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)
x[x < 0] = 0 # Apply half-wave rectification (optional)
# 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)
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 # Plot the power spectrum
axs[1].plot(freq, power) axs[1].plot(freq, power)
axs[1].set_xlim(0, 1500) axs[1].set_xlim(0, 3000)
axs[1].set_ylim(0, 0.05)
axs[1].set_title(f'Power Spectrum (f2={f2[i]} Hz)') axs[1].set_title(f'Power Spectrum (f2={f2[i]} Hz)')
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]')
#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()