Generate synthetic HR GNSS noise

The basic outline for this process is in Melgar et al. (2020, JGR) and it is based on Boore (1983, BSSA) and Graves & Pitarka (2010, BSSA). For more details consult those references. Here I will simply outline with code snippets the process using built in MudPy functions.

First import the relevant mudpy modules (and matplotlib for plotting the results at the end)

from mudpy.hfsims import windowed_gaussian,apply_spectrum
from mudpy.forward import gnss_psd
import matplotlib.pyplot as plt

windowed_gaussain simply generates a white noise time series. gnss_psd reads the different reference power spectra from file, again, see Melgar et al. (2020) for the different noise spectra available. Finally apply_spectrum applies the noise spectra, in the frequency domain, tot he white noise time series.

First some parameters:

#define sample rate of the time series (in seconds)
dt=1.0

#noise model (this is the percentile, here we start with median noise)
percentile=50  #this can be 1,10,20,30,40,50,60,70,80, or 90

#duration of time series (in seconds
duration=2400.

Now build a random Gaussian white noise time series

# get white noise
std=1.0 # this is a dummy parameter give it any value
E_noise=windowed_gaussian(duration,dt,std=std,window_type=None)
N_noise=windowed_gaussian(duration,dt,std=std,window_type=None)
Z_noise=windowed_gaussian(duration,dt,std=std,window_type=None)
noise=windowed_gaussian(duration,dt,std=std,window_type=None)

Now get the power spectra for each component of motion of noise at the specified percentile level

f,Epsd,Npsd,Zpsd=gnss_psd(level=percentile,return_as_frequencies=True,return_as_db=False)

Note we return the PSD in physical units (as opposed to dB as in the paper). This is important for correctly applying the spectra. For plotting, to compare against the results in the paper, you can return the spectra in dB. Also setting return_as_frequencies to False will return the periods at which the PSD is defined, as opposed to the frequencies.

Now convert the power spectra to amplitude spectra, and apply tot he white noise time series.

#Covnert PSDs to amplitude spectrum
Epsd = Epsd**0.5
Npsd = Npsd**0.5
Zpsd = Zpsd**0.5

#apply the spectrum
E_noise=apply_spectrum(E_noise,Epsd,f,dt,is_gnss=True)
N_noise=apply_spectrum(N_noise,Npsd,f,dt,is_gnss=True)
Z_noise=apply_spectrum(Z_noise,Zpsd,f,dt,is_gnss=True)

And done, here's a quick plot of the time series:

plt.figure()
plt.plot(E_noise)
plt.plot(N_noise+0.2)
plt.plot(Z_noise+0.4)
plt.xlabel('Seconds')
plt.ylabel('Position')
plt.legend(['East','North','Up'],loc=0)
plt.ylabel('Position (m)')
plt.title('%dth percentile synthetic GNSS noise\n (waveforms offset for clarity)' % percentile)

Compare to the 10th percentile (less noise) result, by changing to percentile=10:

Happy noising!