From c5222fd9b45beb783062729b39ded0f503bec4e4 Mon Sep 17 00:00:00 2001 From: Jonathan Hogg Date: Mon, 29 Jun 2020 17:25:43 +0100 Subject: [PATCH] Tidying up a bit; work in progress on automatic analysis of results --- analysis.py | 157 ++++++++++++++++++++++++++++++++++++++++++++++++++++ scope.py | 113 ++++++++++++++++++++----------------- streams.py | 18 +++--- 3 files changed, 229 insertions(+), 59 deletions(-) create mode 100644 analysis.py diff --git a/analysis.py b/analysis.py new file mode 100644 index 0000000..d7c173a --- /dev/null +++ b/analysis.py @@ -0,0 +1,157 @@ + +import numpy as np + + +def interpolate_min_x(f, x): + return 0.5 * (f[x-1] - f[x+1]) / (f[x-1] - 2 * f[x] + f[x+1]) + x + + +def rms(f): + return np.sqrt((f ** 2).mean()) + + +def sine_wave(n): + return np.sin(np.linspace(0, 2*np.pi, n, endpoint=False)) + + +def triangle_wave(n): + x = np.linspace(0, 4, n, endpoint=False) + x2 = x % 2 + y = np.where(x2 < 1, x2, 2 - x2) + y = np.where(x // 2 < 1, y, -y) + return y + + +def square_wave(n, duty=0.5): + w = int(n * duty) + return np.hstack([np.ones(w), -np.ones(n - w)]) + + +def sawtooth_wave(n): + return 2 * (np.linspace(0.5, 1.5, n, endpoint=False) % 1) - 1 + + +def moving_average(samples, width, mode='wrap'): + hwidth = width // 2 + samples = np.take(samples, np.arange(-hwidth, len(samples)+width-hwidth), mode=mode) + cumulative = samples.cumsum() + return (cumulative[width:] - cumulative[:-width]) / width + + +def calculate_periodicity(series, window=0.1): + samples = np.array(series.samples) + window = int(len(samples) * window) + errors = np.zeros(len(samples) - window) + for i in range(1, len(errors) + 1): + errors[i-1] = rms(samples[i:] - samples[:-i]) + threshold = errors.max() / 2 + minima = [] + for i in range(window, len(errors) - window): + p = errors[i-window:i+window].argmin() + if p == window and errors[p + i - window] < threshold: + minima.append(interpolate_min_x(errors, i)) + if len(minima) <= 1: + return None + ks = np.polyfit(np.arange(0, len(minima)), minima, 1) + return ks[0] / series.sample_rate + + +def extract_waveform(series, period): + p = int(round(series.sample_rate * period)) + n = len(series.samples) // p + if n <= 2: + return None, None + samples = np.array(series.samples)[:p*n] + cumsum = samples.cumsum() + underlying = (cumsum[p:] - cumsum[:-p]) / p + n -= 1 + samples = samples[p//2:p*n + p//2] - underlying + wave = np.zeros(p) + for i in range(n): + o = i * p + wave += samples[o:o+p] + wave /= n + return wave, p//2, n, underlying + + +def normalize_waveform(samples, smooth=7): + n = len(samples) + smoothed = moving_average(samples, smooth) + scale = (smoothed.max() - smoothed.min()) / 2 + offset = (smoothed.max() + smoothed.min()) / 2 + smoothed -= offset + last_rising = first_falling = None + crossings = [] + for i in range(n): + if smoothed[i-1] < 0 and smoothed[i] > 0: + last_rising = i + elif smoothed[i-1] > 0 and smoothed[i] < 0: + if last_rising is None: + first_falling = i + else: + crossings.append((i - last_rising, last_rising)) + if first_falling is not None: + crossings.append((n + first_falling - last_rising, last_rising)) + width, first = min(crossings) + wave = np.hstack([smoothed[first:], smoothed[:first]]) / scale + return wave, offset, scale, first, sorted((i - first % n, w) for (w, i) in crossings) + + +def characterize_waveform(samples, crossings): + n = len(samples) + possibles = [] + if len(crossings) == 1: + duty_cycle = crossings[0][1] / n + if duty_cycle > 0.45 and duty_cycle < 0.55: + possibles.append((rms(samples - sine_wave(n)), 'sine', None)) + possibles.append((rms(samples - triangle_wave(n)), 'triangle', None)) + possibles.append((rms(samples - sawtooth_wave(n)), 'sawtooth', None)) + possibles.append((rms(samples - square_wave(n, duty_cycle)), 'square', duty_cycle)) + possibles.sort() + return possibles + + +def analyze_series(series): + period = calculate_periodicity(series) + if period is not None: + waveform = DotDict(period=period, frequency=1 / period) + wave, start, count, underlying = extract_waveform(series, period) + wave, offset, scale, first, crossings = normalize_waveform(wave) + waveform.samples = wave + waveform.beginning = start + first + waveform.count = count + waveform.amplitude = scale + waveform.offset = underlying.mean() + offset + possibles = characterize_waveform(wave, crossings) + if possibles: + waveform.error, waveform.shape, waveform.duty_cycle = possibles[0] + series.waveform = waveform + + +# %% + +from pylab import plot +from utils import DotDict + +o = 400 +m = 5 +n = o * m +samples = sine_wave(o) +samples = np.hstack([samples] * m) * 2 + 5 +samples = np.hstack([samples[100:], samples[:100]]) +samples += np.random.normal(size=n) * 0.1 +samples += np.linspace(-0.5, 0.5, n) +series = DotDict(samples=samples, sample_rate=1000000) + +analyze_series(series) + +print(series.waveform.frequency) +print(series.waveform.shape) +print(series.waveform.amplitude, series.waveform.offset) + +plot(series.samples) + +wave = np.hstack([series.waveform.samples[-series.waveform.beginning:]] + + [series.waveform.samples] * series.waveform.count + + [series.waveform.samples[:-series.waveform.beginning]]) +plot(wave * series.waveform.amplitude + series.waveform.offset) diff --git a/scope.py b/scope.py index 9561dfc..3ab7a3a 100755 --- a/scope.py +++ b/scope.py @@ -16,9 +16,8 @@ from utils import DotDict import vm -LOG = logging.getLogger(__name__) - -ANALOG_PARAMETERS_PATH = Path('~/.config/scopething/analog.conf').expanduser() +Log = logging.getLogger(__name__) +AnalogParametersPath = Path('~/.config/scopething/analog.conf').expanduser() class UsageError(Exception): @@ -44,7 +43,7 @@ class Scope(vm.VirtualMachine): break else: raise RuntimeError("No matching serial device found") - LOG.info(f"Connecting to scope at {url}") + Log.info(f"Connecting to scope at {url}") self.close() parts = urlparse(url, scheme='file') if parts.scheme == 'file': @@ -59,7 +58,7 @@ class Scope(vm.VirtualMachine): return self async def reset(self): - LOG.info("Resetting scope") + Log.info("Resetting scope") await self.issue_reset() await self.issue_get_revision() revision = ((await self.read_replies(2))[1]).decode('ascii') @@ -82,31 +81,31 @@ class Scope(vm.VirtualMachine): self._awg_running = False self._clock_running = False self.load_analog_params() - LOG.info(f"Initialised scope, revision: {revision}") + Log.info(f"Initialised scope, revision: {revision}") def load_analog_params(self): config = ConfigParser() - config.read(ANALOG_PARAMETERS_PATH) + config.read(AnalogParametersPath) analog_params = {} for url in config.sections(): if url == self.url: for probes in config[url]: params = self.AnalogParams(*map(float, config[url][probes].split())) analog_params[probes] = params - LOG.debug(f"Loading saved parameters for {probes}: {params!r}") + Log.debug(f"Loading saved parameters for {probes}: {params!r}") if analog_params: self.analog_params.update(analog_params) - LOG.info(f"Loaded analog parameters for probes: {', '.join(analog_params.keys())}") + Log.info(f"Loaded analog parameters for probes: {', '.join(analog_params.keys())}") def save_analog_params(self): - LOG.info("Saving analog parameters") + Log.info("Saving analog parameters") config = ConfigParser() - config.read(ANALOG_PARAMETERS_PATH) + config.read(AnalogParametersPath) config[self.url] = {probes: ' '.join(map(str, self.analog_params[probes])) for probes in self.analog_params} - parent = ANALOG_PARAMETERS_PATH.parent + parent = AnalogParametersPath.parent if not parent.is_dir(): parent.mkdir(parents=True) - with open(ANALOG_PARAMETERS_PATH, 'w') as parameters_file: + with open(AnalogParametersPath, 'w') as parameters_file: config.write(parameters_file) def __enter__(self): @@ -117,7 +116,7 @@ class Scope(vm.VirtualMachine): def close(self): super().close() - LOG.info("Closed scope") + Log.info("Closed scope") def calculate_lo_hi(self, low, high, params): if not isinstance(params, self.AnalogParams): @@ -164,20 +163,20 @@ class Scope(vm.VirtualMachine): ticks = int(round(period / self.master_clock_period / nsamples)) clock_scale = 1 if capture_mode.analog_channels == len(analog_channels) and capture_mode.logic_channels == bool(logic_channels): - LOG.debug(f"Considering trace mode {capture_mode.trace_mode.name}...") + Log.debug(f"Considering trace mode {capture_mode.trace_mode.name}...") if ticks > capture_mode.clock_high and capture_mode.clock_divide > 1: clock_scale = int(math.ceil(period / self.master_clock_period / nsamples / capture_mode.clock_high)) ticks = int(round(period / self.master_clock_period / nsamples / clock_scale)) if ticks in range(capture_mode.clock_low, capture_mode.clock_high+1): - LOG.debug(f"- try with tick count {ticks} x {clock_scale}") + Log.debug(f"- try with tick count {ticks} x {clock_scale}") else: continue elif ticks >= capture_mode.clock_low: if ticks > capture_mode.clock_high: ticks = capture_mode.clock_high - LOG.debug(f"- try with tick count {ticks}") + Log.debug(f"- try with tick count {ticks}") else: - LOG.debug("- mode too slow") + Log.debug("- mode too slow") continue n = int(round(period / self.master_clock_period / ticks / clock_scale)) if len(analog_channels) == 2: @@ -186,16 +185,16 @@ class Scope(vm.VirtualMachine): if logic_channels and analog_channels: buffer_width //= 2 if n <= buffer_width: - LOG.debug(f"- OK; period is {n} samples") + Log.debug(f"- OK; period is {n} samples") nsamples = n break - LOG.debug(f"- insufficient buffer space for necessary {n} samples") + Log.debug(f"- insufficient buffer space for necessary {n} samples") else: raise ConfigurationError("Unable to find appropriate capture mode") sample_period = ticks*clock_scale*self.master_clock_period sample_rate = 1/sample_period if trigger_position and sample_rate > 5e6: - LOG.warn("Pre-trigger capture not supported above 5M samples/s; forcing trigger_position=0") + Log.warn("Pre-trigger capture not supported above 5M samples/s; forcing trigger_position=0") trigger_position = 0 if raw: @@ -206,11 +205,11 @@ class Scope(vm.VirtualMachine): if low is None: low = analog_params.safe_low if analog_channels else self.logic_low elif low < analog_params.safe_low: - LOG.warning(f"Voltage range is below safe minimum: {low} < {analog_params.safe_low}") + Log.warning(f"Voltage range is below safe minimum: {low} < {analog_params.safe_low}") if high is None: high = analog_params.safe_high if analog_channels else self.logic_high elif high > analog_params.safe_high: - LOG.warning(f"Voltage range is above safe maximum: {high} > {analog_params.safe_high}") + Log.warning(f"Voltage range is above safe maximum: {high} > {analog_params.safe_high}") lo, hi = self.calculate_lo_hi(low, high, analog_params) spock_option = vm.SpockOption.TriggerTypeHardwareComparator @@ -225,8 +224,7 @@ class Scope(vm.VirtualMachine): kitchen_sink_b |= vm.KitchenSinkB.AnalogFilterEnable if trigger_level is None: trigger_level = (high + low) / 2 - if not raw: - trigger_level = (trigger_level - analog_params.offset) / analog_params.scale + analog_trigger_level = (trigger_level - analog_params.offset) / analog_params.scale if not raw else trigger_level if trigger == 'A' or trigger == 'B': if trigger == 'A': spock_option |= vm.SpockOption.TriggerSourceA @@ -273,12 +271,12 @@ class Scope(vm.VirtualMachine): else: raise ConfigurationError("Required trigger timeout too long, use a later trigger position") - LOG.info(f"Begin {('mixed' if logic_channels else 'analogue') if analog_channels else 'logic'} signal capture " + Log.info(f"Begin {('mixed' if logic_channels else 'analogue') if analog_channels else 'logic'} signal capture " f"at {sample_rate:,.0f} samples per second (trace mode {capture_mode.trace_mode.name})") async with self.transaction(): await self.set_registers(TraceMode=capture_mode.trace_mode, BufferMode=capture_mode.buffer_mode, SampleAddress=0, ClockTicks=ticks, ClockScale=clock_scale, - TriggerLevel=trigger_level, TriggerLogic=trigger_logic, TriggerMask=trigger_mask, + TriggerLevel=analog_trigger_level, TriggerLogic=trigger_logic, TriggerMask=trigger_mask, TraceIntro=trace_intro, TraceOutro=trace_outro, TraceDelay=0, Timeout=trigger_timeout, TriggerIntro=trigger_intro//2, TriggerOutro=trigger_outro//2, Prelude=0, SpockOption=spock_option, ConverterLo=lo, ConverterHi=hi, @@ -304,7 +302,7 @@ class Scope(vm.VirtualMachine): address -= address % 2 traces = DotDict() - timestamps = array.array('d', (t*self.master_clock_period for t in range(start_timestamp, timestamp, ticks*clock_scale))) + timestamps = array.array('d', (t*self.master_clock_period for t in range(0, timestamp, ticks*clock_scale))) for dump_channel, channel in enumerate(sorted(analog_channels)): asamples = nsamples // len(analog_channels) async with self.transaction(): @@ -315,11 +313,18 @@ class Scope(vm.VirtualMachine): await self.issue_analog_dump_binary() value_multiplier, value_offset = (1, 0) if raw else (high-low, low-analog_params.ab_offset/2*(1 if channel == 'A' else -1)) data = await self.read_analog_samples(asamples, capture_mode.sample_width) - traces[channel] = DotDict({'timestamps': timestamps[dump_channel::len(analog_channels)] if len(analog_channels) > 1 else timestamps, - 'samples': array.array('f', (value*value_multiplier+value_offset for value in data)), - 'sample_period': sample_period*len(analog_channels), - 'sample_rate': sample_rate/len(analog_channels), - 'cause': cause}) + series = DotDict({'channel': channel, + 'start_timestamp': start_timestamp, + 'timestamps': timestamps[dump_channel::len(analog_channels)] if len(analog_channels) > 1 else timestamps, + 'samples': array.array('f', (value*value_multiplier+value_offset for value in data)), + 'sample_period': sample_period*len(analog_channels), + 'sample_rate': sample_rate/len(analog_channels), + 'cause': cause}) + if cause == 'trigger' and channel == trigger: + series.trigger_timestamp = series.timestamps[trigger_samples // len(analog_channels)] + series.trigger_level = trigger_level + series.trigger_type = trigger_type + traces[channel] = series if logic_channels: async with self.transaction(): await self.set_registers(SampleAddress=(address - nsamples) % buffer_width, @@ -329,12 +334,20 @@ class Scope(vm.VirtualMachine): data = await self.read_logic_samples(nsamples) for i in logic_channels: mask = 1 << i - traces[f'L{i}'] = DotDict({'timestamps': timestamps, - 'samples': array.array('B', (1 if value & mask else 0 for value in data)), - 'sample_period': sample_period, - 'sample_rate': sample_rate, - 'cause': cause}) - LOG.info(f"{nsamples} samples captured on {cause}, traces: {', '.join(traces)}") + channel = f'L{i}' + series = DotDict({'channel': channel, + 'start_timestamp': start_timestamp, + 'timestamps': timestamps, + 'samples': array.array('B', (1 if value & mask else 0 for value in data)), + 'sample_period': sample_period, + 'sample_rate': sample_rate, + 'cause': cause}) + if cause == 'trigger' and isinstance(trigger, dict) and i in trigger: + series.trigger_timestamp = series.timestamps[trigger_samples] + series.trigger_level = trigger[i] + series.trigger_type = trigger_type + traces[channel] = series + Log.info(f"{nsamples} samples captured on {cause}, traces: {', '.join(traces)}") return traces async def start_waveform(self, frequency, waveform='sine', ratio=0.5, low=0, high=None, min_samples=50, max_error=1e-4): @@ -355,7 +368,7 @@ class Scope(vm.VirtualMachine): size = int(round(nwaves * width)) actualf = self.master_clock_rate * nwaves / size / clock if actualf == frequency: - LOG.debug(f"Exact solution: size={size} nwaves={nwaves} clock={clock}") + Log.debug(f"Exact solution: size={size} nwaves={nwaves} clock={clock}") break error = abs(frequency - actualf) / frequency if error < max_error and (best_solution is None or error < best_solution[0]): @@ -364,7 +377,7 @@ class Scope(vm.VirtualMachine): if best_solution is None: raise ConfigurationError("No solution to required frequency/min_samples/max_error") error, size, nwaves, clock, actualf = best_solution - LOG.debug(f"Best solution: size={size} nwaves={nwaves} clock={clock} actualf={actualf}") + Log.debug(f"Best solution: size={size} nwaves={nwaves} clock={clock} actualf={actualf}") async with self.transaction(): if isinstance(waveform, str): mode = {'sine': 0, 'triangle': 1, 'exponential': 2, 'square': 3}[waveform.lower()] @@ -391,7 +404,7 @@ class Scope(vm.VirtualMachine): await self.set_registers(KitchenSinkB=vm.KitchenSinkB.WaveformGeneratorEnable) await self.issue_configure_device_hardware() self._awg_running = True - LOG.info(f"Signal generator running at {actualf:0.1f}Hz") + Log.info(f"Signal generator running at {actualf:0.1f}Hz") return actualf async def stop_waveform(self): @@ -402,7 +415,7 @@ class Scope(vm.VirtualMachine): await self.issue_control_clock_generator() await self.set_registers(KitchenSinkB=0) await self.issue_configure_device_hardware() - LOG.info("Signal generator stopped") + Log.info("Signal generator stopped") self._awg_running = False async def start_clock(self, frequency, ratio=0.5, max_error=1e-4): @@ -417,7 +430,7 @@ class Scope(vm.VirtualMachine): await self.set_registers(Map5=0x12, Clock=ticks, Rise=0, Fall=fall, Control=0x80, Cmd=3, Mode=0) await self.issue_control_clock_generator() self._clock_running = True - LOG.info(f"Clock generator running at {actualf:0.1f}Hz, {actualr*100:.0f}% duty cycle") + Log.info(f"Clock generator running at {actualf:0.1f}Hz, {actualr*100:.0f}% duty cycle") return actualf, actualr async def stop_clock(self): @@ -426,7 +439,7 @@ class Scope(vm.VirtualMachine): async with self.transaction(): await self.set_registers(Map5=0, Cmd=1, Mode=0) await self.issue_control_clock_generator() - LOG.info("Clock generator stopped") + Log.info("Clock generator stopped") self._clock_running = False async def calibrate(self, probes='x1', n=32, save=True): @@ -477,7 +490,7 @@ class Scope(vm.VirtualMachine): full = (full + 1) / 3 analog_scale = self.clock_voltage / (full - zero) analog_offset = -zero * analog_scale - LOG.info(f"Analog full range = {analog_scale:.2f}V, zero offset = {analog_offset:.2f}V") + Log.info(f"Analog full range = {analog_scale:.2f}V, zero offset = {analog_offset:.2f}V") for lo in np.linspace(self.analog_lo_min, 0.5, n, endpoint=False): for hi in np.linspace(self.analog_hi_max, 0.5, n): zero, full, offset = await measure(lo, hi, 2e-3 if len(items) % 4 < 2 else 1e-3, len(items) % 2 == 0) @@ -497,7 +510,7 @@ class Scope(vm.VirtualMachine): constraints=[{'type': 'eq', 'fun': lambda x: x[0]*1/3 + x[1]*2/3 + x[2] - 1/3}, {'type': 'eq', 'fun': lambda x: x[3]*2/3 + x[4]*1/3 + x[5] - 2/3}]) if result.success: - LOG.info(f"Calibration succeeded: {result.message}") + Log.info(f"Calibration succeeded: {result.message}") params = self.AnalogParams(*result.x, analog_scale, analog_offset, None, None, None) def f(x): @@ -507,15 +520,15 @@ class Scope(vm.VirtualMachine): safe_low, safe_high = minimize(f, (low[0], high[0])).x offset_mean = offset.mean() params = self.analog_params[probes] = self.AnalogParams(*result.x, analog_scale, analog_offset, safe_low, safe_high, offset_mean) - LOG.info(f"{params!r} ±{100*offset.std()/offset_mean:.1f}%)") + Log.info(f"{params!r} ±{100*offset.std()/offset_mean:.1f}%)") clo, chi = self.calculate_lo_hi(low, high, params) lo_error = np.sqrt((((clo-lo)/(hi-lo))**2).mean()) hi_error = np.sqrt((((chi-hi)/(hi-lo))**2).mean()) - LOG.info(f"Mean error: lo={lo_error*10000:.1f}bps hi={hi_error*10000:.1f}bps") + Log.info(f"Mean error: lo={lo_error*10000:.1f}bps hi={hi_error*10000:.1f}bps") if save: self.save_analog_params() else: - LOG.warning(f"Calibration failed: {result.message}") + Log.warning(f"Calibration failed: {result.message}") return result.success def __repr__(self): diff --git a/streams.py b/streams.py index fbd999a..5f13b8f 100644 --- a/streams.py +++ b/streams.py @@ -14,7 +14,7 @@ import serial from serial.tools.list_ports import comports -LOG = logging.getLogger(__name__) +Log = logging.getLogger(__name__) class SerialStream: @@ -36,7 +36,7 @@ class SerialStream: self._use_threads = sys.platform == 'win32' if use_threads is None else use_threads self._connection = serial.Serial(self._device, **kwargs) if self._use_threads else \ serial.Serial(self._device, timeout=0, write_timeout=0, **kwargs) - LOG.debug(f"Opened SerialStream on {device}") + Log.debug(f"Opened SerialStream on {device}") self._loop = loop if loop is not None else asyncio.get_event_loop() self._output_buffer = bytes() self._output_buffer_empty = None @@ -63,10 +63,10 @@ class SerialStream: except serial.SerialTimeoutException: n = 0 except Exception: - LOG.exception("Error writing to stream") + Log.exception("Error writing to stream") raise if n: - LOG.debug(f"Write {data[:n]!r}") + Log.debug(f"Write {data[:n]!r}") self._output_buffer = data[n:] else: self._output_buffer += data @@ -84,11 +84,11 @@ class SerialStream: except serial.SerialTimeoutException: n = 0 except Exception as e: - LOG.exception("Error writing to stream") + Log.exception("Error writing to stream") self._output_buffer_empty.set_exception(e) self._loop.remove_writer(self._connection) if n: - LOG.debug(f"Write {self._output_buffer[:n]!r}") + Log.debug(f"Write {self._output_buffer[:n]!r}") self._output_buffer = self._output_buffer[n:] if not self._output_buffer: self._loop.remove_writer(self._connection) @@ -104,7 +104,7 @@ class SerialStream: n = self._connection.write(data) finally: self._output_buffer_lock.acquire() - LOG.debug(f"Write {self._output_buffer[:n]!r}") + Log.debug(f"Write {self._output_buffer[:n]!r}") self._output_buffer = self._output_buffer[n:] self._output_buffer_empty = None @@ -115,7 +115,7 @@ class SerialStream: w = self._connection.in_waiting if w: data = self._connection.read(w if n is None else min(n, w)) - LOG.debug(f"Read {data!r}") + Log.debug(f"Read {data!r}") return data else: future = self._loop.create_future() @@ -130,7 +130,7 @@ class SerialStream: w = self._connection.in_waiting if w and (n is None or n > 1): data += self._connection.read(w if n is None else min(n-1, w)) - LOG.debug(f"Read {data!r}") + Log.debug(f"Read {data!r}") return data async def readexactly(self, n):