Add these Python routines

c2_noise.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Filename: c2_noise.py
+# Program to extract the noise level from the 'wsprd -c' C2 format file
+## V1 by Christoph Mayer. This version V1.1 by Gwyn Griffiths to output a single value
+## being the total power (dB arbitary scale) in the lowest 30% of the Fourier coefficients
+## between 1369.5 and 1630.5 Hz where the passband is flat.
+
+import struct
+import sys
+import numpy as np
+
+fn = sys.argv[1] ## '000000_0001.c2'
+
+with open(fn, 'rb') as fp:
+     ## decode the header:
+     filename,wspr_type,wspr_freq = struct.unpack('<14sid', fp.read(14+4+8))
+
+     ## extract I/Q samples
+     samples = np.fromfile(fp, dtype=np.float32)
+     z = samples[0::2]+1j*samples[1::2]
+     #print(filename,wspr_type,wspr_freq,samples[:100], len(samples), z[:10])
+
+     ## z contains 45000 I/Q samples
+     ## we perform 180 FFTs, each 250 samples long
+     a     = z.reshape(180,250)
+     a    *= np.hanning(250)
+     freqs = np.arange(-125,125, dtype=np.float32)/250*375 ## was just np.abs, square to get power
+     w     = np.square(np.abs(np.fft.fftshift(np.fft.fft(a, axis=1), axes=1)))
+     ## these expressions first trim the frequency range to 1369.5 to 1630.5 Hz to ensure
+     ## a flat passband without bias from the shoulders of the bandpass filter
+     ## i.e. array indices 38:213
+     w_bandpass=w[0:179,38:213]
+     ## partitioning is done on the flattened array of coefficients
+     w_flat_sorted=np.partition(w_bandpass, 9345, axis=None)
+     noise_level_flat=10*np.log10(np.sum(w_flat_sorted[0:9344]))
+     print(' %6.2f' % (noise_level_flat))

noise_plot.py

@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Filename: noise_plot.py
+# April-May  2019  Gwyn Griffiths G3ZIL
+# Use matplotlib to plot noise levels recorded by wsprdaemon by the sox stats RMS and sox stat -freq methods
+# V0 Testing prototype 
+
+# Import the required Python modules and methods some may need downloading 
+from __future__ import print_function
+import math
+import datetime
+#import scipy
+import numpy as np
+from numpy import genfromtxt
+import csv
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+#from matplotlib import cm
+import matplotlib.dates as mdates
+import sys
+
+# Get cmd line args
+reporter=sys.argv[1]
+maidenhead=sys.argv[2]
+output_png_filepath=sys.argv[3]
+calibration_file_path=sys.argv[4]
+csv_file_path_list=sys.argv[5].split()    ## noise_plot.py KPH "/home/pi/.../2200 /home/pi/.../630 ..."
+
+# read in the reporter-specific calibration file and print out
+# if one didn't exist the bash script would have created one
+# the user can of course manually edit the specific noise_cal_vals.csv file if need be
+cal_vals=genfromtxt(calibration_file_path, delimiter=',')
+nom_bw=cal_vals[0]
+ne_bw=cal_vals[1]
+rms_offset=cal_vals[2]
+freq_offset=cal_vals[3]
+fft_band=cal_vals[4]
+threshold=cal_vals[5]
+
+# need to set the noise equiv bw for the -freq method. It is 322 Hz if nom bw is 500Hz else it is ne_bw as set
+if nom_bw==500:
+    freq_ne_bw=322
+else:
+    freq_ne_bw=ne_bw
+
+x_pixel=40
+y_pixel=30
+my_dpi=50         # set dpi and size for plot - these values are largest I can get on Pi window, resolution is good
+fig = plt.figure(figsize=(x_pixel, y_pixel), dpi=my_dpi)
+fig.subplots_adjust(hspace=0.4, wspace=0.4)
+plt.rcParams.update({'font.size': 18})
+
+# get, then set, start and stop time in UTC for use in overall title of charts
+stop_t=datetime.datetime.utcnow()
+start_t=stop_t-datetime.timedelta(days=1)   ### Plot last 24 hours
+stop_time=stop_t.strftime('%Y-%m-%d %H:%M')
+start_time=start_t.strftime('%Y-%m-%d %H:%M')
+
+fig.suptitle("Site: '%s' Maidenhead: '%s'\n Calibrated noise (dBm in 1Hz, Temperature in K) red=RMS blue=FFT\n24 hour time span from '%s' to '%s' UTC" % (reporter, maidenhead, start_time, stop_time), x=0.5, y=0.99, fontsize=24)
+
+# Process the list of csv  noise files
+j=1
+# get number of csv files to plot then divide by three and round up to get number of rows
+plot_rows=int(math.ceil((len(csv_file_path_list)/3.0)))
+for csv_file_path in csv_file_path_list:
+    # matplotlib x axes with time not straightforward, get timestamp in separate 1D array as string
+    timestamp  = genfromtxt(csv_file_path, delimiter=',', usecols=0, dtype=str)
+    noise_vals = genfromtxt(csv_file_path, delimiter=',')[:,1:]  
+
+    n_recs=int((noise_vals.size)/15)              # there are 15 comma separated fields in each row, all in one dimensional array as read
+    noise_vals=noise_vals.reshape(n_recs,15)      # reshape to 2D array with n_recs rows and 15 columns
+
+    # now  extract the freq method data and calibrate
+    freq_noise_vals=noise_vals[:,13]  ### +freq_offset+10*np.log10(1/freq_ne_bw)+fft_band+threshold
+    rms_trough_start=noise_vals[:,3]
+    rms_trough_end=noise_vals[:,11]
+    rms_noise_vals=np.minimum(rms_trough_start, rms_trough_end)
+    rms_noise_vals=rms_noise_vals     #### +rms_offset+10*np.log10(1/ne_bw)
+    ov_vals=noise_vals[:,14]          ### The OV (overload counts) reported by Kiwis have been added in V2.9
+
+    # generate x axis with time
+    fmt = mdates.DateFormatter('%H')          # fmt line sets the format that will be printed on the x axis
+    timeArray = [datetime.datetime.strptime(k, '%d/%m/%y %H:%M') for k in timestamp]     # here we extract the fields from our original .csv timestamp
+
+    ax1 = fig.add_subplot(plot_rows, 3, j)
+    ax1.plot(timeArray, freq_noise_vals, 'b.', ms=2)
+    ax1.plot(timeArray, rms_noise_vals, 'r.', ms=2)
+    # ax1.plot(timeArray, ov_vals, 'g.', ms=2)       # OV values will need to be scaled if they are to appear on the graph along with noise levels
+
+    ax1.xaxis.set_major_formatter(fmt)
+
+    path_elements=csv_file_path.split('/')
+    plt.title("Receiver %s   Band:%s" % (path_elements[len(path_elements)-3], path_elements[len(path_elements)-2]), fontsize=24)
+
+    #axes = plt.gca()
+    # GG chart start and stop UTC time as end now and start 1 day earlier, same time as the x axis limits
+    ax1.set_xlim([datetime.datetime.utcnow()-datetime.timedelta(days=1), datetime.datetime.utcnow()])
+    # first get 'loc' for the hour tick marks at an interval of 2 hours then use 'loc' to set the major tick marks and grid
+    loc=mpl.dates.HourLocator(byhour=None, interval=2, tz=None)
+    ax1.xaxis.set_major_locator(loc)
+
+    #   set y axes lower and upper limits
+    y_dB_lo=-175
+    y_dB_hi=-105
+    y_K_lo=10**((y_dB_lo-30)/10.)*1e23/1.38
+    y_K_hi=10**((y_dB_hi-30)/10.)*1e23/1.38
+    ax1.set_ylim([y_dB_lo, y_dB_hi])
+    ax1.grid()
+
+    # set up secondary y axis
+    ax2 = ax1.twinx()
+    # automatically set its limits to be equivalent to the dBm limits
+    ax2.set_ylim([y_K_lo, y_K_hi])
+    ax2.set_yscale("log")
+
+    j=j+1  
+fig.savefig(output_png_filepath)

