This page was generated by nbsphinx from docs/examples/Filtering.ipynb.
Interactive online version:

Scipy filters compatibility

Only for Google Colab users:

%pip install --upgrade ipympl speasy
try:
    from google.colab import output

    output.enable_custom_widget_manager()
except:
    print("Not running inside Google Collab")

For all users:

import speasy as spz
from speasy.core.codecs import get_codec
from speasy.signal.filtering import sosfiltfilt
import numpy as np
%matplotlib widget
mms1_products = spz.inventories.tree.cda.MMS.MMS1
# Use this instead if you are not using jupyterlab yet
#%matplotlib notebook
import matplotlib.pyplot as plt
from scipy import signal

Loading data

Let’s use the MMS1 EDP burst data from the CDAWeb to test the filtering.

Here we read the CDF file directly through speasy.core.codecs.get_codec, Speasy’s internal codec layer — the machinery spz.get_data(...) uses behind the scenes to turn files into SpeasyVariables. You normally don’t need it: spz.get_data(...) with a product from the inventory is the public, recommended entry point. We use the codec here only to demonstrate that any CDF file — even one not indexed in an archive — can be loaded into a SpeasyVariable and filtered.

mms1_edp_hmfe_par_epar_brst_l2 = get_codec("cdf").load_variable(variable="mms1_edp_hmfe_par_epar_brst_l2", file="https://cdaweb.gsfc.nasa.gov/pub/data/mms/mms1/edp/brst/l2/hmfe/2021/07/mms1_edp_brst_l2_hmfe_20210701020013_v2.0.0.cdf")

sampling_frequency = 1e9/float(mms1_edp_hmfe_par_epar_brst_l2.time[1]-mms1_edp_hmfe_par_epar_brst_l2.time[0])
sampling_frequency

65535.09404285995

Filtering with scipy.signal

Let’s design a time domain bandpass filter with scipy.signal cutting off frequencies below 100 Hz and above 10 kHz.

sos = signal.iirfilter(N=8, Wn=[100, 10000], rs=80, btype='bandpass',output='sos', fs=sampling_frequency)
w, h = signal.sosfreqz(sos, worN=150000, fs=sampling_frequency)
plt.figure()
plt.subplot(2, 1, 1)
db = 20*np.log10(np.maximum(np.abs(h), 1e-5))
plt.plot(w/np.pi, db)
plt.ylim(-75, 5)
plt.grid(True)
plt.yticks([0, -20, -40, -60, -80])
plt.semilogx()
plt.ylabel('Gain [dB]')
plt.title('Frequency Response')
plt.subplot(2, 1, 2)
plt.plot(w/np.pi, np.angle(h))
plt.grid(True)
plt.yticks([-np.pi, -0.5*np.pi, 0, 0.5*np.pi, np.pi],
           [r'$-\pi$', r'$-\pi/2$', '0', r'$\pi/2$', r'$\pi$'])
plt.ylabel('Phase [rad]')
plt.xlabel('Normalized frequency (1.0 = Nyquist)')
plt.semilogx()
plt.show()

Speasy variables filtering

We can directly filter the data using the sos filter designed with scipy.signal.

filtered_mms1_edp_hmfe_par_epar_brst_l2 = sosfiltfilt(sos, mms1_edp_hmfe_par_epar_brst_l2)

Plotting the results

The oddly precise zoom window below was picked by hand: it is the small burst of wave activity visible around 02:00:46.9 in the full-interval plot above — an interval short enough to show individual oscillations.

start_time = np.datetime64("2021-07-01T02:00:46.9378", "ns")
stop_time = np.datetime64("2021-07-01T02:00:46.9826", "ns")
filtered_mms1_edp_hmfe_par_epar_brst_l2 = filtered_mms1_edp_hmfe_par_epar_brst_l2[start_time:stop_time]
mms1_edp_hmfe_par_epar_brst_l2 = mms1_edp_hmfe_par_epar_brst_l2[start_time:stop_time]

plt.figure()
plt.plot(mms1_edp_hmfe_par_epar_brst_l2.time, mms1_edp_hmfe_par_epar_brst_l2.values, label="Original data")
plt.plot(filtered_mms1_edp_hmfe_par_epar_brst_l2.time, filtered_mms1_edp_hmfe_par_epar_brst_l2.values, label="Filtered data")
plt.legend()
plt.xticks(rotation=45)
plt.tight_layout()