From bbf6e684527eb3222f77b1c24f75d961a10c1935 Mon Sep 17 00:00:00 2001 From: Rob Robinett Date: Sun, 19 Jul 2020 15:00:48 -1000 Subject: [PATCH] Add these Python routines --- c2_noise.py | 37 ++++++++++++++++ noise_plot.py | 118 ++++++++++++++++++++++++++++++++++++++++++++++++++ suntimes.py | 11 +++++ wav_window.py | 41 ++++++++++++++++++ 4 files changed, 207 insertions(+) create mode 100755 c2_noise.py create mode 100755 noise_plot.py create mode 100755 suntimes.py create mode 100755 wav_window.py diff --git a/c2_noise.py b/c2_noise.py new file mode 100755 index 0000000..daeb280 --- /dev/null +++ b/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)) diff --git a/noise_plot.py b/noise_plot.py new file mode 100755 index 0000000..1d6c021 --- /dev/null +++ b/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) diff --git a/suntimes.py b/suntimes.py new file mode 100755 index 0000000..90780e7 --- /dev/null +++ b/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] ) diff --git a/wav_window.py b/wav_window.py new file mode 100755 index 0000000..af8b951 --- /dev/null +++ b/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)