Imagine a beam of light racing around Earth’s equator at 299,702,547 meters per second. In just 0.134 seconds, it would complete the journey and hit us in the back. This calculation is straightforward:
\(T = \dfrac{40075000}{299702547} = 0.13371591399\) seconds
But here’s where things get interesting: while light travels in straight lines and disappears beyond the horizon due to Earth’s curvature, radio waves in the High Frequency (HF) band (3-30 MHz) have a special trick up their sleeve. They can bounce off the ionosphere—a layer of charged particles in Earth’s upper atmosphere—allowing them to travel far beyond the horizon.
This remarkable property of HF radio waves enables long-distance communication without the need for wires or satellites. But it gets even more fascinating: because these waves follow predictable paths through the atmosphere, we can actually calculate how far away a radio signal originated. All we need is a reference time to compare the signal against, and we can determine the source’s distance with surprising accuracy.
This ability to “bend” around Earth’s curvature makes HF radio waves not just a means of communication, but also a tool for measuring distances across the globe.
import numpy as np
import soundfile as sf
import matplotlib.pyplot as plt
from scipy.io import wavfile
from wavfile import read, write
from ipywidgets import widgets
import plotly.graph_objects as go
import plotly.io as pio
import plotly.express as px
import pandas as pd
import matplotlib.patches as patches
import pymc as pm
import arviz as azstuff = read('wav_files/20240325T112459Z_15000000_iq.wav', readmarkers=True,
readmarkerlabels=True, readmarkerslist=True, readpitch=True, readloops=True)# 5 msec of samples at 12000 samples per second is 60 samples
d = 60
# calculate the number of 5 msec chunks
nc = stuff[1].shape[0] // d
# put into a 2D-matrix of shape (nc, d)
chan_I = stuff[1][:d*nc,0].reshape((nc, d))
chan_Q = stuff[1][:d*nc,1].reshape((nc, d))
# do Fourier transform of complex signal
fft = np.abs(np.fft.fft(chan_I +1.0j*chan_Q))
fft = np.fft.fftshift(fft)f, ax = plt.subplots(nrows=1, ncols=3, figsize=(15,5), dpi=300)
ax[0].imshow(fft[0:60000,22:39], norm='log', aspect='auto')
ax[0].set_xticks(np.arange(17))
ax[0].set_xticklabels((np.arange(17)-8)*200, rotation=90)
ax[0].set_yticks(np.linspace(0,60000,21))
ax[0].set_yticklabels([f'{a:.2f}' for a in (np.linspace(0,60000,21)*0.005/60 - 15/60)])
ax[0].set_ylabel('Minutes')
ax[0].set_xlabel('Frequency from Carrier (Hz)')
ax[2].imshow(fft[3052:3235,34:38], norm='log', aspect='auto')
ax[2].set_xticks(np.arange(4))
ax[2].set_xticklabels((np.arange(4)+4)*200, rotation=90)
ax[1].imshow(fft[3052:3235,23:27], norm='log', aspect='auto')
ax[1].set_xticks(np.arange(4))
ax[1].set_xticklabels((np.arange(4)[::-1]+4)*-200, rotation=90)
rect1 = patches.Rectangle((.5, 2400), 4, 1500, linewidth=1.5, edgecolor='k', facecolor='none', ls="-")
rect2 = patches.Rectangle((11.5, 2400), 4, 1500, linewidth=1.5, edgecolor='k', facecolor='none')
ax[0].add_patch(rect1)
ax[0].add_patch(rect2)
plt.show()
## Plot All Tracesoffset = 3039
endset = 200
repeat = 12000
df = pd.DataFrame()
for i in range(40):
df['wwv1_'+str(i)] = fft[i*repeat+offset:i*repeat+offset+endset, 25]
df['wwv2_'+str(i)] = fft[i*repeat+offset:i*repeat+offset+endset, 35]
df['wwvh1_'+str(i)] = fft[i*repeat+offset:i*repeat+offset+endset, 24]
df['wwvh2_'+str(i)] = fft[i*repeat+offset:i*repeat+offset+endset, 36]
df = df.copy()
df_max = df.max().max()
# f, ax = plt.subplots(ncols=5, nrows=2, figsize=(12.5, 4), dpi=300)
# for i in range(10):
# ax[i//5][i%5].plot(df.iloc[:,4*i]/df_max, c=(1, 0, 0))
# ax[i//5][i%5].plot(df.iloc[:,4*i+1]/df_max, c=(0, 0, 1))
# ax[i//5][i%5].plot(df.iloc[:,4*i+2]/df_max, c=(0, 180/255, 0))
# ax[i//5][i%5].plot(df.iloc[:,4*i+3]/df_max, c=(1, 200/255, 0))
# ax[i//5][i%5].set_title(f'Minute {i+1}')
# ax[i//5][i%5].set_yticks([])
# plt.tight_layout()offset = 3039
endset = 200
repeat = 12000
fours = np.linspace(0, 36, 10)
fig = px.line(data_frame=df, x=np.arange(0,1000,5), y=[df.iloc[:,fours+0].mean(axis=1),
df.iloc[:,fours+1].mean(axis=1),
df.iloc[:,fours+2].mean(axis=1),
df.iloc[:,fours+3].mean(axis=1)],
color_discrete_sequence=["rgb(255, 0, 0)", "rgb(0, 0, 255)", "rgb(0, 180, 0)", "rgb(255, 200, 0)"],
)
series_names = ["WWV Lower", "WWV Upper", "WWVH Lower", "WWVH Upper"]
for idx, name in enumerate(series_names):
fig.data[idx].name = name
fig.data[idx].hovertemplate = name
fig.update_layout(legend={"title":"Trace"})
fig.update_layout(width=1200, height=500,
updatemenus= [
{'buttons': [
{
'method': 'restyle',
'label': '0-9',
'args': [{'y': [df.iloc[:,fours+0].mean(axis=1),
df.iloc[:,fours+1].mean(axis=1),
df.iloc[:,fours+2].mean(axis=1),
df.iloc[:,fours+3].mean(axis=1)]}]
},
{
'method': 'restyle',
'label': '10-19',
'args': [{'y': [df.iloc[:,fours+40].mean(axis=1),
df.iloc[:,fours+41].mean(axis=1),
df.iloc[:,fours+42].mean(axis=1),
df.iloc[:,fours+43].mean(axis=1)]}]
},
{
'method': 'restyle',
'label': '20-29',
'args': [{'y': [df.iloc[:,fours+80].mean(axis=1),
df.iloc[:,fours+81].mean(axis=1),
df.iloc[:,fours+82].mean(axis=1),
df.iloc[:,fours+83].mean(axis=1)]}]
},
{
'method': 'restyle',
'label': '30-39',
'args': [{'y': [df.iloc[:,fours+120].mean(axis=1),
df.iloc[:,fours+121].mean(axis=1),
df.iloc[:,fours+122].mean(axis=1),
df.iloc[:,fours+123].mean(axis=1)]}]
}
],
'direction': 'down',
'showactive': True,
}
] )
fig.update_traces(line=dict(width=4))
fig.update_layout(
plot_bgcolor="White",
xaxis = dict(
tickmode = 'linear',
tick0 = 0,
dtick = 100,
#tickvals = np.arange(-601,1001,100),
#ticktext = [f"{i}" for i in np.arange(-601,1001,100)]
),
)
fig.update_xaxes(
mirror=True,
ticks='outside',
showline=True,
linecolor='black',
gridcolor='lightgrey',
#zerolinecolor='darkgrey',
)
fig.update_yaxes(
mirror=True,
ticks='outside',
showline=True,
linecolor='black',
gridcolor='lightgrey',
zerolinecolor='black',
)
fig.add_vline(x=85, line_width=4, line_dash="dash", line_color="black")
# Save figure as standalone HTML
fig.write_html("plotly_figure.html")
# For embedding in Quarto, you can use:
from IPython.display import HTML
HTML(filename="plotly_figure.html")
# Save to HTML string
html_str = fig.to_html(include_plotlyjs='cdn')
# Save to file for embedding
with open('figure.html', 'w') as f:
f.write(html_str)m = 100 # will be fitted
d =75 # gonna be fixed
w1 = 1 # probably gonna be fixed
w2 = .25 # probably gonna be fixed
A = 1 # will be fitted
x = np.arange(0, 200)
s1 = (1/(1+np.exp(w1*(-x+(m-d)))))
s2 = (1/(1+np.exp(-w2*(-x+(m+d)))))
s = (s1 + s2)-1
#plt.plot(s1)
#plt.plot(s2)
# f, ax = plt.subplots(nrows=1, ncols=3, figsize=(12,4), dpi=300)
# ax[0].plot(s1, label=r'$\dfrac{1}{1+e^{w(m-d-x)}}$', c='k', lw=1.5)
# ax[1].plot(s2, label=r'$\dfrac{1}{1+e^{-w(m-d-x)}}$', c='k', lw=1.5)
# ax[2].plot(s, label=r'$\dfrac{1}{1+e^{w(m-d-x)}} + \dfrac{1}{1+e^{-w(m-d-x)}}-1$', c='k', lw=1.5)
# ax[0].legend(loc=4, fontsize=8)
# ax[1].legend(loc=3, fontsize=8)
# ax[2].legend(loc=8, fontsize=8)
fig = plt.figure(figsize=(8, 8), layout="constrained", dpi=300)
spec = fig.add_gridspec(3, 2)
#annotate_axes(ax0, 'ax0')
ax00 = fig.add_subplot(spec[0, 0])
ax00.plot(s1, label=r'$\dfrac{1}{1+e^{w_f(m-d-x)}}$', c='k', lw=1.5)
ax00.legend(loc=4, fontsize=16)
ax01 = fig.add_subplot(spec[0, 1])
ax01.plot(s2, label=r'$\dfrac{1}{1+e^{-w_b(m-d-x)}}$', c='k', lw=1.5)
ax01.legend(loc=3, fontsize=16)
ax11 = fig.add_subplot(spec[1, :])
ax11.plot(s, label=r'$\dfrac{1}{1+e^{w_f(m-d-x)}} + \dfrac{1}{1+e^{-w_b(m-d-x)}}-1$', c='k', lw=1.5)
ax11.legend(loc=8, fontsize=16)
x = np.arange(0, 200)
shift = -10
ns1 = (1/(1+np.exp(w1*(-x-shift+(m-d)))))
ns2 = (1/(1+np.exp(-w2*(-x-shift+(m+d)))))
ns = (ns1 + ns2)-1
ax22 = fig.add_subplot(spec[2, :])
ax22.plot(np.arange(len(ns)), ns*1.6 + np.random.randn(len(ns))*.1000 + .2000, c='r', lw=1.5)
ax22.plot(np.arange(len(ns)), s*1 , c='k', lw=1.5, ls="-")
#ax22.legend(loc=8, fontsize=16)
ax11.vlines(100, 0.95, 1.05, color='r', ls="-", lw=1.5)
ax11.vlines(175, 0.0, 0.6, color='r', ls="-", lw=1.5)
ax11.vlines(25, 0.0, 0.6, color='r', ls="-", lw=1.5)
def burst_model(t, A, m, off, shift, w1, w2, d, σ):
ns1 = (1/(1+np.exp(w1*(-t-shift+(m-d)))))
ns2 = (1/(1+np.exp(-w2*(-t-shift+(m+d)))))
return off + A*(ns1 + ns2 -1) + np.random.randn(len(t))*σ
t = np.arange(0, 200)
m = 100
shift = -5
σ = 0.5
A = 10
off = -0.5
w1 = 0.5
w2 = 1
d = 80
burst = burst_model(t, A, m, off, shift, w1, w2, d, σ)
#%matplotlib widget
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,4))
ax.plot(t, burst, c='r')
plt.show()
x = np.arange(200)
wwv = np.zeros((8,200))
wwvh = np.zeros((8,200))
for l, i in enumerate(np.linspace(0, 120, 4)):
for j in range(2):
wwv[int(l)*2+j,:] = df.iloc[:,fours+j+i].mean(axis=1).to_numpy()
wwvh[int(l)*2+j,:] = df.iloc[:,fours+j+i+2].mean(axis=1).to_numpy()
# f, ax = plt.subplots(nrows=2, ncols=1, figsize=(8,16))
# ax[0].plot(wwv.T)
# ax[1].plot(wwvh.T)
# plt.show()
wwv.shape(8, 200)with pm.Model() as model:
# Define priors
wwv_σ = pm.Exponential("wwv_σ", 1/3000, shape=(wwv.shape[0], 1))
wwv_m = pm.Normal("wwv_m", mu=99.6, sigma=1, shape=(wwv.shape[0], 1))
d = 80
#w = 0.18
off_wwv = pm.Exponential("off_wwv", 1/5000, shape=(wwv.shape[0], 1))
wwv_A = pm.Exponential("wwv_A", 1/100000, shape=(wwv.shape[0], 1))
wwvh_σ = pm.Exponential("wwvh_σ", 1/3000, shape=(wwv.shape[0], 1))
wwvh_m = pm.Normal("wwvh_m", mu=103.6, sigma=1, shape=(wwv.shape[0], 1))
d = 80
#w = 0.18
#w_f = pm.Uniform("w_f", lower=0.010, upper=0.50)
#w_b = pm.Uniform("w_b", lower=0.010, upper=0.50)
w_f = pm.Normal("w_f", mu=0.13, sigma=0.025, shape=(wwv.shape[0], 1))
w_b = pm.Normal("w_b", mu=0.26, sigma=0.025, shape=(wwv.shape[0], 1))
w_f_wwvh = pm.Normal("w_f_wwvh", mu=0.13, sigma=0.025, shape=(wwv.shape[0], 1))
w_b_wwvh = pm.Normal("w_b_wwvh", mu=0.26, sigma=0.025, shape=(wwv.shape[0], 1))
off_wwvh = pm.Exponential("off_wwvh", 1/5000, shape=(wwv.shape[0], 1))
wwvh_A = pm.Exponential("wwvh_A", 1/100000, shape=(wwv.shape[0], 1))
# model
wwv_s = wwv_A * (((1/(1+np.exp((-x+(wwv_m-d))/w_f)))+(1/(1+np.exp(-1*(-x+(wwv_m+d))/w_b))))-1) + off_wwv
wwvh_s = wwvh_A * (((1/(1+np.exp((-x+(wwvh_m-d))/w_f_wwvh)))+(1/(1+np.exp(-1*(-x+(wwvh_m+d))/w_b_wwvh))))-1) + off_wwvh
# Define likelihood
wwv_lh = pm.Normal("wwv_lh", mu=wwv_s, sigma=wwv_σ, observed=wwv)
wwvh_lh = pm.Normal("wwvh_lh", mu=wwvh_s, sigma=wwvh_σ, observed=wwvh)
wwv_m_ms = pm.Deterministic('wwv_m_ms', 5*wwv_m)
wwvh_m_ms = pm.Deterministic('wwvh_m_ms', 5*wwvh_m)
delay_ms = pm.Deterministic('delay_ms', wwvh_m_ms-wwv_m_ms)
# Inference!
model = pm.sample( tune=4000, draws=4000, target_accept=0.95, idata_kwargs = {'log_likelihood': True}, nuts_sampler="pymc")
model.extend(pm.sample_posterior_predictive(model))Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [wwv_σ, wwv_m, off_wwv, wwv_A, wwvh_σ, wwvh_m, w_f, w_b, w_f_wwvh, w_b_wwvh, off_wwvh, wwvh_A]Sampling 4 chains for 4_000 tune and 4_000 draw iterations (16_000 + 16_000 draws total) took 52 seconds.
Sampling: [wwv_lh, wwvh_lh]#%matplotlib widget
f, ax = plt.subplots(nrows=15, ncols=2, figsize=(10,30))#, dpi=300)
az.plot_trace(model, axes=ax)#, var_names=['wwv_m_ms', 'wwvh_m_ms', 'delay_ms', 'w_f', 'w_b'], axes=ax)#, 'A_0', 'Rxy'])
plt.tight_layout()
plt.show()
#%matplotlib widget
f, ax = plt.subplots(nrows=5, ncols=2, figsize=(10,10), dpi=300)
az.plot_trace(model, var_names=['wwv_m_ms', 'wwvh_m_ms', 'delay_ms', 'w_f', 'w_b'], axes=ax)#, 'A_0', 'Rxy'])
plt.tight_layout()
plt.show()
means = model.posterior.delay_ms.mean(axis=(0, 1, 3)).to_numpy()
stds = model.posterior.delay_ms.std(axis=(0, 1, 3)).to_numpy()print(means)
print(stds)[18.66503136 19.01079225 19.67016744 19.33521841 20.01496352 19.60958308
20.97870987 19.88572992]
[0.79124762 0.54791341 0.79155523 0.4705195 0.34865943 0.53787142
0.70727115 0.63355319]w_f_mean = model.posterior.w_f.mean().to_numpy()
w_f_std = model.posterior.w_f.std().to_numpy()
w_b_mean = model.posterior.w_b.mean().to_numpy()
w_b_std = model.posterior.w_b.std().to_numpy()w_f_meanarray(0.13082639)w_f_stdarray(0.02443936)wwvh_m_means = model.posterior.wwvh_m_ms.mean(axis=(0, 1, 3)).to_numpy()
wwv_m_means = model.posterior.wwv_m_ms.mean(axis=(0, 1, 3)).to_numpy()wwv_m_meansarray([498.69768125, 498.82687994, 498.39402808, 499.02650424,
499.18136159, 499.19634433, 498.36120488, 498.68941852])wwvh_m_meansarray([517.36271261, 517.83767218, 518.06419551, 518.36172265,
519.19632512, 518.80592742, 519.33991474, 518.57514844])foursarray([ 0., 4., 8., 12., 16., 20., 24., 28., 32., 36.])wwv1 = df.iloc[:,fours+0].mean(axis=1).to_numpy()
wwv2 = df.iloc[:,fours+1].mean(axis=1).to_numpy()
wwvh1 = df.iloc[:,fours+2].mean(axis=1).to_numpy()
wwvh2 = df.iloc[:,fours+3].mean(axis=1).to_numpy()type(wwv)numpy.ndarraywwv_pp = model.posterior_predictive.wwv_lh.to_numpy()
wwvh_pp = model.posterior_predictive.wwvh_lh.to_numpy()#%matplotlib widget
f, ax = plt.subplots(nrows=8, ncols=2, figsize=(16,24))
for i in range(8):
for j in range(10):
ax[i][0].scatter(x, wwv_pp[0, j, i, :], c='k', s=10, alpha=0.5, edgecolors=None, linewidths=0)
ax[i][1].scatter(x, wwvh_pp[0, j, i, :], c='k', s=10, alpha=0.5, edgecolors=None, linewidths=0)
ax[i][0].plot(x, wwv[i, :], c='r', lw=2)
ax[i][1].plot(x, wwvh[i, :], c='r', lw=2)
ax[i][0].set_title(f'WWV {i}')
ax[i][1].set_title(f'WWVH {i}')
ax[i][0].set_xlabel('Time (ms)')
ax[i][1].set_xlabel('Time (ms)')
ax[i][0].set_ylabel('Amplitude')
ax[i][1].set_ylabel('Amplitude')
ax[i][0].set_xlim(0, 200)
ax[i][1].set_xlim(0, 200)
ax[i][0].set_ylim(-25000, 200000)
ax[i][1].set_ylim(-25000, 200000)
plt.tight_layout()
plt.show()