jaseg/master-thesis
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Repo re-org

Browse source
This commit is contained in:
jaseg 2021-04-09 18:38:02 +02:00
parent 312fee491c
commit 50998fcfb9
270 changed files with 9 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

1215
tools/BER Plots.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

1803
tools/Phase Measurement Prototype.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

235
tools/ROCOF test data generator.ipynb Normal file
View file

@ -0,0 +1,235 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ROCOF test waveform library\n",
"\n",
"This is a re-implementation of the ROCOF test waveforms described in https://zenodo.org/record/3559798\n",
"\n",
"**This file is exported as a python module and loaded from other notebooks here, so please make sure to re-export when changing it.**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import itertools\n",
"\n",
"import numpy as np\n",
"from scipy import signal\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def sample_waveform(generator, duration:\"s\"=10, sampling_rate:\"sp/s\"=10000, frequency:\"Hz\"=50):\n",
" samples = int(duration*sampling_rate)\n",
" phases = np.linspace(0, 2*np.pi, 6, endpoint=False)\n",
" omega_t = np.linspace(phases, phases + 2*np.pi*duration*frequency, samples)\n",
" fundamental = np.sin(omega_t)\n",
" return generator(omega_t, fundamental, sampling_rate=sampling_rate, duration=duration, frequency=frequency).swapaxes(0, 1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def gen_harmonics(amplitudes, phases=[]):\n",
" return lambda omega_t, fundamental, **_: fundamental + np.sum([\n",
" a*np.sin((p if p else 0) + i*omega_t)\n",
" for i, (a, p) in enumerate(itertools.zip_longest(amplitudes, phases), start=2)\n",
" ], axis=0)\n",
"\n",
"def test_harmonics():\n",
" return gen_harmonics([0.02, 0.05, 0.01, 0.06, 0.005, 0.05, 0.005, 0.015, 0.005, 0.035, 0.005, 0.003])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def gen_interharmonic(amplitudes, ih=[], ih_phase=[]):\n",
" def gen(omega_t, fundamental, **_):\n",
" return fundamental + np.sum([\n",
" a*np.sin(omega_t * ih + (p if p else 0))\n",
" for a, ih, p in itertools.zip_longest(amplitudes, ih, ih_phase)\n",
" ], axis=0)\n",
" return gen\n",
"\n",
"def test_interharmonics():\n",
" return gen_interharmonic([0.1], [15.01401], [np.pi])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def gen_noise(amplitude=0.2, fmax:'Hz'=4.9e3, fmin:'Hz'=100, filter_order=6):\n",
" def gen(omega_t, fundamental, sampling_rate, **_):\n",
" noise = np.random.normal(0, amplitude, fundamental.shape)\n",
" b, a = signal.butter(filter_order,\n",
" [fmin, min(fmax, sampling_rate//2-1)],\n",
" btype='bandpass',\n",
" fs=sampling_rate)\n",
" return fundamental + signal.lfilter(b, a, noise, axis=0)\n",
" return gen\n",
"\n",
"def test_noise():\n",
" return gen_noise()\n",
"\n",
"def test_noise_loud():\n",
" return gen_noise(amplitude=0.5, fmin=10)"
]
},
{
"cell_type": "code",
"execution_count": 406,
"metadata": {},
"outputs": [],
"source": [
"def gen_steps(size_amplitude=0.1, size_phase=0.1*np.pi, steps_per_sec=1):\n",
" def gen(omega_t, fundamental, duration, **_):\n",
" n = int(steps_per_sec * duration)\n",
" indices = np.random.randint(0, len(omega_t), n)\n",
" amplitudes = np.random.normal(1, size_amplitude, (n, 6))\n",
" phases = np.random.normal(0, size_phase, (n, 6))\n",
" amplitude = np.ones(omega_t.shape)\n",
" for start, end, a, p in zip(indices, indices[1:], amplitudes, phases):\n",
" omega_t[start:end] += p\n",
" amplitude[start:end] = a\n",
" return amplitude*np.sin(omega_t)\n",
" return gen\n",
"\n",
"def test_amplitude_steps():\n",
" return gen_steps(size_amplitude=0.4, size_phase=0)\n",
"\n",
"def test_phase_steps():\n",
" return gen_steps(size_amplitude=0, size_phase=0.1)\n",
"\n",
"def test_amplitude_and_phase_steps():\n",
" return gen_steps(size_amplitude=0.2, size_phase=0.07)"
]
},
{
"cell_type": "code",
"execution_count": 418,
"metadata": {},
"outputs": [],
"source": [
"def step_gen(shape, stdev, duration, steps_per_sec=1.0, mean=0.0):\n",
" samples, channels = shape\n",
" n = int(steps_per_sec * duration)\n",
" indices = np.random.randint(0, samples, n)\n",
" phases = np.random.normal(mean, stdev, (n, 6))\n",
" amplitude = np.ones((samples, channels))\n",
" out = np.zeros(shape)\n",
" for start, end, a in zip(indices, indices[1:], amplitude):\n",
" out[start:end] = a\n",
" return out\n",
"\n",
"def gen_chirp(fmin, fmax, period, dwell_time=1.0, amplitude=None, phase_steps=None):\n",
" def gen(omega_t, fundamental, sampling_rate, duration, **_):\n",
" samples = int(duration*sampling_rate)\n",
" phases = np.linspace(0, 2*np.pi, 6, endpoint=False)\n",
" \n",
" c = (fmax-fmin)/period\n",
" t = np.linspace(0, duration, samples)\n",
" \n",
" x = np.repeat(np.reshape(2*np.pi*fmin*t, (-1,1)), 6, axis=1)\n",
" data = (phases + x)[:int(sampling_rate*dwell_time)]\n",
" current_phase = 2*np.pi*fmin*dwell_time\n",
" direction = 'up'\n",
" \n",
" for idx in range(int(dwell_time*sampling_rate), samples, int(2*period*sampling_rate)):\n",
" t1 = np.linspace(0, period, int(period*sampling_rate))\n",
" t2 = np.linspace(0, period, int(period*sampling_rate))\n",
" chirp_phase = np.hstack((\n",
" 2*np.pi*(c/2 * t1**2 + fmin * t1),\n",
" 2*np.pi*(-c/2 * t2**2 + fmax * t2 - (c/2 * period**2 + fmin * period))\n",
" ))\n",
" chirp_phase = np.repeat(np.reshape(chirp_phase, (-1, 1)), 6, axis=1)\n",
" new = phases + chirp_phase + current_phase\n",
" current_phase = chirp_phase[-1]\n",
" data = np.vstack((data, new))\n",
" \n",
" data = data[:len(fundamental)]\n",
" \n",
" if phase_steps:\n",
" (step_amplitude, steps_per_sec) = phase_steps\n",
" steps = step_gen(data.shape, step_amplitude, duration, steps_per_sec)\n",
" data += steps\n",
" \n",
" if amplitude is None:\n",
" return np.sin(data)\n",
" else:\n",
" return fundamental + amplitude*np.sin(data)\n",
" return gen\n",
"\n",
"def test_close_interharmonics_and_flicker():\n",
" return gen_chirp(90.0, 150.0, 10, 1, amplitude=0.1)\n",
"\n",
"def test_off_frequency():\n",
"# return gen_chirp(48.0, 52.0, 0.25, 1)\n",
" return gen_chirp(48.0, 52.0, 10, 1)\n",
"\n",
"def test_sweep_phase_steps():\n",
" return gen_chirp(48.0, 52.0, 10, 1, phase_steps=(0.1, 1))\n",
"# return gen_chirp(48.0, 52.0, 0.25, 1, phase_steps=(0.1, 1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"all_tests = [test_harmonics, test_interharmonics, test_noise, test_noise_loud, test_amplitude_steps, test_phase_steps, test_amplitude_and_phase_steps, test_close_interharmonics_and_flicker, test_off_frequency, test_sweep_phase_steps]"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

220
tools/dec_proto_am_ber_top.py Executable file
View file

@ -0,0 +1,220 @@
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Dec Proto Am Ber Top
# GNU Radio version: 3.7.13.5
##################################################
from gnuradio import analog
from gnuradio import blocks
from gnuradio import digital
from gnuradio import eng_notation
from gnuradio import fec
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser
import pmt
class dec_proto_am_ber_top(gr.top_block):
def __init__(self, ber_file='0', carrier=1, mod_depth=0.8, signal_strength=2):
gr.top_block.__init__(self, "Dec Proto Am Ber Top")
##################################################
# Parameters
##################################################
self.ber_file = ber_file
self.carrier = carrier
self.mod_depth = mod_depth
self.signal_strength = signal_strength
##################################################
# Variables
##################################################
self.sim_mul = sim_mul = 1e4
self.actual_sampling_rate = actual_sampling_rate = 10
self.sync_tag = sync_tag = gr.tag_utils.python_to_tag((0, pmt.intern("sync"), pmt.from_double(0.0), pmt.intern("correlate_access_code")))
self.samp_rate = samp_rate = actual_sampling_rate*sim_mul
self.pi = pi = 3.141592653589793
self.packet_time_est_tag = packet_time_est_tag = gr.tag_utils.python_to_tag((0, pmt.intern("start"), pmt.from_double(0.0), pmt.intern("packet_vector_source")))
##################################################
# Blocks
##################################################
self.low_pass_filter_0 = filter.fir_filter_ccf(1, firdes.low_pass(
1.0/signal_strength*1000 * 2, samp_rate, 0.2 * sim_mul, 0.05 * sim_mul, firdes.WIN_HAMMING, 6.76))
self.fec_ber_bf_0 = fec.ber_bf(False, 0, -7.0)
self.digital_clock_recovery_mm_xx_0 = digital.clock_recovery_mm_ff(50, 0.001, 0, 0.01, 0.01)
self.digital_binary_slicer_fb_0 = digital.binary_slicer_fb()
self.blocks_vector_source_x_0_0_1_0 = blocks.vector_source_f([1,0]*(4*12)+[1,1,0,1,0,1,0,1]*12+[1,0,1,1,1,1,1,0,0,1]+[1,1,1,1,0,1,1,0,0,1]+[1,0,1,1,1,1,1,0,0,1]+[0,1,1,1,0,1,1,0,1,0]+[0,0,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0]+[0]*128, True, 1, [packet_time_est_tag])
self.blocks_throttle_0 = blocks.throttle(gr.sizeof_float*1, samp_rate,True)
self.blocks_repeat_0 = blocks.repeat(gr.sizeof_float*1, 10*5)
self.blocks_null_source_1 = blocks.null_source(gr.sizeof_float*1)
self.blocks_multiply_xx_0_0 = blocks.multiply_vff(1)
self.blocks_multiply_xx_0 = blocks.multiply_vcc(1)
self.blocks_multiply_const_vxx_1 = blocks.multiply_const_vff((mod_depth, ))
self.blocks_multiply_const_vxx_0 = blocks.multiply_const_vcc((1.0/mod_depth, ))
self.blocks_float_to_complex_0 = blocks.float_to_complex(1)
self.blocks_float_to_char_0 = blocks.float_to_char(1, 1)
self.blocks_file_source_0 = blocks.file_source(gr.sizeof_float*1, '/home/user/research/smart_meter_reset/gm_platform/fw/raw_freq.bin', True)
self.blocks_file_source_0.set_begin_tag(pmt.PMT_NIL)
self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, ber_file, False)
self.blocks_file_sink_0.set_unbuffered(False)
self.blocks_delay_0 = blocks.delay(gr.sizeof_float*1, 5)
self.blocks_complex_to_mag_0 = blocks.complex_to_mag(1)
self.blocks_add_xx_0 = blocks.add_vff(1)
self.blocks_add_const_vxx_2_0 = blocks.add_const_vff((-0.5, ))
self.blocks_add_const_vxx_2 = blocks.add_const_vcc((-(1-mod_depth), ))
self.blocks_add_const_vxx_1 = blocks.add_const_vff((1-mod_depth, ))
self.blocks_add_const_vxx_0 = blocks.add_const_vff((-50, ))
self.analog_sig_source_x_0_0 = analog.sig_source_f(samp_rate, analog.GR_COS_WAVE, carrier * sim_mul, signal_strength*1e-3, 0)
self.analog_sig_source_x_0 = analog.sig_source_c(samp_rate, analog.GR_COS_WAVE, -carrier * sim_mul, 1, 0)
##################################################
# Connections
##################################################
self.connect((self.analog_sig_source_x_0, 0), (self.blocks_multiply_xx_0, 1))
self.connect((self.analog_sig_source_x_0_0, 0), (self.blocks_multiply_xx_0_0, 1))
self.connect((self.blocks_add_const_vxx_0, 0), (self.blocks_add_xx_0, 0))
self.connect((self.blocks_add_const_vxx_1, 0), (self.blocks_multiply_xx_0_0, 0))
self.connect((self.blocks_add_const_vxx_2, 0), (self.blocks_multiply_const_vxx_0, 0))
self.connect((self.blocks_add_const_vxx_2_0, 0), (self.digital_clock_recovery_mm_xx_0, 0))
self.connect((self.blocks_add_xx_0, 0), (self.blocks_float_to_complex_0, 0))
self.connect((self.blocks_complex_to_mag_0, 0), (self.blocks_add_const_vxx_2_0, 0))
self.connect((self.blocks_delay_0, 0), (self.blocks_float_to_char_0, 0))
self.connect((self.blocks_file_source_0, 0), (self.blocks_throttle_0, 0))
self.connect((self.blocks_float_to_char_0, 0), (self.fec_ber_bf_0, 1))
self.connect((self.blocks_float_to_complex_0, 0), (self.blocks_multiply_xx_0, 0))
self.connect((self.blocks_multiply_const_vxx_0, 0), (self.blocks_complex_to_mag_0, 0))
self.connect((self.blocks_multiply_const_vxx_1, 0), (self.blocks_add_const_vxx_1, 0))
self.connect((self.blocks_multiply_xx_0, 0), (self.low_pass_filter_0, 0))
self.connect((self.blocks_multiply_xx_0_0, 0), (self.blocks_add_xx_0, 1))
self.connect((self.blocks_null_source_1, 0), (self.blocks_float_to_complex_0, 1))
self.connect((self.blocks_repeat_0, 0), (self.blocks_multiply_const_vxx_1, 0))
self.connect((self.blocks_throttle_0, 0), (self.blocks_add_const_vxx_0, 0))
self.connect((self.blocks_vector_source_x_0_0_1_0, 0), (self.blocks_delay_0, 0))
self.connect((self.blocks_vector_source_x_0_0_1_0, 0), (self.blocks_repeat_0, 0))
self.connect((self.digital_binary_slicer_fb_0, 0), (self.fec_ber_bf_0, 0))
self.connect((self.digital_clock_recovery_mm_xx_0, 0), (self.digital_binary_slicer_fb_0, 0))
self.connect((self.fec_ber_bf_0, 0), (self.blocks_file_sink_0, 0))
self.connect((self.low_pass_filter_0, 0), (self.blocks_add_const_vxx_2, 0))
def get_ber_file(self):
return self.ber_file
def set_ber_file(self, ber_file):
self.ber_file = ber_file
self.blocks_file_sink_0.open(self.ber_file)
def get_carrier(self):
return self.carrier
def set_carrier(self, carrier):
self.carrier = carrier
self.analog_sig_source_x_0_0.set_frequency(self.carrier * self.sim_mul)
self.analog_sig_source_x_0.set_frequency(-self.carrier * self.sim_mul)
def get_mod_depth(self):
return self.mod_depth
def set_mod_depth(self, mod_depth):
self.mod_depth = mod_depth
self.blocks_multiply_const_vxx_1.set_k((self.mod_depth, ))
self.blocks_multiply_const_vxx_0.set_k((1.0/self.mod_depth, ))
self.blocks_add_const_vxx_2.set_k((-(1-self.mod_depth), ))
self.blocks_add_const_vxx_1.set_k((1-self.mod_depth, ))
def get_signal_strength(self):
return self.signal_strength
def set_signal_strength(self, signal_strength):
self.signal_strength = signal_strength
self.low_pass_filter_0.set_taps(firdes.low_pass(1.0/self.signal_strength*1000 * 2, self.samp_rate, 0.2 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.analog_sig_source_x_0_0.set_amplitude(self.signal_strength*1e-3)
def get_sim_mul(self):
return self.sim_mul
def set_sim_mul(self, sim_mul):
self.sim_mul = sim_mul
self.set_samp_rate(self.actual_sampling_rate*self.sim_mul)
self.low_pass_filter_0.set_taps(firdes.low_pass(1.0/self.signal_strength*1000 * 2, self.samp_rate, 0.2 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.analog_sig_source_x_0_0.set_frequency(self.carrier * self.sim_mul)
self.analog_sig_source_x_0.set_frequency(-self.carrier * self.sim_mul)
def get_actual_sampling_rate(self):
return self.actual_sampling_rate
def set_actual_sampling_rate(self, actual_sampling_rate):
self.actual_sampling_rate = actual_sampling_rate
self.set_samp_rate(self.actual_sampling_rate*self.sim_mul)
def get_sync_tag(self):
return self.sync_tag
def set_sync_tag(self, sync_tag):
self.sync_tag = sync_tag
def get_samp_rate(self):
return self.samp_rate
def set_samp_rate(self, samp_rate):
self.samp_rate = samp_rate
self.low_pass_filter_0.set_taps(firdes.low_pass(1.0/self.signal_strength*1000 * 2, self.samp_rate, 0.2 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.blocks_throttle_0.set_sample_rate(self.samp_rate)
self.analog_sig_source_x_0_0.set_sampling_freq(self.samp_rate)
self.analog_sig_source_x_0.set_sampling_freq(self.samp_rate)
def get_pi(self):
return self.pi
def set_pi(self, pi):
self.pi = pi
def get_packet_time_est_tag(self):
return self.packet_time_est_tag
def set_packet_time_est_tag(self, packet_time_est_tag):
self.packet_time_est_tag = packet_time_est_tag
self.blocks_vector_source_x_0_0_1_0.set_data([1,0]*(4*12)+[1,1,0,1,0,1,0,1]*12+[1,0,1,1,1,1,1,0,0,1]+[1,1,1,1,0,1,1,0,0,1]+[1,0,1,1,1,1,1,0,0,1]+[0,1,1,1,0,1,1,0,1,0]+[0,0,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0]+[0]*128, [self.packet_time_est_tag])
def argument_parser():
parser = OptionParser(usage="%prog: [options]", option_class=eng_option)
parser.add_option(
"", "--ber-file", dest="ber_file", type="string", default='0',
help="Set BER data output file [default=%default]")
parser.add_option(
"", "--carrier", dest="carrier", type="eng_float", default=eng_notation.num_to_str(1),
help="Set Carrier frequency in Hz [default=%default]")
parser.add_option(
"", "--mod-depth", dest="mod_depth", type="eng_float", default=eng_notation.num_to_str(0.8),
help="Set Modulation depth (0-1) [default=%default]")
parser.add_option(
"", "--signal-strength", dest="signal_strength", type="eng_float", default=eng_notation.num_to_str(2),
help="Set signal strength in mHz [default=%default]")
return parser
def main(top_block_cls=dec_proto_am_ber_top, options=None):
if options is None:
options, _ = argument_parser().parse_args()
tb = top_block_cls(ber_file=options.ber_file, carrier=options.carrier, mod_depth=options.mod_depth, signal_strength=options.signal_strength)
tb.start()
try:
raw_input('Press Enter to quit: ')
except EOFError:
pass
tb.stop()
tb.wait()
if __name__ == '__main__':
main()

191
tools/dec_proto_am_dc_ber_top.py Executable file
View file

@ -0,0 +1,191 @@
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Dec Proto Am Dc Ber Top
# GNU Radio version: 3.7.13.5
##################################################
from gnuradio import blocks
from gnuradio import digital
from gnuradio import eng_notation
from gnuradio import fec
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser
import pmt
class dec_proto_am_dc_ber_top(gr.top_block):
def __init__(self, ber_file='0', signal_strength=50):
gr.top_block.__init__(self, "Dec Proto Am Dc Ber Top")
##################################################
# Parameters
##################################################
self.ber_file = ber_file
self.signal_strength = signal_strength
##################################################
# Variables
##################################################
self.sim_mul = sim_mul = 1e4
self.actual_sampling_rate = actual_sampling_rate = 10
self.sync_tag = sync_tag = gr.tag_utils.python_to_tag((0, pmt.intern("sync"), pmt.from_double(0.0), pmt.intern("correlate_access_code")))
self.samp_rate = samp_rate = actual_sampling_rate*sim_mul
self.pi = pi = 3.141592653589793
self.packet_time_est_tag = packet_time_est_tag = gr.tag_utils.python_to_tag((0, pmt.intern("start"), pmt.from_double(0.0), pmt.intern("packet_vector_source")))
self.ber_delay = ber_delay = 198
##################################################
# Blocks
##################################################
self.low_pass_filter_0 = filter.fir_filter_ccf(1, firdes.low_pass(
1, samp_rate, 0.1 * sim_mul, 0.05 * sim_mul, firdes.WIN_HAMMING, 6.76))
self.high_pass_filter_0 = filter.fir_filter_ccf(1, firdes.high_pass(
1, samp_rate, sim_mul/200, sim_mul/800, firdes.WIN_HAMMING, 6.76))
self.fec_ber_bf_0 = fec.ber_bf(False, 0, -7.0)
self.digital_clock_recovery_mm_xx_0 = digital.clock_recovery_mm_ff(50, 0.001, 0, 0.01, 0.01)
self.digital_binary_slicer_fb_0 = digital.binary_slicer_fb()
self.blocks_vector_source_x_0_0_1_0 = blocks.vector_source_f([1,0]*(4*12)+[1,1,0,1,0,1,0,1]*12+[1,0,1,1,1,1,1,0,0,1]+[1,1,1,1,0,1,1,0,0,1]+[1,0,1,1,1,1,1,0,0,1]+[0,1,1,1,0,1,1,0,1,0]+[0,0,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0]+[0]*8, True, 1, [packet_time_est_tag])
self.blocks_throttle_0 = blocks.throttle(gr.sizeof_float*1, samp_rate,True)
self.blocks_repeat_0 = blocks.repeat(gr.sizeof_float*1, 10*5)
self.blocks_null_source_1 = blocks.null_source(gr.sizeof_float*1)
self.blocks_multiply_const_vxx_1 = blocks.multiply_const_vff((signal_strength * 0.001, ))
self.blocks_multiply_const_vxx_0 = blocks.multiply_const_vff((1000.0/signal_strength, ))
self.blocks_float_to_complex_0 = blocks.float_to_complex(1)
self.blocks_float_to_char_0 = blocks.float_to_char(1, 1)
self.blocks_file_source_0 = blocks.file_source(gr.sizeof_float*1, '/home/user/research/smart_meter_reset/gm_platform/fw/raw_freq.bin', True)
self.blocks_file_source_0.set_begin_tag(pmt.PMT_NIL)
self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, ber_file, False)
self.blocks_file_sink_0.set_unbuffered(False)
self.blocks_delay_0 = blocks.delay(gr.sizeof_float*1, ber_delay)
self.blocks_complex_to_real_0 = blocks.complex_to_real(1)
self.blocks_add_xx_0 = blocks.add_vff(1)
self.blocks_add_const_vxx_1 = blocks.add_const_vff((-signal_strength*0.001/2.0, ))
self.blocks_add_const_vxx_0 = blocks.add_const_vff((-50, ))
##################################################
# Connections
##################################################
self.connect((self.blocks_add_const_vxx_0, 0), (self.blocks_add_xx_0, 0))
self.connect((self.blocks_add_const_vxx_1, 0), (self.blocks_add_xx_0, 1))
self.connect((self.blocks_add_xx_0, 0), (self.blocks_float_to_complex_0, 0))
self.connect((self.blocks_complex_to_real_0, 0), (self.blocks_multiply_const_vxx_0, 0))
self.connect((self.blocks_delay_0, 0), (self.blocks_float_to_char_0, 0))
self.connect((self.blocks_file_source_0, 0), (self.blocks_throttle_0, 0))
self.connect((self.blocks_float_to_char_0, 0), (self.fec_ber_bf_0, 1))
self.connect((self.blocks_float_to_complex_0, 0), (self.high_pass_filter_0, 0))
self.connect((self.blocks_multiply_const_vxx_0, 0), (self.digital_clock_recovery_mm_xx_0, 0))
self.connect((self.blocks_multiply_const_vxx_1, 0), (self.blocks_add_const_vxx_1, 0))
self.connect((self.blocks_null_source_1, 0), (self.blocks_float_to_complex_0, 1))
self.connect((self.blocks_repeat_0, 0), (self.blocks_multiply_const_vxx_1, 0))
self.connect((self.blocks_throttle_0, 0), (self.blocks_add_const_vxx_0, 0))
self.connect((self.blocks_vector_source_x_0_0_1_0, 0), (self.blocks_delay_0, 0))
self.connect((self.blocks_vector_source_x_0_0_1_0, 0), (self.blocks_repeat_0, 0))
self.connect((self.digital_binary_slicer_fb_0, 0), (self.fec_ber_bf_0, 0))
self.connect((self.digital_clock_recovery_mm_xx_0, 0), (self.digital_binary_slicer_fb_0, 0))
self.connect((self.fec_ber_bf_0, 0), (self.blocks_file_sink_0, 0))
self.connect((self.high_pass_filter_0, 0), (self.low_pass_filter_0, 0))
self.connect((self.low_pass_filter_0, 0), (self.blocks_complex_to_real_0, 0))
def get_ber_file(self):
return self.ber_file
def set_ber_file(self, ber_file):
self.ber_file = ber_file
self.blocks_file_sink_0.open(self.ber_file)
def get_signal_strength(self):
return self.signal_strength
def set_signal_strength(self, signal_strength):
self.signal_strength = signal_strength
self.blocks_multiply_const_vxx_1.set_k((self.signal_strength * 0.001, ))
self.blocks_multiply_const_vxx_0.set_k((1000.0/self.signal_strength, ))
self.blocks_add_const_vxx_1.set_k((-self.signal_strength*0.001/2.0, ))
def get_sim_mul(self):
return self.sim_mul
def set_sim_mul(self, sim_mul):
self.sim_mul = sim_mul
self.set_samp_rate(self.actual_sampling_rate*self.sim_mul)
self.low_pass_filter_0.set_taps(firdes.low_pass(1, self.samp_rate, 0.1 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.high_pass_filter_0.set_taps(firdes.high_pass(1, self.samp_rate, self.sim_mul/200, self.sim_mul/800, firdes.WIN_HAMMING, 6.76))
def get_actual_sampling_rate(self):
return self.actual_sampling_rate
def set_actual_sampling_rate(self, actual_sampling_rate):
self.actual_sampling_rate = actual_sampling_rate
self.set_samp_rate(self.actual_sampling_rate*self.sim_mul)
def get_sync_tag(self):
return self.sync_tag
def set_sync_tag(self, sync_tag):
self.sync_tag = sync_tag
def get_samp_rate(self):
return self.samp_rate
def set_samp_rate(self, samp_rate):
self.samp_rate = samp_rate
self.low_pass_filter_0.set_taps(firdes.low_pass(1, self.samp_rate, 0.1 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.high_pass_filter_0.set_taps(firdes.high_pass(1, self.samp_rate, self.sim_mul/200, self.sim_mul/800, firdes.WIN_HAMMING, 6.76))
self.blocks_throttle_0.set_sample_rate(self.samp_rate)
def get_pi(self):
return self.pi
def set_pi(self, pi):
self.pi = pi
def get_packet_time_est_tag(self):
return self.packet_time_est_tag
def set_packet_time_est_tag(self, packet_time_est_tag):
self.packet_time_est_tag = packet_time_est_tag
self.blocks_vector_source_x_0_0_1_0.set_data([1,0]*(4*12)+[1,1,0,1,0,1,0,1]*12+[1,0,1,1,1,1,1,0,0,1]+[1,1,1,1,0,1,1,0,0,1]+[1,0,1,1,1,1,1,0,0,1]+[0,1,1,1,0,1,1,0,1,0]+[0,0,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0]+[0]*8, [self.packet_time_est_tag])
def get_ber_delay(self):
return self.ber_delay
def set_ber_delay(self, ber_delay):
self.ber_delay = ber_delay
self.blocks_delay_0.set_dly(self.ber_delay)
def argument_parser():
parser = OptionParser(usage="%prog: [options]", option_class=eng_option)
parser.add_option(
"", "--ber-file", dest="ber_file", type="string", default='0',
help="Set BER data output file [default=%default]")
parser.add_option(
"", "--signal-strength", dest="signal_strength", type="eng_float", default=eng_notation.num_to_str(50),
help="Set signal strength in mHz [default=%default]")
return parser
def main(top_block_cls=dec_proto_am_dc_ber_top, options=None):
if options is None:
options, _ = argument_parser().parse_args()
tb = top_block_cls(ber_file=options.ber_file, signal_strength=options.signal_strength)
tb.start()
try:
raw_input('Press Enter to quit: ')
except EOFError:
pass
tb.stop()
tb.wait()
if __name__ == '__main__':
main()

183
tools/dec_proto_fm_ber_top.py Executable file
View file

@ -0,0 +1,183 @@
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Dec Proto Fm Ber Top
# GNU Radio version: 3.7.13.5
##################################################
from gnuradio import analog
from gnuradio import blocks
from gnuradio import digital
from gnuradio import eng_notation
from gnuradio import fec
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser
import math
import pmt
class dec_proto_fm_ber_top(gr.top_block):
def __init__(self, ber_file='0', signal_strength=1):
gr.top_block.__init__(self, "Dec Proto Fm Ber Top")
##################################################
# Parameters
##################################################
self.ber_file = ber_file
self.signal_strength = signal_strength
##################################################
# Variables
##################################################
self.sim_mul = sim_mul = 1e4
self.actual_sampling_rate = actual_sampling_rate = 10
self.sync_tag = sync_tag = gr.tag_utils.python_to_tag((0, pmt.intern("sync"), pmt.from_double(0.0), pmt.intern("correlate_access_code")))
self.samp_rate = samp_rate = actual_sampling_rate*sim_mul
self.pi = pi = 3.141592653589793
self.packet_time_est_tag = packet_time_est_tag = gr.tag_utils.python_to_tag((0, pmt.intern("start"), pmt.from_double(0.0), pmt.intern("packet_vector_source")))
##################################################
# Blocks
##################################################
self.low_pass_filter_0 = filter.fir_filter_ccf(1, firdes.low_pass(
1, samp_rate, 0.1 * sim_mul, 0.05 * sim_mul, firdes.WIN_HAMMING, 6.76))
self.fec_ber_bf_0 = fec.ber_bf(False, 0, -7.0)
self.digital_clock_recovery_mm_xx_0 = digital.clock_recovery_mm_ff(50, 0.001, 0, 0.01, 0.01)
self.digital_binary_slicer_fb_0 = digital.binary_slicer_fb()
self.blocks_vector_source_x_0_0_1_0 = blocks.vector_source_f([1,0]*(4*12)+[1,1,0,1,0,1,0,1]*12+[1,0,1,1,1,1,1,0,0,1]+[1,1,1,1,0,1,1,0,0,1]+[1,0,1,1,1,1,1,0,0,1]+[0,1,1,1,0,1,1,0,1,0]+[0,0,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0]+[0]*128, True, 1, [packet_time_est_tag])
self.blocks_vco_f_0 = blocks.vco_f(samp_rate, sim_mul*2*pi, signal_strength*1e-3)
self.blocks_throttle_0 = blocks.throttle(gr.sizeof_float*1, samp_rate,True)
self.blocks_repeat_0 = blocks.repeat(gr.sizeof_float*1, 10*5)
self.blocks_null_source_1 = blocks.null_source(gr.sizeof_float*1)
self.blocks_multiply_xx_0 = blocks.multiply_vcc(1)
self.blocks_multiply_const_vxx_1 = blocks.multiply_const_vff((0.2, ))
self.blocks_float_to_complex_0 = blocks.float_to_complex(1)
self.blocks_float_to_char_0 = blocks.float_to_char(1, 1)
self.blocks_file_source_0 = blocks.file_source(gr.sizeof_float*1, '/home/user/research/smart_meter_reset/gm_platform/fw/raw_freq.bin', True)
self.blocks_file_source_0.set_begin_tag(pmt.PMT_NIL)
self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, ber_file, False)
self.blocks_file_sink_0.set_unbuffered(False)
self.blocks_delay_0 = blocks.delay(gr.sizeof_float*1, 5)
self.blocks_add_xx_0 = blocks.add_vff(1)
self.blocks_add_const_vxx_1 = blocks.add_const_vff((1.0, ))
self.blocks_add_const_vxx_0 = blocks.add_const_vff((-50, ))
self.analog_sig_source_x_0 = analog.sig_source_c(samp_rate, analog.GR_COS_WAVE, -1.1 * sim_mul, 1, 0)
self.analog_quadrature_demod_cf_0 = analog.quadrature_demod_cf(8)
##################################################
# Connections
##################################################
self.connect((self.analog_quadrature_demod_cf_0, 0), (self.digital_clock_recovery_mm_xx_0, 0))
self.connect((self.analog_sig_source_x_0, 0), (self.blocks_multiply_xx_0, 1))
self.connect((self.blocks_add_const_vxx_0, 0), (self.blocks_add_xx_0, 0))
self.connect((self.blocks_add_const_vxx_1, 0), (self.blocks_vco_f_0, 0))
self.connect((self.blocks_add_xx_0, 0), (self.blocks_float_to_complex_0, 0))
self.connect((self.blocks_delay_0, 0), (self.blocks_float_to_char_0, 0))
self.connect((self.blocks_file_source_0, 0), (self.blocks_throttle_0, 0))
self.connect((self.blocks_float_to_char_0, 0), (self.fec_ber_bf_0, 1))
self.connect((self.blocks_float_to_complex_0, 0), (self.blocks_multiply_xx_0, 0))
self.connect((self.blocks_multiply_const_vxx_1, 0), (self.blocks_add_const_vxx_1, 0))
self.connect((self.blocks_multiply_xx_0, 0), (self.low_pass_filter_0, 0))
self.connect((self.blocks_null_source_1, 0), (self.blocks_float_to_complex_0, 1))
self.connect((self.blocks_repeat_0, 0), (self.blocks_multiply_const_vxx_1, 0))
self.connect((self.blocks_throttle_0, 0), (self.blocks_add_const_vxx_0, 0))
self.connect((self.blocks_vco_f_0, 0), (self.blocks_add_xx_0, 1))
self.connect((self.blocks_vector_source_x_0_0_1_0, 0), (self.blocks_delay_0, 0))
self.connect((self.blocks_vector_source_x_0_0_1_0, 0), (self.blocks_repeat_0, 0))
self.connect((self.digital_binary_slicer_fb_0, 0), (self.fec_ber_bf_0, 0))
self.connect((self.digital_clock_recovery_mm_xx_0, 0), (self.digital_binary_slicer_fb_0, 0))
self.connect((self.fec_ber_bf_0, 0), (self.blocks_file_sink_0, 0))
self.connect((self.low_pass_filter_0, 0), (self.analog_quadrature_demod_cf_0, 0))
def get_ber_file(self):
return self.ber_file
def set_ber_file(self, ber_file):
self.ber_file = ber_file
self.blocks_file_sink_0.open(self.ber_file)
def get_signal_strength(self):
return self.signal_strength
def set_signal_strength(self, signal_strength):
self.signal_strength = signal_strength
def get_sim_mul(self):
return self.sim_mul
def set_sim_mul(self, sim_mul):
self.sim_mul = sim_mul
self.set_samp_rate(self.actual_sampling_rate*self.sim_mul)
self.low_pass_filter_0.set_taps(firdes.low_pass(1, self.samp_rate, 0.1 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.analog_sig_source_x_0.set_frequency(-1.1 * self.sim_mul)
def get_actual_sampling_rate(self):
return self.actual_sampling_rate
def set_actual_sampling_rate(self, actual_sampling_rate):
self.actual_sampling_rate = actual_sampling_rate
self.set_samp_rate(self.actual_sampling_rate*self.sim_mul)
def get_sync_tag(self):
return self.sync_tag
def set_sync_tag(self, sync_tag):
self.sync_tag = sync_tag
def get_samp_rate(self):
return self.samp_rate
def set_samp_rate(self, samp_rate):
self.samp_rate = samp_rate
self.low_pass_filter_0.set_taps(firdes.low_pass(1, self.samp_rate, 0.1 * self.sim_mul, 0.05 * self.sim_mul, firdes.WIN_HAMMING, 6.76))
self.blocks_throttle_0.set_sample_rate(self.samp_rate)
self.analog_sig_source_x_0.set_sampling_freq(self.samp_rate)
def get_pi(self):
return self.pi
def set_pi(self, pi):
self.pi = pi
def get_packet_time_est_tag(self):
return self.packet_time_est_tag
def set_packet_time_est_tag(self, packet_time_est_tag):
self.packet_time_est_tag = packet_time_est_tag
self.blocks_vector_source_x_0_0_1_0.set_data([1,0]*(4*12)+[1,1,0,1,0,1,0,1]*12+[1,0,1,1,1,1,1,0,0,1]+[1,1,1,1,0,1,1,0,0,1]+[1,0,1,1,1,1,1,0,0,1]+[0,1,1,1,0,1,1,0,1,0]+[0,0,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0]+[0]*128, [self.packet_time_est_tag])
def argument_parser():
parser = OptionParser(usage="%prog: [options]", option_class=eng_option)
parser.add_option(
"", "--ber-file", dest="ber_file", type="string", default='0',
help="Set BER data output file [default=%default]")
parser.add_option(
"", "--signal-strength", dest="signal_strength", type="eng_float", default=eng_notation.num_to_str(1),
help="Set signal strength in mHz [default=%default]")
return parser
def main(top_block_cls=dec_proto_fm_ber_top, options=None):
if options is None:
options, _ = argument_parser().parse_args()
tb = top_block_cls(ber_file=options.ber_file, signal_strength=options.signal_strength)
tb.start()
try:
raw_input('Press Enter to quit: ')
except EOFError:
pass
tb.stop()
tb.wait()
if __name__ == '__main__':
main()

4206
tools/grid_frequency_spectra.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

180
tools/rocof_test_data.py Normal file
View file

@ -0,0 +1,180 @@
#!/usr/bin/env python
# coding: utf-8
# # ROCOF test waveform library
#
# This is a re-implementation of the ROCOF test waveforms described in https://zenodo.org/record/3559798
#
# **This file is exported as a python module and loaded from other notebooks here, so please make sure to re-export when changing it.**
# In[ ]:
import math
import itertools
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt
# In[ ]:
get_ipython().run_line_magic('matplotlib', 'notebook')
# In[ ]:
def sample_waveform(generator, duration:"s"=10, sampling_rate:"sp/s"=10000, frequency:"Hz"=50):
samples = int(duration*sampling_rate)
phases = np.linspace(0, 2*np.pi, 6, endpoint=False)
omega_t = np.linspace(phases, phases + 2*np.pi*duration*frequency, samples)
fundamental = np.sin(omega_t)
return generator(omega_t, fundamental, sampling_rate=sampling_rate, duration=duration, frequency=frequency).swapaxes(0, 1)
# In[ ]:
def gen_harmonics(amplitudes, phases=[]):
return lambda omega_t, fundamental, **_: fundamental + np.sum([
a*np.sin((p if p else 0) + i*omega_t)
for i, (a, p) in enumerate(itertools.zip_longest(amplitudes, phases), start=2)
], axis=0)
def test_harmonics():
return gen_harmonics([0.02, 0.05, 0.01, 0.06, 0.005, 0.05, 0.005, 0.015, 0.005, 0.035, 0.005, 0.003])
# In[ ]:
def gen_interharmonic(amplitudes, ih=[], ih_phase=[]):
def gen(omega_t, fundamental, **_):
return fundamental + np.sum([
a*np.sin(omega_t * ih + (p if p else 0))
for a, ih, p in itertools.zip_longest(amplitudes, ih, ih_phase)
], axis=0)
return gen
def test_interharmonics():
return gen_interharmonic([0.1], [15.01401], [np.pi])
# In[ ]:
def gen_noise(amplitude=0.2, fmax:'Hz'=4.9e3, fmin:'Hz'=100, filter_order=6):
def gen(omega_t, fundamental, sampling_rate, **_):
noise = np.random.normal(0, amplitude, fundamental.shape)
b, a = signal.butter(filter_order,
[fmin, min(fmax, sampling_rate//2-1)],
btype='bandpass',
fs=sampling_rate)
return fundamental + signal.lfilter(b, a, noise, axis=0)
return gen
def test_noise():
return gen_noise()
def test_noise_loud():
return gen_noise(amplitude=0.5, fmin=10)
# In[406]:
def gen_steps(size_amplitude=0.1, size_phase=0.1*np.pi, steps_per_sec=1):
def gen(omega_t, fundamental, duration, **_):
n = int(steps_per_sec * duration)
indices = np.random.randint(0, len(omega_t), n)
amplitudes = np.random.normal(1, size_amplitude, (n, 6))
phases = np.random.normal(0, size_phase, (n, 6))
amplitude = np.ones(omega_t.shape)
for start, end, a, p in zip(indices, indices[1:], amplitudes, phases):
omega_t[start:end] += p
amplitude[start:end] = a
return amplitude*np.sin(omega_t)
return gen
def test_amplitude_steps():
return gen_steps(size_amplitude=0.4, size_phase=0)
def test_phase_steps():
return gen_steps(size_amplitude=0, size_phase=0.1)
def test_amplitude_and_phase_steps():
return gen_steps(size_amplitude=0.2, size_phase=0.07)
# In[418]:
def step_gen(shape, stdev, duration, steps_per_sec=1.0, mean=0.0):
samples, channels = shape
n = int(steps_per_sec * duration)
indices = np.random.randint(0, samples, n)
phases = np.random.normal(mean, stdev, (n, 6))
amplitude = np.ones((samples, channels))
out = np.zeros(shape)
for start, end, a in zip(indices, indices[1:], amplitude):
out[start:end] = a
return out
def gen_chirp(fmin, fmax, period, dwell_time=1.0, amplitude=None, phase_steps=None):
def gen(omega_t, fundamental, sampling_rate, duration, **_):
samples = int(duration*sampling_rate)
phases = np.linspace(0, 2*np.pi, 6, endpoint=False)
c = (fmax-fmin)/period
t = np.linspace(0, duration, samples)
x = np.repeat(np.reshape(2*np.pi*fmin*t, (-1,1)), 6, axis=1)
data = (phases + x)[:int(sampling_rate*dwell_time)]
current_phase = 2*np.pi*fmin*dwell_time
direction = 'up'
for idx in range(int(dwell_time*sampling_rate), samples, int(2*period*sampling_rate)):
t1 = np.linspace(0, period, int(period*sampling_rate))
t2 = np.linspace(0, period, int(period*sampling_rate))
chirp_phase = np.hstack((
2*np.pi*(c/2 * t1**2 + fmin * t1),
2*np.pi*(-c/2 * t2**2 + fmax * t2 - (c/2 * period**2 + fmin * period))
))
chirp_phase = np.repeat(np.reshape(chirp_phase, (-1, 1)), 6, axis=1)
new = phases + chirp_phase + current_phase
current_phase = chirp_phase[-1]
data = np.vstack((data, new))
data = data[:len(fundamental)]
if phase_steps:
(step_amplitude, steps_per_sec) = phase_steps
steps = step_gen(data.shape, step_amplitude, duration, steps_per_sec)
data += steps
if amplitude is None:
return np.sin(data)
else:
return fundamental + amplitude*np.sin(data)
return gen
def test_close_interharmonics_and_flicker():
return gen_chirp(90.0, 150.0, 10, 1, amplitude=0.1)
def test_off_frequency():
# return gen_chirp(48.0, 52.0, 0.25, 1)
return gen_chirp(48.0, 52.0, 10, 1)
def test_sweep_phase_steps():
return gen_chirp(48.0, 52.0, 10, 1, phase_steps=(0.1, 1))
# return gen_chirp(48.0, 52.0, 0.25, 1, phase_steps=(0.1, 1))
# In[ ]:
all_tests = [test_harmonics, test_interharmonics, test_noise, test_noise_loud, test_amplitude_steps, test_phase_steps, test_amplitude_and_phase_steps, test_close_interharmonics_and_flicker, test_off_frequency, test_sweep_phase_steps]

93
tools/sweep_gr_sims.py Normal file
View file

@ -0,0 +1,93 @@
#!/usr/bin/env python3
import subprocess
import time
import struct
import random
import math
import statistics
import tempfile
import itertools
from collections import defaultdict
from os import path
import tqdm
SIMS = [ 'dec_proto_am_ber_top.py',
'dec_proto_am_dc_ber_top.py',
'dec_proto_fm_ber_top.py',
]
E12_SERIES = [1.0, 1.2, 1.5, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6, 6.8, 8.2]
AMPLITUDES_MILLIHERTZ = [ x*y for y in [1, 10, 100] for x in E12_SERIES ]
SIMULATION_DURATION = 30.0 # seconds realtime
MAX_CONCURRENT_PROCESSES = 8
def run_simulation(
duration = SIMULATION_DURATION,
simulations = SIMS,
forklimit = MAX_CONCURRENT_PROCESSES,
amplitudes:'list(millihertz)' = AMPLITUDES_MILLIHERTZ,
terminate_timeout:'s' = 5.0,
communicate_timeout:'s' = 10.0,
repeat_runs = 1,
shuffle = False,
tqdm = tqdm.tqdm
):
with tempfile.TemporaryDirectory() as tmpdir:
jobs = list(enumerate(itertools.product(simulations, amplitudes * repeat_runs)))
if shuffle:
random.shuffle(jobs)
nchunks = int(math.ceil(len(jobs)/forklimit))
def start_processes(jobs):
for i, (sim, ampl_mhz) in jobs:
berfile = path.join(tmpdir, f'berfile_{i}')
proc = subprocess.Popen(['/usr/bin/python2', sim,
'--signal-strength', str(ampl_mhz),
'--ber-file', berfile],
stdin=subprocess.PIPE, stdout=subprocess.DEVNULL)
yield proc, sim, ampl_mhz, berfile
results = { sim: defaultdict(lambda: ([], [])) for sim in simulations }
with tqdm(total = nchunks * duration, bar_format='{l_bar}{bar}| {elapsed}<{remaining}') as tq:
tq.write(f'Will launch {len(jobs)} simulation jobs in {nchunks} batches of {forklimit}')
for n, i in enumerate(range(0, len(jobs), forklimit)):
batch = jobs[i:][:forklimit]
tq.write(f'Starting batch {n+1}/{nchunks}...')
processes = list(start_processes(batch))
tq.write('done.')
tq.write('Waiting for simulation:')
for _ in range(100):
time.sleep(duration/100)
tq.update(duration/100)
tq.write('Terminating processes...')
for proc, *_ in processes:
proc.communicate(b'\n', timeout=communicate_timeout)
for proc, *_ in processes:
proc.wait(terminate_timeout)
tq.write('done.')
tq.write('Processing simulation results')
for _proc, sim, ampl_mhz, berfile in processes:
with open(berfile, 'rb') as f:
data = f.read()
floats = struct.unpack(f'{len(data)//4}f', data)
ber = statistics.mean(floats[-256:])
stdev = statistics.stdev(floats[-256:])
bers, stdevs = results[sim][ampl_mhz]
bers.append(ber)
stdevs.append(stdev)
return results
if __name__ == '__main__':
print(run_simulation())