<jupyter_start><jupyter_text>Analysis of lightning discharges as seen on a Software Defined Radio  Experimental Setup:RTL-SDR connected to a Linux laptop running Ubuntu. Radio tuned to 30MHz and ~4s of time series (I,Q) data collected at 2Ms/s. Radio is connected to a "rabbit ears" dipole antenna inside my home. 30MHz was choses as it is fairly quiet and as low frequency as the device will go without additional hardware. Noise is around .3 units so only data containing peaks about 0.5 units are saved using a numpy save file. These files are saved and timestamped for future analysis. This notebook shows what we can see during a lightning storm on a very simple (~$35) SDR set up. There is a lot we can do to improve this! Our goal is to equip Sage nodes with SDRs in the future to test these as a affordable ubiquitios lightbing detection network. Future work includes using am upconvertor to access lower frequencies.<jupyter_code>#Import the goodness
import numpy as np
from matplotlib import pyplot as plt
from scipy import signal
%matplotlib inline
#load a time series file. Collected at 1626 UTC on the 19th of July
ts = np.load('/users/scollis/data/sage_lightning/data/event_200719_1626_00_0p6896.npy', allow_pickle=True).item()
#Create a time array using the known sampling frequency
dt = 1./2.048e6
xtime = np.arange(len(ts['sig'])) * dt

#Quick and dirty look at the data, note we use np.abs as the signal is complex

myfi = plt.figure(figsize=[15,5])
plt.plot(xtime, np.abs(ts['sig']))
plt.xlabel('Signal (Arb)')
plt.xlabel('Time since start (s)')
#Quick code to find the element where the maximum of the signal is

maxme = np.where(np.abs(ts['sig']) == np.max(np.abs(ts['sig'])))[0][0]
t_max = xtime[maxme]
print(t_max)
#What time window do we want to look at?
window = 4e-3
#Zoom in to the max signal times
myfi = plt.figure(figsize=[15,5])
plt.plot(xtime, np.abs(ts['sig']))
plt.ylabel('Signal (Arb)')
plt.xlabel('Time since start (s)')
t0 = t_max - window/2.
t1 =  t_max + window/2.
plt.xlim([t0,t1])
#Code Courtesy of <NAME> TTU
#Plot the time series and rolling FFT of the signal

fs = 2.048e6
sub = slice(int(t0*fs),int(t1*fs))

dt = 1./fs
xtime = np.arange(len(ts['sig'])) * dt
myfi, axs = plt.subplots(2,1, figsize=[15,10], sharex=True)
axs[0].plot(xtime[sub]-t0, np.abs(ts['sig'])[sub])

_,_,_,sgimg = axs[1].specgram(ts['sig'][sub], Fs=fs, vmin=-100, vmax=-70)
axs[0].set_ylabel('Signal (Arb)')
axs[1].set_ylabel('Spectrum (Hz)')
axs[1].set_xlabel('Time since start (s)')
plt.colorbar(sgimg, orientation='horizontal', ax=axs[1])<jupyter_output><empty_output><jupyter_text>You can see some coherent signals from radio etc... But the lightning impulse is broad band filling the whole bandpass window.<jupyter_code>#The next but of code is designed to look for multiple peaks in the data. 


#First, Smooth the data using a Savgol filter.. otherwise we get too many peaks
filtered_data = signal.savgol_filter(np.abs(ts['sig']),5 ,2)

#Now use SciPy's peak finder to find peaks.. These myst be at least .25 units above noise and 
# at least 10,000 elements apart

locs, props = signal.find_peaks(filtered_data, height = .25, distance = 10000, width = 1)
locs<jupyter_output><empty_output><jupyter_text>So we have three distinct peaks in this time series.<jupyter_code>#Lets make Eric's code into a function

def plot_pulse(ts, time, window):
    t0 = time - window/2.
    t1 =  time + window/2.
    fs = 2.048e6
    sub = slice(int(t0*fs),int(t1*fs))
    dt = 1./fs
    xtime = np.arange(len(ts['sig'])) * dt
    myfi, axs = plt.subplots(2,1, figsize=[15,10], sharex=True)
    axs[0].plot(xtime[sub]-t0, np.abs(ts['sig'])[sub])
    _,_,_,sgimg = axs[1].specgram(ts['sig'][sub], Fs=fs, vmin=-100, vmax=-70)
    axs[0].set_ylabel('Signal (Arb)')
    axs[1].set_ylabel('Spectrum (Hz)')
    axs[1].set_xlabel('Time since start (s)')
    plt.colorbar(sgimg, orientation='horizontal', ax=axs[1])
#Now lets loop over all peaks and plot them

for loc in locs:
    plot_pulse(ts, xtime[loc], 1e-3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12790697674418605, 'fraction_numerical': 0.03240969816922316, 'mean_word_length': 4.775714285714286, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '267027', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1519, 'n_tokens_neox': 1421, 'n_words': 489, 'path': 'lightning/analyze_lightning_signal.ipynb', 'repo_name': 'jcorner1/SageScienceExamples', 'size': '4048'}