[Python] How to reduce deviation when calculating quadratic value from interpolated data set?

I have written a Python program that calculates the received signal strength of an antenna communicating with a moving transmitter. The trajectory data for the transmitter isn't evenly spaced, and project requirements require me to interpolate it out. The program is written such that the user can choose the sample rate of the interpolation.

There is a squared term in the signal strength formula that is calculated from the transmission position. I have noticed that this term varies substantially when the sample rate is varied, which I assume is because small deviations in the interpolation lead to exponential changes in this term.

I don't have much data science/scipy experience, so I am not sure how to smooth this data out - frankly, I'm not even sure what keywords to Google for this.

My interpolation code is as follows - note that the unique time stuff is because the user also chooses the time slice of the trajectory that they want to calculate to improve performance:

def interp_position(times, sample_rate, position_vector, start_time=0, end_time= None, coord_system = "spherical",new_times=None):
#times with parameterized sample rate
if end_time is None:
end_time = times[-1]
if new_times is None:
new_times = np.arange(times[0], times[-1],1/sample_rate)

position_vector[:, 1] = np.radians(position_vector[:, 1])
position_vector[:, 2] = np.radians(position_vector[:, 2])

# Find the indices corresponding to the time interval
unique_indices = np.unique(times, return_index=True)[1]
times = times[unique_indices]
position_vector = position_vector[unique_indices]

interval_indices = np.where((new_times >= start_time) & (new_times <= end_time))[0]


# Interpolate position functions for the specified time interval
pos_interp_funcs = [CubicSpline(times, position_vector[:,i], extrapolate=True) for i in range(position_vector.shape[1])]

pos_interp = np.column_stack([func(new_times[interval_indices]) for func in pos_interp_funcs])

# Get cartesian versions for the specified time interval
pos_interp_cart = np.apply_along_axis(s_c_vec_conversion, 1, np.copy(pos_interp))


new_times = new_times[interval_indices]
if coord_system == "spherical":
return pos_interp
elif coord_system == "cartesian":
return pos_interp_cart
else:
raise ValueError("Invalid coordinate system specified. Must be 'spherical' or 'cartesian'.")


The term that is squared in the final formula is calculated like this:

def clip_norm_dots(vec1, vec2):
return np.clip(np.einsum('ij,ij->i',vec1, vec2),0,0.9999)
def unit_vector(vector):
return vector / np.linalg.norm(vector)
def get_polarization_loss(receiver, source, separation):

receiver_hat = unit_vector(receiver)
source_hat = unit_vector(source)
separation_hat = unit_vector(separation)
dot_rs = clip_norm_dots(receiver_hat,source_hat)
dot_rsep = clip_norm_dots(receiver_hat,separation_hat)
dot_ssep = clip_norm_dots(source_hat,separation_hat)
denominator = np.sqrt(1-dot_rsep**2)*np.sqrt(1-dot_ssep**2)
denominator = np.maximum(denominator,1e-6)
ploss = (dot_rs - (dot_rsep*dot_ssep))/denominator
return ploss


Where the separation vector is the position vector from interp_position at a given time step. As an MRE, the receiver and source vectors can be set to:

receiver = np.array([1,0,0])
source = np.array([1,1,0])


On my machine, the maximum ploss squared factors calculated with a sample rate of 30 deviated from the one with a sample rate of 60 by a factor of 10. I am looking for the best method to minimize this deviation.

As a note, the interpolated time is also returned, just by a different function for various clarity reasons. The code is identical to the code in interp_position().

Updates @NickODell - here how the functions are called, and additional functions that are used.

traj_arrays = ny.read_traj_data("Traj_Right.txt")
times = ny.get_times(traj_arrays)
positions = np.zeros((len(times),len(Receiver),3))
for i in range(len(Receiver)):
positions[:,i] = ny.get_spherical_position(coords[0], coords[1], traj_arrays)
losses = np.zeros((len(times_interp),len(Receiver),2))
for i in range(len(Receiver)):
losses[:,i,0] = ut.get_polarization_loss(receivers_ew,transmitters_aligned,positions_interp_cartesian[:,i])
losses[:,i,1] = ut.get_polarization_loss(receivers_ns,transmitters_aligned,positions_interp_cartesian[:,i])


Additional relevant functions:

def s_c_vec_conversion(spherical_vec):
x = spherical_vec[0] * np.sin(spherical_vec[1]) * np.cos(spherical_vec[2])
y = spherical_vec[0] * np.sin(spherical_vec[1]) * np.sin(spherical_vec[2])
z = spherical_vec[0] * np.cos(spherical_vec[1])
return np.array([x, y, z])
# cartesian to spherical
def c_s_vec_conversion(cartesian_vec):
r = np.sqrt(cartesian_vec[0]**2+cartesian_vec[1]**2+cartesian_vec[2]**2)
theta = np.arccos(cartesian_vec[2]/r)
phi = np.arctan(cartesian_vec[1]/cartesian_vec[0])
return np.array([r,theta,phi])
def interp_time(times, sample_rate, start_time=0, end_time=None):
if end_time is None:
end_time = times[-1]
new_times = np.arange(times[0], times[-1],1/sample_rate)
interval_indices = np.where((new_times >= start_time) & (new_times <= end_time))[0]
return new_times[interval_indices]

# Translate a point in spherical coordinates relative to a reference point
def translate_point_spherical(r1, theta1, phi1, r2, theta2, phi2):
# Convert reference point and satellite point to Cartesian
p1_cartesian = spherical_to_cartesian(r1, theta1, phi1)
p2_cartesian = spherical_to_cartesian(r2, theta2, phi2)
# Calculate the translated Cartesian coordinates
translated_cartesian = p2_cartesian - p1_cartesian
# Convert back to spherical coordinates
return cartesian_to_spherical(*translated_cartesian)

def get_spherical_position(receiver_lat, receiver_lon, traj_arrays):
pos_vec = np.zeros((len(traj_arrays["Time"]),3))
times = traj_arrays["Time"]
# Process trajectory data
for i, time in enumerate(traj_arrays["Time"]):
latitude = traj_arrays["Latgd"]
longitude = traj_arrays["Long"]
altitude = traj_arrays["Altkm"] * 1000 # Convert altitude to meters
time = traj_arrays["Time"]

# Poker Flat as origin
r, theta, phi = translate_point_spherical(
R_EARTH, np.pi/2 - np.radians(receiver_lat), np.radians(receiver_lon),
R_EARTH + altitude, np.pi/2 - np.radians(latitude), np.radians(longitude)
)
pos_vec[0] = r
pos_vec[1] = (np.degrees(np.arccos(altitude/r)))
phi = np.degrees(phi)
if phi < 0:
phi += 180

pos_vec[2] = (phi)
return pos_vec
def calc_received_power(rocket_pos, gains_rx, gains_tx, ploss):
result_watts = (txPwr * np.multiply(np.multiply(gains_tx,gains_rx),ploss**2))/path_loss(rocket_pos)
result_watts[result_watts<=0]=1e-100
result_watts = result_watts.astype(np.float64)
result_dBm = 10*np.log10(result_watts)+30
return result_dBm

Continue reading...