64 lines
1.7 KiB
Python
64 lines
1.7 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Tue Oct 15 16:46:07 2024
|
|
|
|
@author: ernesto311
|
|
"""
|
|
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from scipy.signal import welch
|
|
from matplotlib.animation import FuncAnimation, PillowWriter
|
|
|
|
# Generate distances and corresponding frequencies
|
|
distances = np.arange(-400, 451, 1)
|
|
f1 = 800
|
|
f2 = f1 + distances
|
|
|
|
# Time parameters
|
|
dt = 0.00001
|
|
t = np.arange(0, 2, dt)
|
|
|
|
# Set up the figure and axis for the animation
|
|
fig, axs = plt.subplots(1, 2, layout='constrained')
|
|
axs[0].set_xlabel('Time [s]')
|
|
axs[0].set_ylabel('Amplitude')
|
|
axs[1].set_xlabel('Frequency [Hz]')
|
|
axs[1].set_ylabel('Power [1/Hz]')
|
|
axs[1].set_xlim(0, 1500)
|
|
|
|
# Function to compute and plot the power spectrum
|
|
def plot_powerspectrum(i):
|
|
# Clear the previous plots
|
|
axs[0].cla()
|
|
axs[1].cla()
|
|
|
|
# Generate the signal
|
|
x = np.sin(2*np.pi*f1*t) + 0.2 * np.sin(2*np.pi*f2[i]*t)
|
|
x[x < 0] = 0 # Apply half-wave rectification
|
|
|
|
# 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(0, 1.2)
|
|
|
|
# Compute power spectrum
|
|
freq, power = welch(x, fs=1/dt, nperseg=2**16)
|
|
|
|
# Plot the power spectrum
|
|
axs[1].plot(freq, power)
|
|
axs[1].set_xlim(0, 1500)
|
|
axs[1].set_ylim(0, 0.05)
|
|
axs[1].set_title(f'Power Spectrum (f2={f2[i]} Hz)')
|
|
axs[1].set_xlabel('Frequency [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()
|