diff options
Diffstat (limited to 'bin/mimosa-etv')
| -rwxr-xr-x | bin/mimosa-etv | 200 |
1 files changed, 200 insertions, 0 deletions
diff --git a/bin/mimosa-etv b/bin/mimosa-etv new file mode 100755 index 0000000..9d98592 --- /dev/null +++ b/bin/mimosa-etv @@ -0,0 +1,200 @@ +#!/usr/bin/env python3 +# vim:tabstop=4:softtabstop=4:shiftwidth=4:textwidth=160:smarttab:expandtab + +import getopt +import itertools +import matplotlib.pyplot as plt +import numpy as np +import os +import re +import sys +from dfatool import aggregate_measures, running_mean, MIMOSA + +opt = dict() + +def show_help(): + print('''mimosa-etv - MIMOSA Analyzer and Visualizer + +USAGE + +mimosa-etv [--skip <count>] [--threshold <power>] [--plot] [--stat] <file> + +DESCRIPTION + +mimosa-etv analyzes measurements taken via MIMOSA. Data can be plotted or aggregated on stdout. + +OPTIONS + + --skip <count> + Skip the first <count> data samples. + --threshold <watts>|mean + Partition data into points with mean power >= <watts> and points with + mean power < <watts>, and print some statistics. higher power is handled + as peaks, whereas low-power measurements constitute the baseline. + If the threshold is set to "mean", the mean power of all measurements + will be used + --threshold-peakcount <num> + Automatically determine threshold so that there are exactly <num> peaks. + A peaks is a group of consecutive measurements with mean power >= threshold + --plot + Show power/time plot + --stat + Show mean voltage, current, and power as well as total energy consumption. + ''') + +def peak_search(data, lower, upper, direction_function): + while upper - lower > 1e-6: + bs_test = np.mean([lower, upper]) + peakcount = itertools.groupby(data, lambda x: x >= bs_test) + peakcount = filter(lambda x: x[0] == True, peakcount) + peakcount = sum(1 for i in peakcount) + direction = direction_function(peakcount, bs_test) + if direction == 0: + return bs_test + elif direction == 1: + lower = bs_test + else: + upper = bs_test + return None + +def peak_search2(data, lower, upper, check_function): + for power in np.arange(lower, upper, 1e-6): + peakcount = itertools.groupby(data, lambda x: x >= power) + peakcount = filter(lambda x: x[0] == True, peakcount) + peakcount = sum(1 for i in peakcount) + if check_function(peakcount, power) == 0: + return power + return None + +if __name__ == '__main__': + try: + optspec = ('help skip= threshold= threshold-peakcount= plot stat') + raw_opts, args = getopt.getopt(sys.argv[1:], "", optspec.split(' ')) + + for option, parameter in raw_opts: + optname = re.sub(r'^--', '', option) + opt[optname] = parameter + + if 'help' in opt: + show_help() + sys.exit(0) + + if 'skip' in opt: + opt['skip'] = int(opt['skip']) + else: + opt['skip'] = 0 + + if 'threshold' in opt and opt['threshold'] != 'mean': + opt['threshold'] = float(opt['threshold']) + + if 'threshold-peakcount' in opt: + opt['threshold-peakcount'] = int(opt['threshold-peakcount']) + + except getopt.GetoptError as err: + print(err) + sys.exit(2) + except IndexError: + print('Usage: mimosa-etv <duration>') + sys.exit(2) + except ValueError: + print('Error: duration or skip is not a number') + sys.exit(2) + + voltage, shunt, inputfile = args + voltage = float(voltage) + shunt = int(shunt) + mim = MIMOSA(voltage, shunt) + charges, triggers = mim.load_file(inputfile) + + currents = mim.charge_to_current_nocal(charges) * 1e-6 + powers = currents * voltage + + if 'threshold-peakcount' in opt: + bs_mean = np.mean(powers) + + # Finding the correct threshold is tricky. If #peaks < peakcont, our + # current threshold may be too low (extreme case: a single peak + # containing all measurements), but it may also be too high (extreme + # case: a single peak containing just one data point). Similarly, + # #peaks > peakcount may be due to baseline noise causing lots of + # small peaks, or due to peak noise (if the threshold is already rather + # high). + # For now, we first try a simple binary search: + # The threshold is probably somewhere around the mean, so if + # #peaks != peakcount and threshold < mean, we go up, and if + # #peaks != peakcount and threshold >= mean, we go down. + # If that doesn't work, we fall back to a linear search in 1 µW steps + def direction_function(peakcount, power): + if peakcount == opt['threshold-peakcount']: + return 0 + if power < bs_mean: + return 1 + return -1 + threshold = peak_search(power, np.min(power), np.max(power), direction_function) + if threshold == None: + threshold = peak_search2(power, np.min(power), np.max(power), direction_function) + + if threshold != None: + print('Threshold set to {:.0f} µW : {:.9f}'.format(threshold * 1e6, threshold)) + opt['threshold'] = threshold + else: + print('Found no working threshold') + + if 'threshold' in opt: + if opt['threshold'] == 'mean': + opt['threshold'] = np.mean(powers) + print('Threshold set to {:.0f} µW : {:.9f}'.format(opt['threshold'] * 1e6, opt['threshold'])) + + baseline_mean = 0 + if np.any(powers < opt['threshold']): + baseline_mean = np.mean(powers[powers < opt['threshold']]) + print('Baseline mean: {:.0f} µW : {:.9f}'.format( + baseline_mean * 1e6, baseline_mean)) + if np.any(powers >= opt['threshold']): + print('Peak mean: {:.0f} µW : {:.9f}'.format( + np.mean(powers[powers >= opt['threshold']]) * 1e6, + np.mean(powers[powers >= opt['threshold']]))) + + peaks = [] + peak_start = -1 + for i, dp in enumerate(powers): + if dp >= opt['threshold'] and peak_start == -1: + peak_start = i + elif dp < opt['threshold'] and peak_start != -1: + peaks.append((peak_start, i)) + peak_start = -1 + + total_energy = 0 + delta_energy = 0 + for peak in peaks: + duration = (peak[1] - peak[0]) * 1e-5 + total_energy += np.mean(powers[peak[0] : peak[1]]) * duration + delta_energy += (np.mean(powers[peak[0] : peak[1]]) - baseline_mean) * duration + delta_powers = powers[peak[0] : peak[1]] - baseline_mean + print('{:.2f}ms peak ({:f} -> {:f})'.format(duration * 1000, + peak[0], peak[1])) + print(' {:f} µJ / mean {:f} µW'.format( + np.mean(powers[peak[0] : peak[1]]) * duration * 1e6, + np.mean(powers[peak[0] : peak[1]]) * 1e6 )) + measures = aggregate_measures(np.mean(delta_powers), delta_powers) + print(' {:f} µW delta mean = {:0.1f}% / {:f} µW error'.format(np.mean(delta_powers) * 1e6, measures['smape'], measures['rmsd'] * 1e6 )) + print('Peak energy mean: {:.0f} µJ : {:.9f}'.format( + total_energy * 1e6 / len(peaks), total_energy / len(peaks))) + print('Average per-peak energy (delta over baseline): {:.0f} µJ : {:.9f}'.format( + delta_energy * 1e6 / len(peaks), delta_energy / len(peaks))) + + + if 'stat' in opt: + mean_current = np.mean(currents) + mean_power = np.mean(powers) + print('Mean current: {:.0f} µA : {:.9f}'.format(mean_current * 1e6, mean_current)) + print('Mean power: {:.0f} µW : {:.9f}'.format(mean_power * 1e6, mean_power)) + + if 'plot' in opt: + timestamps = np.arange(len(powers)) * 1e-5 + pwrhandle, = plt.plot(timestamps, powers, 'b-', label='U*I', markersize=1) + plt.legend(handles=[pwrhandle]) + plt.xlabel('Time [s]') + plt.ylabel('Power [W]') + plt.grid(True) + plt.show() |