suntimes.py

@@ -0,0 +1,11 @@
+import datetime, sys
+from astral import Astral, Location
+from datetime import date
+lat=float(sys.argv[1])
+lon=float(sys.argv[2])
+zone=sys.argv[3]
+l = Location(('wsprep', 'local', lat, lon, zone, 0))
+l.sun()
+d = date.today()
+sun = l.sun(local=True, date=d)
+print( str(sun['sunrise'])[11:16] + " " + str(sun['sunset'])[11:16] )

wav_window.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Filename: wav_window_v1.py
+# January  2020  Gwyn Griffiths
+# Program to apply a Hann window to a wsprdaemon wav file for subsequent processing by sox stat -freq (initially at least)
+
+from __future__ import print_function
+import math
+import scipy
+import scipy.io.wavfile as wavfile
+import numpy as np
+import wave
+import sys
+
+WAV_INPUT_FILENAME=sys.argv[1]
+WAV_OUTPUT_FILENAME=sys.argv[2]
+
+# Set up the audio file parameters for windowing
+# fs_rate is passed to the output file
+fs_rate, signal = wavfile.read(WAV_INPUT_FILENAME)   # returns sample rate as int and data as numpy array
+# set some constants
+N_FFT=352                                   # this being the number expected
+N_FFT_POINTS=4096                           # number of input samples in each sox stat -freq FFT (fixed)
+                                            # so N_FFT * N_FFT_POINTS = 1441792 samples, which at 12000 samples per second is 120.15 seconds
+                                            # while we have only 120 seconds, so for now operate with N_FFT-1 to have all filled
+                                            # may decide all 352 are overkill anyway
+N=N_FFT*N_FFT_POINTS
+w=np.zeros(N_FFT_POINTS)
+
+output=np.zeros(N, dtype=np.int16)          # declaring as dtype=np.int16 is critical as the wav file needs to be 16 bit integers
+
+# create a N_FFT_POINTS array with the Hann weighting function
+for i in range (0, N_FFT_POINTS):
+  x=(math.pi*float(i))/float(N_FFT_POINTS)
+  w[i]=np.sin(x)**2
+
+for j in range (0, N_FFT-1):
+  offset=j*N_FFT_POINTS
+  for i in range (0, N_FFT_POINTS):
+     output[i+offset]=int(w[i]*signal[i+offset])
+wavfile.write(WAV_OUTPUT_FILENAME, fs_rate, output)