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

Fix frequency measurement simulation

Browse source
This commit is contained in:
jaseg 2020-03-04 17:10:31 +01:00
parent 4d80111cad
commit d9b26d16c0
22 changed files with 1201 additions and 164 deletions
Download patch file Download diff file Expand all files Collapse all files

29
controller/fw/Makefile
View file

@ -46,17 +46,17 @@ FMEAS_ADC_SAMPLING_RATE ?= 1000
FMEAS_ADC_MAX ?= 4096
FMEAS_FFT_LEN ?= 256
FMEAS_FFT_WINDOW ?= gaussian
FMEAS_FFT_WINDOW_SIGMA ?= 8.0
FMEAS_FFT_WINDOW_SIGMA ?= 16.0
CC ?= $(PREFIX)gcc
CXX ?= $(PREFIX)g++
LD ?= $(PREFIX)gcc
AR ?= $(PREFIX)ar
AS ?= $(PREFIX)as
OBJCOPY ?= $(PREFIX)objcopy
OBJDUMP ?= $(PREFIX)objdump
GDB ?= $(PREFIX)gdb
CC := $(PREFIX)gcc
CXX := $(PREFIX)g++
LD := $(PREFIX)gcc
AR := $(PREFIX)ar
AS := $(PREFIX)as
OBJCOPY := $(PREFIX)objcopy
OBJDUMP := $(PREFIX)objdump
GDB := $(PREFIX)gdb
HOST_CC ?= $(HOST_PREFIX)gcc
HOST_CXX ?= $(HOST_PREFIX)g++
@ -75,9 +75,9 @@ LIBSODIUM_DIR_ABS := $(abspath $(LIBSODIUM_DIR))
TINYAES_DIR_ABS := $(abspath $(TINYAES_DIR))
MUSL_DIR_ABS := $(abspath $(MUSL_DIR))
COMMON_CFLAGS += -I$(OPENCM3_DIR_ABS)/include -Imspdebug/util -Imspdebug/drivers
COMMON_CFLAGS += -I$(OPENCM3_DIR_ABS)/include -Imspdebug/util -Imspdebug/drivers -Ilevmarq
COMMON_CFLAGS += -I$(CMSIS_DIR_ABS)/CMSIS/DSP/Include -I$(CMSIS_DIR_ABS)/CMSIS/Core/Include
COMMON_CFLAGS += -I$(abspath musl_include_shims) -Ilevmarq
CFLAGS += -I$(abspath musl_include_shims)
COMMON_CFLAGS += -Os -std=gnu11 -g -DSTM32F4
CFLAGS += -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16
@ -91,7 +91,7 @@ COMMON_CFLAGS += -DFMEAS_FFT_WINDOW_SIGMA=$(FMEAS_FFT_WINDOW_SIGMA)
# for musl
CFLAGS += -Dhidden=
SIM_CFLAGS += -Isrc
SIM_CFLAGS += -Isrc -lm -DSIMULATION
INT_CFLAGS += -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef
INT_CFLAGS += -Wredundant-decls -Wmissing-prototypes -Wstrict-prototypes
@ -112,7 +112,7 @@ LDFLAGS += -L$(OPENCM3_DIR_ABS)/lib -l$(OPENCM3_LIB) $(shell $(CC) -print-libg
all: $(BUILDDIR)/$(BINARY)
tests: $(BUILDDIR)/tools/freq_meas_test.elf
tests: $(BUILDDIR)/tools/freq_meas_test
OBJS := $(addprefix $(BUILDDIR)/,$(C_SOURCES:.c=.o) $(CXX_SOURCES:.cpp=.o))
@ -127,7 +127,7 @@ ALL_OBJS += $(BUILDDIR)/generated/fmeas_fft_window.o
$(BUILDDIR)/$(BINARY): $(ALL_OBJS)
$(LD) -T$(LDSCRIPT) $(COMMON_LDFLAGS) $(LDFLAGS) -o $@ -Wl,-Map=$(BUILDDIR)/src/$*.map $^
$(BUILDDIR)/tools/freq_meas_test.elf: tools/freq_meas_test.c src/freq_meas.c levmarq/levmarq.c $(BUILDDIR)/generated/fmeas_fft_window.c $(CMSIS_SOURCES)
$(BUILDDIR)/tools/freq_meas_test: tools/freq_meas_test.c src/freq_meas.c levmarq/levmarq.c $(BUILDDIR)/generated/fmeas_fft_window.c $(CMSIS_SOURCES)
mkdir -p $(@D)
$(HOST_CC) $(COMMON_CFLAGS) $(SIM_CFLAGS) -o $@ $^
@ -177,6 +177,7 @@ clean:
-rm -r $(BUILDDIR)/src
-rm -r $(BUILDDIR)/generated
-rm $(BUILDDIR)/$(BINARY)
-rm $(BUILDDIR)/tools/freq_meas_test
mrproper: clean
-rm -r build

250
controller/fw/levmarq/levmarq.c
View file

@ -20,6 +20,7 @@
#include <arm_math.h>
#include "levmarq.h"
#include "simulation.h"
#define TOL 1e-20f /* smallest value allowed in cholesky_decomp() */
@ -28,18 +29,20 @@
/* set parameters required by levmarq() to default values */
void levmarq_init(LMstat *lmstat)
{
lmstat->max_it = 10000;
lmstat->init_lambda = 0.0001f;
lmstat->up_factor = 10.0f;
lmstat->down_factor = 10.0f;
lmstat->target_derr = 1e-12f;
lmstat->max_it = 10000;
lmstat->init_lambda = 0.0001f;
lmstat->up_factor = 10.0f;
lmstat->down_factor = 10.0f;
lmstat->target_derr = 1e-12f;
}
#ifndef SIMULATION
float sqrtf(float arg) {
float out=NAN;
arm_sqrt_f32(arg, &out);
return out;
}
#endif
/* perform least-squares minimization using the Levenberg-Marquardt
algorithm. The arguments are as follows:
@ -49,140 +52,147 @@ float sqrtf(float arg) {
ny number of measurements to be fit
y array of measurements
dysq array of error in measurements, squared
(set dysq=NULL for unweighted least-squares)
(set dysq=NULL for unweighted least-squares)
func function to be fit
grad gradient of "func" with respect to the input parameters
fdata pointer to any additional data required by the function
lmstat pointer to the "status" structure, where minimization parameters
are set and the final status is returned.
are set and the final status is returned.
Before calling levmarq, several of the parameters in lmstat must be set.
For default values, call levmarq_init(lmstat).
*/
*/
int levmarq(int npar, float *par, int ny, float *y, float *dysq,
float (*func)(float *, int, void *),
void (*grad)(float *, float *, int, void *),
void *fdata, LMstat *lmstat)
float (*func)(float *, int, void *),
void (*grad)(float *, float *, int, void *),
void *fdata, LMstat *lmstat)
{
int x,i,j,it,nit,ill;
float lambda,up,down,mult,weight,err,newerr,derr,target_derr;
float h[npar][npar],ch[npar][npar];
float g[npar],d[npar],delta[npar],newpar[npar];
int x,i,j,it,nit,ill;
float lambda,up,down,mult,weight,err,newerr,derr,target_derr;
float h[npar][npar],ch[npar][npar];
float g[npar],d[npar],delta[npar],newpar[npar];
nit = lmstat->max_it;
lambda = lmstat->init_lambda;
up = lmstat->up_factor;
down = 1/lmstat->down_factor;
target_derr = lmstat->target_derr;
weight = 1;
derr = newerr = 0; /* to avoid compiler warnings */
nit = lmstat->max_it;
lambda = lmstat->init_lambda;
up = lmstat->up_factor;
down = 1/lmstat->down_factor;
target_derr = lmstat->target_derr;
weight = 1;
derr = newerr = 0; /* to avoid compiler warnings */
/* calculate the initial error ("chi-squared") */
err = error_func(par, ny, y, dysq, func, fdata);
/* calculate the initial error ("chi-squared") */
err = error_func(par, ny, y, dysq, func, fdata);
/* main iteration */
for (it=0; it<nit; it++) {
/* main iteration */
for (it=0; it<nit; it++) {
//DEBUG_PRINT("iteration %d", it);
/* calculate the approximation to the Hessian and the "derivative" d */
for (i=0; i<npar; i++) {
d[i] = 0;
for (j=0; j<=i; j++)
h[i][j] = 0;
}
for (x=0; x<ny; x++) {
if (dysq) weight = 1/dysq[x]; /* for weighted least-squares */
grad(g, par, x, fdata);
for (i=0; i<npar; i++) {
d[i] += (y[x] - func(par, x, fdata))*g[i]*weight;
for (j=0; j<=i; j++)
h[i][j] += g[i]*g[j]*weight;
}
/* calculate the approximation to the Hessian and the "derivative" d */
for (i=0; i<npar; i++) {
d[i] = 0;
for (j=0; j<=i; j++)
h[i][j] = 0;
}
for (x=0; x<ny; x++) {
if (dysq) weight = 1/dysq[x]; /* for weighted least-squares */
grad(g, par, x, fdata);
for (i=0; i<npar; i++) {
d[i] += (y[x] - func(par, x, fdata))*g[i]*weight;
for (j=0; j<=i; j++)
h[i][j] += g[i]*g[j]*weight;
}
}
/* make a step "delta." If the step is rejected, increase
lambda and try again */
mult = 1 + lambda;
ill = 1; /* ill-conditioned? */
while (ill && (it<nit)) {
for (i=0; i<npar; i++)
h[i][i] = h[i][i]*mult;
ill = cholesky_decomp(npar, ch, h);
if (!ill) {
solve_axb_cholesky(npar, ch, delta, d);
for (i=0; i<npar; i++)
newpar[i] = par[i] + delta[i];
newerr = error_func(newpar, ny, y, dysq, func, fdata);
derr = newerr - err;
ill = (derr > 0);
}
if (ill) {
mult = (1 + lambda*up)/(1 + lambda);
lambda *= up;
it++;
}
}
for (i=0; i<npar; i++)
par[i] = newpar[i];
err = newerr;
lambda *= down;
if ((!ill)&&(-derr<target_derr)) break;
}
/* make a step "delta." If the step is rejected, increase
lambda and try again */
mult = 1 + lambda;
ill = 1; /* ill-conditioned? */
while (ill && (it<nit)) {
for (i=0; i<npar; i++)
h[i][i] = h[i][i]*mult;
lmstat->final_it = it;
lmstat->final_err = err;
lmstat->final_derr = derr;
ill = cholesky_decomp(npar, ch, h);
if (!ill) {
solve_axb_cholesky(npar, ch, delta, d);
for (i=0; i<npar; i++)
newpar[i] = par[i] + delta[i];
newerr = error_func(newpar, ny, y, dysq, func, fdata);
derr = newerr - err;
ill = (derr > 0);
}
if (ill) {
mult = (1 + lambda*up)/(1 + lambda);
lambda *= up;
it++;
}
if (it == nit) {
DEBUG_PRINT("did not converge");
return -1;
}
for (i=0; i<npar; i++)
par[i] = newpar[i];
err = newerr;
lambda *= down;
if ((!ill)&&(-derr<target_derr)) break;
}
lmstat->final_it = it;
lmstat->final_err = err;
lmstat->final_derr = derr;
return (it==nit);
return it;
}
/* calculate the error function (chi-squared) */
float error_func(float *par, int ny, float *y, float *dysq,
float (*func)(float *, int, void *), void *fdata)
float (*func)(float *, int, void *), void *fdata)
{
int x;
float res,e=0;
int x;
float res,e=0;
for (x=0; x<ny; x++) {
res = func(par, x, fdata) - y[x];
if (dysq) /* weighted least-squares */
e += res*res/dysq[x];
else
e += res*res;
}
return e;
for (x=0; x<ny; x++) {
res = func(par, x, fdata) - y[x];
if (dysq) /* weighted least-squares */
e += res*res/dysq[x];
else
e += res*res;
}
return e;
}
/* solve the equation Ax=b for a symmetric positive-definite matrix A,
using the Cholesky decomposition A=LL^T. The matrix L is passed in "l".
Elements above the diagonal are ignored.
*/
*/
void solve_axb_cholesky(int n, float l[n][n], float x[n], float b[n])
{
int i,j;
float sum;
int i,j;
float sum;
/* solve L*y = b for y (where x[] is used to store y) */
/* solve L*y = b for y (where x[] is used to store y) */
for (i=0; i<n; i++) {
sum = 0;
for (j=0; j<i; j++)
sum += l[i][j] * x[j];
x[i] = (b[i] - sum)/l[i][i];
}
for (i=0; i<n; i++) {
sum = 0;
for (j=0; j<i; j++)
sum += l[i][j] * x[j];
x[i] = (b[i] - sum)/l[i][i];
}
/* solve L^T*x = y for x (where x[] is used to store both y and x) */
/* solve L^T*x = y for x (where x[] is used to store both y and x) */
for (i=n-1; i>=0; i--) {
sum = 0;
for (j=i+1; j<n; j++)
sum += l[j][i] * x[j];
x[i] = (x[i] - sum)/l[i][i];
}
for (i=n-1; i>=0; i--) {
sum = 0;
for (j=i+1; j<n; j++)
sum += l[j][i] * x[j];
x[i] = (x[i] - sum)/l[i][i];
}
}
@ -190,26 +200,26 @@ void solve_axb_cholesky(int n, float l[n][n], float x[n], float b[n])
its (lower-triangular) Cholesky factor in "l". Elements above the
diagonal are neither used nor modified. The same array may be passed
as both l and a, in which case the decomposition is performed in place.
*/
*/
int cholesky_decomp(int n, float l[n][n], float a[n][n])
{
int i,j,k;
float sum;
int i,j,k;
float sum;
for (i=0; i<n; i++) {
for (j=0; j<i; j++) {
sum = 0;
for (k=0; k<j; k++)
sum += l[i][k] * l[j][k];
l[i][j] = (a[i][j] - sum)/l[j][j];
for (i=0; i<n; i++) {
for (j=0; j<i; j++) {
sum = 0;
for (k=0; k<j; k++)
sum += l[i][k] * l[j][k];
l[i][j] = (a[i][j] - sum)/l[j][j];
}
sum = 0;
for (k=0; k<i; k++)
sum += l[i][k] * l[i][k];
sum = a[i][i] - sum;
if (sum<TOL) return 1; /* not positive-definite */
l[i][i] = sqrtf(sum);
}
sum = 0;
for (k=0; k<i; k++)
sum += l[i][k] * l[i][k];
sum = a[i][i] - sum;
if (sum<TOL) return 1; /* not positive-definite */
l[i][i] = sqrtf(sum);
}
return 0;
return 0;
}

71
controller/fw/src/freq_meas.c
View file

@ -7,6 +7,7 @@
#include "freq_meas.h"
#include "sr_global.h"
#include "simulation.h"
/* FTT window lookup table defined in generated/fmeas_fft_window.c */
@ -29,17 +30,33 @@ extern arm_status arm_rfft_4096_fast_init_f32(arm_rfft_fast_instance_f32 * S);
void func_gauss_grad(float *out, float *params, int x, void *userdata);
float func_gauss(float *params, int x, void *userdata);
int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out) {
int adc_buf_measure_freq(int16_t adc_buf[FMEAS_FFT_LEN], float *out) {
int rc;
float in_buf[FMEAS_FFT_LEN];
float out_buf[FMEAS_FFT_LEN];
/*
DEBUG_PRINTN(" [emulated adc buf] ");
for (size_t i=0; i<FMEAS_FFT_LEN; i++)
DEBUG_PRINTN("%5d, ", adc_buf[i]);
DEBUG_PRINTN("\n");
*/
//DEBUG_PRINT("Applying window function");
for (size_t i=0; i<FMEAS_FFT_LEN; i++)
in_buf[i] = (float)adc_buf[i] / (float)FMEAS_ADC_MAX * fmeas_fft_window_table[i];
//DEBUG_PRINT("Running FFT");
arm_rfft_fast_instance_f32 fft_inst;
if ((rc = arm_rfft_init_name(FMEAS_FFT_LEN)(&fft_inst)) != ARM_MATH_SUCCESS)
if ((rc = arm_rfft_init_name(FMEAS_FFT_LEN)(&fft_inst)) != ARM_MATH_SUCCESS) {
*out = NAN;
return rc;
}
/*
DEBUG_PRINTN(" [input] ");
for (size_t i=0; i<FMEAS_FFT_LEN; i++)
DEBUG_PRINTN("%010f, ", in_buf[i]);
DEBUG_PRINTN("\n");
*/
arm_rfft_fast_f32(&fft_inst, in_buf, out_buf, 0);
#define FMEAS_FFT_WINDOW_MIN_F_HZ 30.0f
@ -49,10 +66,24 @@ int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out) {
const size_t last_bin = (int)(FMEAS_FFT_WINDOW_MAX_F_HZ / binsize_hz + 0.5f);
const size_t nbins = last_bin - first_bin + 1;
/* Copy real values of target data to front of output buffer */
for (size_t i=0; i<nbins; i++)
out_buf[i] = out_buf[2 * (first_bin + i)];
DEBUG_PRINT("binsize_hz=%f first_bin=%zd last_bin=%zd nbins=%zd", binsize_hz, first_bin, last_bin, nbins);
DEBUG_PRINTN(" [bins real] ");
for (size_t i=0; i<FMEAS_FFT_LEN/2; i+=2)
DEBUG_PRINTN("%010f, ", out_buf[i]);
DEBUG_PRINTN("\n [bins imag] ");
for (size_t i=1; i<FMEAS_FFT_LEN/2; i+=2)
DEBUG_PRINTN("%010f, ", out_buf[i]);
DEBUG_PRINT("\n");
DEBUG_PRINT("Repacking FFT results");
/* Copy real values of target data to front of output buffer */
for (size_t i=0; i<nbins; i++) {
float real = out_buf[2 * (first_bin + i)];
float imag = out_buf[2 * (first_bin + i) + 1];
out_buf[i] = sqrtf(real*real + imag*imag);
}
DEBUG_PRINT("Running Levenberg-Marquardt");
LMstat lmstat;
levmarq_init(&lmstat);
@ -68,27 +99,37 @@ int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out) {
float par[3] = {
a_max, i_max, 1.0f
};
DEBUG_PRINT(" par_pre={%010f, %010f, %010f}", par[0], par[1], par[2]);
if (levmarq(3, par, nbins, out_buf, NULL, func_gauss, func_gauss_grad, NULL, &lmstat))
if (levmarq(3, par, nbins, out_buf, NULL, func_gauss, func_gauss_grad, NULL, &lmstat) < 0) {
*out = NAN;
return -1;
}
DEBUG_PRINT(" par_post={%010f, %010f, %010f}", par[0], par[1], par[2]);
DEBUG_PRINT("done.");
*out = (par[1] + first_bin) * binsize_hz;
return 0;
}
float func_gauss(float *params, int x, void *userdata) {
UNUSED(userdata);
float a = params[0];
float mu = params[1];
float sigma = params[2];
return a*expf(-powf((x-mu), 2.0f/(2.0f*(sigma*sigma))));
float a = params[0], b = params[1], c = params[2];
float n = x-b;
return a*expf(-n*n / (2.0f* c*c));
}
void func_gauss_grad(float *out, float *params, int x, void *userdata) {
UNUSED(userdata);
float a = params[0];
float mu = params[1];
float sigma = params[2];
*out = -(x-mu) / ( sigma*sigma*sigma * 2.5066282746310002f) * a*expf(-powf((x-mu), 2.0f/(2.0f*(sigma*sigma))));
float a = params[0], b = params[1], c = params[2];
float n = x-b;
float e = expf(-n*n / (2.0f * c*c));
/* d/da */
out[0] = e;
/* d/db */
out[1] = a*n/(c*c) * e;
/* d/dc */
out[2] = a*n*n/(c*c*c) * e;
}

2
controller/fw/src/freq_meas.h
View file

@ -2,6 +2,6 @@
#ifndef __FREQ_MEAS_H__
#define __FREQ_MEAS_H__
int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out);
int adc_buf_measure_freq(int16_t adc_buf[FMEAS_FFT_LEN], float *out);
#endif /* __FREQ_MEAS_H__ */

14
controller/fw/src/simulation.h Normal file
View file

@ -0,0 +1,14 @@
#ifndef __SIMULATION_H__
#define __SIMULATION_H__
#ifdef SIMULATION
#include <stdio.h>
#define DEBUG_PRINTN(...) fprintf(stderr, __VA_ARGS__)
#define DEBUG_PRINTNF(fmt, ...) DEBUG_PRINTN("%s:%d: " fmt, __FILE__, __LINE__, ##__VA_ARGS__)
#define DEBUG_PRINT(fmt, ...) DEBUG_PRINTNF(fmt "\n", ##__VA_ARGS__)
#else
#define DEBUG_PRINT(...) ((void)0)
#define DEBUG_PRINTN(...) ((void)0)
#endif
#endif /* __SIMULATION_H__ */

38
controller/fw/tools/freq_meas_test.c
View file

@ -1,4 +1,6 @@
#include <stdint.h>
#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
@ -13,7 +15,7 @@
void print_usage(void);
void print_usage() {
fprintf(stderr, "Usage: freq_meas_test [test_data.bin]");
fprintf(stderr, "Usage: freq_meas_test [test_data.bin]\n");
}
int main(int argc, char **argv) {
@ -46,6 +48,7 @@ int main(int argc, char **argv) {
return 2;
}
fprintf(stderr, "Reading %zd samples test data...", st.st_size/sizeof(float));
size_t nread = 0;
while (nread < st.st_size) {
ssize_t rc = read(fd, buf, st.st_size - nread);
@ -54,41 +57,48 @@ int main(int argc, char **argv) {
continue;
if (rc < 0) {
fprintf(stderr, "Error reading test data: %s\n", strerror(errno));
fprintf(stderr, "\nError reading test data: %s\n", strerror(errno));
return 2;
}
if (rc == 0) {
fprintf(stderr, "Error reading test data: Unexpected end of file\n");
fprintf(stderr, "\nError reading test data: Unexpected end of file\n");
return 2;
}
nread += rc;
}
fprintf(stderr, " done.\n");
size_t n_samples = st.st_size / sizeof(float);
float *buf_f = (float *)buf;
uint16_t *sim_adc_buf = calloc(sizeof(uint16_t), n_samples);
int16_t *sim_adc_buf = calloc(sizeof(int16_t), n_samples);
if (!sim_adc_buf) {
fprintf(stderr, "Error allocating memory\n");
return 2;
}
fprintf(stderr, "Converting and truncating test data...");
for (size_t i=0; i<n_samples; i++)
sim_adc_buf[i] = 2048 + buf_f[i] * 2047;
/* Note on scaling: We can't simply scale by 0x8000 (1/2 full range) here. Our test data is nominally 1Vp-p but
* certain tests such as the interharmonics one can have some samples exceeding that range. */
sim_adc_buf[i] = buf_f[i] * (0x4000-1);
fprintf(stderr, " done.\n");
for (size_t i=0; i<n_samples; i+=FMEAS_FFT_LEN) {
fprintf(stderr, "Starting simulation.\n");
float out;
int rc = adc_buf_measure_freq(sim_adc_buf + i, &out);
if (rc) {
fprintf(stderr, "Simulation error in iteration %zd at position %zd: %d\n", i/FMEAS_FFT_LEN, i, rc);
return 3;
}
size_t iterations = (n_samples-FMEAS_FFT_LEN)/(FMEAS_FFT_LEN/2);
for (size_t i=0; i<iterations; i++) {
printf("%09zd %015f\n", i, out);
fprintf(stderr, "Iteration %zd/%zd\n", i, iterations);
float res = NAN;
int rc = adc_buf_measure_freq(sim_adc_buf + i*(FMEAS_FFT_LEN/2), &res);
if (rc)
printf("ERROR: Simulation error in iteration %zd at position %zd: %d\n", i, i*(FMEAS_FFT_LEN/2), rc);
printf("%09zd %12f\n", i, res);
}
return 0;
}

39
controller/fw/tools/freq_meas_test_runner.py Normal file
View file

@ -0,0 +1,39 @@
#!/usr/bin/env python3
import os
from os import path
import subprocess
import json
import numpy as np
np.set_printoptions(linewidth=240)
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument(metavar='test_data_directory', dest='dir', help='Directory with test data .bin files')
default_binary = path.abspath(path.join(path.dirname(__file__), '../build/tools/freq_meas_test'))
parser.add_argument(metavar='test_binary', dest='binary', nargs='?', default=default_binary)
parser.add_argument('-d', '--dump', help='Write raw measurements to JSON file')
args = parser.parse_args()
bin_files = [ path.join(args.dir, d) for d in os.listdir(args.dir) if d.lower().endswith('.bin') ]
savedata = {}
for p in bin_files:
output = subprocess.check_output([args.binary, p], stderr=subprocess.DEVNULL)
measurements = np.array([ float(value) for _offset, value in [ line.split() for line in output.splitlines() ] ])
savedata[p] = list(measurements)
# Cut off first and last sample for mean and RMS calculations as these show boundary effects.
measurements = measurements[1:-1]
mean = np.mean(measurements)
rms = np.sqrt(np.mean(np.square(measurements - mean)))
print(f'{path.basename(p):<60}: mean={mean:<8.4f}Hz rms={rms*1000:.3f}mHz')
if args.dump:
with open(args.dump, 'w') as f:
json.dump(savedata, f)

396
lab-windows/Phase Measurement Prototype.ipynb Normal file
View file

@ -0,0 +1,396 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import struct\n",
"\n",
"import numpy as np\n",
"from scipy import signal, optimize\n",
"from matplotlib import pyplot as plt\n",
"\n",
"import rocof_test_data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"fs = 1000 # Hz\n",
"ff = 50 # Hz\n",
"duration = 60 # seconds\n",
"# test_data = rocof_test_data.sample_waveform(rocof_test_data.test_close_interharmonics_and_flicker(),\n",
"# duration=20,\n",
"# sampling_rate=fs,\n",
"# frequency=ff)[0]\n",
"# test_data = rocof_test_data.sample_waveform(rocof_test_data.gen_noise(fmin=10, amplitude=1),\n",
"# duration=20,\n",
"# sampling_rate=fs,\n",
"# frequency=ff)[0]\n",
"\n",
"\n",
"#gen = rocof_test_data.gen_noise(fmin=10, amplitude=1)\n",
"# gen = rocof_test_data.gen_noise(fmin=60, amplitude=0.2)\n",
"# gen = rocof_test_data.test_harmonics()\n",
"# gen = rocof_test_data.gen_interharmonic(*rocof_test_data.test_interharmonics)\n",
"# gen = rocof_test_data.test_amplitude_steps()\n",
"# gen = rocof_test_data.test_amplitude_and_phase_steps()\n",
"test_data = []\n",
"test_labels = [ fun.__name__.replace('test_', '') for fun in rocof_test_data.all_tests ]\n",
"for gen in rocof_test_data.all_tests:\n",
" test_data.append(rocof_test_data.sample_waveform(gen(),\n",
" duration=duration,\n",
" sampling_rate=fs,\n",
" frequency=ff)[0])\n",
"# d = 10 # seconds\n",
"# test_data = np.sin(2*np.pi * ff * np.linspace(0, d, int(d*fs)))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"spr_fmt = f'{fs}Hz' if fs<1000 else f'{fs/1e3:f}'.rstrip('.0') + 'kHz'\n",
"for label, data in zip(test_labels, test_data):\n",
" with open(f'rocof_test_data/rocof_test_{label}_{spr_fmt}.bin', 'wb') as f:\n",
" for sample in data:\n",
" f.write(struct.pack('<f', sample))"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [],
"source": [
"analysis_periods = 10\n",
"window_len = 256 # fs * analysis_periods/ff\n",
"nfft_factor = 1\n",
"sigma = window_len/8 # samples\n",
"quantization_bits = 14\n",
"\n",
"ffts = []\n",
"for item in test_data:\n",
" f, t, Zxx = signal.stft((item * (2**(quantization_bits-1) - 1)).round().astype(np.int16).astype(float),\n",
" fs = fs,\n",
" window=('gaussian', sigma),\n",
" nperseg = window_len,\n",
" nfft = window_len * nfft_factor)\n",
" #boundary = 'zeros')\n",
" ffts.append((f, t, Zxx))"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(129, 470)"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Zxx.shape"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.90625"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1000/256"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-61-530955947ba4>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig, ax = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "30b7e529e3f94c1b8241c6b6760b7e69",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n",
"fig.tight_layout(pad=2, h_pad=0.1)\n",
"\n",
"for fft, ax, label in zip(test_data, ax.flatten(), test_labels):\n",
" ax.plot((item * (2**(quantization_bits-1) - 1)).round())\n",
" \n",
" ax.set_title(label, pad=-20, color='white', bbox=dict(boxstyle=\"square\", ec=(0,0,0,0), fc=(0,0,0,0.8)))\n",
" ax.grid()\n",
" ax.set_ylabel('f [Hz]')\n",
"ax.set_xlabel('simulation time t [s]')\n",
"ax.set_xlim([5000, 5200])\n",
"None"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-76-31c82486a777>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig, ax = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "9e0976ade81c4990a10fd1182f0f20d5",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n",
"fig.tight_layout(pad=2, h_pad=0.1)\n",
"\n",
"for fft, ax, label in zip(ffts, ax.flatten(), test_labels):\n",
" f, t, Zxx = fft\n",
" ax.pcolormesh(t[1:], f[:250], np.abs(Zxx[:250,1:]))\n",
" ax.set_title(label, pad=-20, color='white')\n",
" ax.grid()\n",
" ax.set_ylabel('f [Hz]')\n",
" ax.set_ylim([30, 75]) # Hz\n",
"ax.set_xlabel('simulation time t [s]')\n",
"None"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0. , 3.90625, 7.8125 , 11.71875, 15.625 , 19.53125,\n",
" 23.4375 , 27.34375, 31.25 , 35.15625, 39.0625 , 42.96875,\n",
" 46.875 , 50.78125, 54.6875 , 58.59375, 62.5 , 66.40625,\n",
" 70.3125 , 74.21875, 78.125 , 82.03125, 85.9375 , 89.84375,\n",
" 93.75 , 97.65625, 101.5625 , 105.46875, 109.375 , 113.28125,\n",
" 117.1875 , 121.09375, 125. , 128.90625, 132.8125 , 136.71875,\n",
" 140.625 , 144.53125, 148.4375 , 152.34375, 156.25 , 160.15625,\n",
" 164.0625 , 167.96875, 171.875 , 175.78125, 179.6875 , 183.59375,\n",
" 187.5 , 191.40625, 195.3125 , 199.21875, 203.125 , 207.03125,\n",
" 210.9375 , 214.84375, 218.75 , 222.65625, 226.5625 , 230.46875,\n",
" 234.375 , 238.28125, 242.1875 , 246.09375, 250. , 253.90625,\n",
" 257.8125 , 261.71875, 265.625 , 269.53125, 273.4375 , 277.34375,\n",
" 281.25 , 285.15625, 289.0625 , 292.96875, 296.875 , 300.78125,\n",
" 304.6875 , 308.59375, 312.5 , 316.40625, 320.3125 , 324.21875,\n",
" 328.125 , 332.03125, 335.9375 , 339.84375, 343.75 , 347.65625,\n",
" 351.5625 , 355.46875, 359.375 , 363.28125, 367.1875 , 371.09375,\n",
" 375. , 378.90625, 382.8125 , 386.71875, 390.625 , 394.53125,\n",
" 398.4375 , 402.34375, 406.25 , 410.15625, 414.0625 , 417.96875,\n",
" 421.875 , 425.78125, 429.6875 , 433.59375, 437.5 , 441.40625,\n",
" 445.3125 , 449.21875, 453.125 , 457.03125, 460.9375 , 464.84375,\n",
" 468.75 , 472.65625, 476.5625 , 480.46875, 484.375 , 488.28125,\n",
" 492.1875 , 496.09375, 500. ])"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-63-888d30b0f7d6>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
" fig, axs = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "cef58f5753084c54a648418e134775d8",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axs = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n",
"fig.tight_layout(pad=2.2, h_pad=0, w_pad=1)\n",
"\n",
"for fft, ax, label in zip(ffts, axs.flatten(), test_labels):\n",
" f, f_t, Zxx = fft\n",
" \n",
" n_f, n_t = Zxx.shape\n",
" # start, stop = 180, 220\n",
" # start, stop = 90, 110\n",
" # start, stop = 15, 35\n",
" # bounds_f = slice(start // 4 * nfft_factor, stop // 4 * nfft_factor)\n",
" f_min, f_max = 30, 70 # Hz\n",
" bounds_f = slice(np.argmax(f > f_min), np.argmin(f < f_max))\n",
" \n",
"\n",
" f_mean = np.zeros(Zxx.shape[1])\n",
" for t in range(1, Zxx.shape[1] - 1):\n",
" frame_f = f[bounds_f]\n",
" frame_step = frame_f[1] - frame_f[0]\n",
" time_step = f_t[1] - f_t[0]\n",
" #if t == 10:\n",
" # axs[-1].plot(frame_f, frame_Z)\n",
" frame_Z = np.abs(Zxx[bounds_f, t])\n",
" # frame_f = f[180:220]\n",
" # frame_Z = np.abs(Zxx[180:220, 40])\n",
" # frame_f = f[15:35]\n",
" # frame_Z = np.abs(Zxx[15:35, 40])\n",
" # plt.plot(frame_f, frame_Z)\n",
"\n",
" # peak_f = frame_f[np.argmax(frame)]\n",
" # plt.axvline(peak_f, color='red')\n",
"\n",
"# def gauss(x, *p):\n",
"# A, mu, sigma, o = p\n",
"# return A*np.exp(-(x-mu)**2/(2.*sigma**2)) + o\n",
" \n",
" def gauss(x, *p):\n",
" A, mu, sigma = p\n",
" return A*np.exp(-(x-mu)**2/(2.*sigma**2))\n",
" \n",
" f_start = frame_f[np.argmax(frame_Z)]\n",
" A_start = np.max(frame_Z)\n",
" p0 = [A_start, f_start, 1.]\n",
" try:\n",
" coeff, var = optimize.curve_fit(gauss, frame_f, frame_Z, p0=p0)\n",
" # plt.plot(frame_f, gauss(frame_f, *coeff))\n",
" #print(coeff)\n",
" A, mu, sigma, *_ = coeff\n",
" f_mean[t] = mu\n",
" except RuntimeError:\n",
" f_mean[t] = np.nan\n",
" ax.plot(f_t[1:-1], f_mean[1:-1])\n",
" \n",
"# b, a = signal.butter(3,\n",
"# 1/5, # Hz\n",
"# btype='lowpass',\n",
"# fs=1/time_step)\n",
"# filtered = signal.lfilter(b, a, f_mean[1:-1], axis=0)\n",
"# ax.plot(f_t[1:-1], filtered)\n",
" \n",
" ax.set_title(label, pad=-20)\n",
" ax.set_ylabel('f [Hz]')\n",
" ax.grid()\n",
" if not label in ['off_frequency', 'sweep_phase_steps']:\n",
" ax.set_ylim([49.90, 50.10])\n",
" var = np.var(f_mean[1:-1])\n",
" ax.text(0.5, 0.1, f'σ²={var * 1e3:.3g} mHz²', transform=ax.transAxes, ha='center')\n",
" ax.text(0.5, 0.25, f'σ={np.sqrt(var) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n",
"# ax.text(0.5, 0.2, f'filt. σ²={np.var(filtered) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n",
" else:\n",
" f_min, f_max = min(f_mean[1:-1]), max(f_mean[1:-1])\n",
" delta = f_max - f_min\n",
" ax.set_ylim(f_min - delta * 0.1, f_max + delta * 0.3)\n",
" \n",
"ax.set_xlabel('simulation time t [s]')\n",
"None"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "labenv",
"language": "python",
"name": "labenv"
},
"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.8.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}

184
lab-windows/gauss_derivations.ipynb Normal file
View file

@ -0,0 +1,184 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{2}$}f\\left(a , \\mu , \\sigma , x\\right):=a\\,\\exp \\left(\\frac{-\\left(x-\\mu\\right)^2}{2\\,\\sigma^2}\\right)\\]"
],
"text/plain": [
" 2\n",
" - (x - mu)\n",
"(%o2) f(a, mu, sigma, x) := a exp(-----------)\n",
" 2\n",
" 2 sigma"
],
"text/x-maxima": [
"f(a,mu,sigma,x):=a*exp((-(x-mu)^2)/(2*sigma^2))"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(a, mu, sigma, x) := a*exp(-(x-mu)^2/(2*sigma^2));"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{11}$}-\\frac{a\\,\\left(x-\\mu\\right)\\,e^ {- \\frac{\\left(x-\\mu\\right)^2}{2\\,\\sigma^2} }}{\\sigma^2}\\]"
],
"text/plain": [
" 2\n",
" (x - mu)\n",
" - ---------\n",
" 2\n",
" 2 sigma\n",
" a (x - mu) %e\n",
"(%o11) - ------------------------\n",
" 2\n",
" sigma"
],
"text/x-maxima": [
"-(a*(x-mu)*%e^-((x-mu)^2/(2*sigma^2)))/sigma^2"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(f(a, mu, sigma, x), x);"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{12}$}e^ {- \\frac{\\left(x-\\mu\\right)^2}{2\\,\\sigma^2} }\\]"
],
"text/plain": [
" 2\n",
" (x - mu)\n",
" - ---------\n",
" 2\n",
" 2 sigma\n",
"(%o12) %e"
],
"text/x-maxima": [
"%e^-((x-mu)^2/(2*sigma^2))"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(f(a, mu, sigma, x), a);"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{13}$}\\frac{a\\,\\left(x-\\mu\\right)\\,e^ {- \\frac{\\left(x-\\mu\\right)^2}{2\\,\\sigma^2} }}{\\sigma^2}\\]"
],
"text/plain": [
" 2\n",
" (x - mu)\n",
" - ---------\n",
" 2\n",
" 2 sigma\n",
" a (x - mu) %e\n",
"(%o13) ------------------------\n",
" 2\n",
" sigma"
],
"text/x-maxima": [
"(a*(x-mu)*%e^-((x-mu)^2/(2*sigma^2)))/sigma^2"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(f(a, mu, sigma, x), mu);"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{14}$}\\frac{a\\,\\left(x-\\mu\\right)^2\\,e^ {- \\frac{\\left(x-\\mu\\right)^2}{2\\,\\sigma^2} }}{\\sigma^3}\\]"
],
"text/plain": [
" 2\n",
" (x - mu)\n",
" - ---------\n",
" 2\n",
" 2 2 sigma\n",
" a (x - mu) %e\n",
"(%o14) -------------------------\n",
" 3\n",
" sigma"
],
"text/x-maxima": [
"(a*(x-mu)^2*%e^-((x-mu)^2/(2*sigma^2)))/sigma^3"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(f(a, mu, sigma, x), sigma);"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Maxima",
"language": "maxima",
"name": "maxima"
},
"language_info": {
"codemirror_mode": "maxima",
"file_extension": ".mac",
"mimetype": "text/x-maxima",
"name": "maxima",
"pygments_lexer": "maxima",
"version": "5.43.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}

1
lab-windows/impl_test_out.json Normal file
View file

File diff suppressed because one or more lines are too long

174
lab-windows/rocof_test_data.py Normal file
View file

@ -0,0 +1,174 @@
#!/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[ ]:
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]

BIN
lab-windows/rocof_test_data/rocof_test_amplitude_and_phase_steps_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_amplitude_steps_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_close_interharmonics_and_flicker_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_harmonics_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_interharmonics_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_noise_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_noise_loud_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_off_frequency_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_phase_steps_1kHz.bin Normal file
View file

Binary file not shown.

BIN
lab-windows/rocof_test_data/rocof_test_sweep_phase_steps_1kHz.bin Normal file
View file

Binary file not shown.

167
lab-windows/scratch.ipynb Normal file
View file

@ -0,0 +1,167 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import json\n",
"\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"import matplotlib"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "712481c28d9d4e1d874a66d31c3e8bff",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"(0, 64)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"a = np.array([-00.000732, -00.000352, -00.000666, -00.000202, -00.000706, -00.000006, -00.000597, -00.002039, 000.050663, -00.644566, 004.456614, -16.817095, 034.654587, -39.021217, 024.007816, -08.070650, 001.478795, -00.150260, 000.006110, -00.002328, -00.002322, -00.002426, -00.002177, -00.002452, -00.002333, -00.002438, -00.002342, -00.002396, -00.001979, -00.003049, -00.001720, -00.002686, -00.002168, -00.002507, -00.001868, -00.002899, -00.002017, -00.001952, -00.003255, -00.001080, -00.003335, -00.001575, -00.002704, -00.001872, -00.002735, -00.001983, -00.002191, -00.002478, -00.002155, -00.002203, -00.002328, -00.002206, -00.002443, -00.001770, -00.002718, -00.002004, -00.002378, -00.002112, -00.002122, -00.002691, -00.001679, -00.002690, -00.001946, -00.002232])\n",
"b = np.array([-00.002734, -00.001325, -00.002220, -00.003693, -00.004907, -00.006454, -00.007737, 000.004823, -00.363143, 004.688968, -33.795303, 130.992630, -274.092651, 309.377991, -188.427826, 061.912941, -10.974002, 001.053608, -00.048927, 000.007710, 000.007010, 000.006493, 000.007234, 000.006725, 000.006938, 000.006694, 000.006356, 000.006173, 000.006333, 000.005684, 000.005697, 000.005575, 000.005101, 000.005693, 000.004319, 000.005344, 000.004673, 000.003566, 000.006213, 000.002719, 000.004850, 000.003755, 000.004243, 000.003419, 000.003960, 000.003498, 000.003297, 000.003877, 000.002836, 000.003487, 000.003144, 000.002824, 000.003355, 000.002528, 000.002975, 000.003012, 000.002137, 000.003112, 000.002416, 000.002512, 000.002084, 000.003008, 000.001837, 000.002351])\n",
"ax.plot([3.906250*i for i in range(len(a))], np.sqrt(a**2 + b**2))\n",
"a2 = ax.twiny()\n",
"a2.set_xlim([0, len(a)])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "d4024377df494eac935fd487026edc8b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f2eb410f280>]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"d = [50.000839,50.000839,50.000832,50.000824,50.000839,50.000832,50.000839,50.000824,50.000847,50.000824,50.000824,50.000839,50.000832,50.000839,50.000824,50.000824,50.000839,50.000824,50.000824,50.000835,50.000816,50.000832,50.000847,50.000832,50.000835,50.000824,50.000824,50.000832,50.000832,50.000843,50.000824,50.000832,50.000832,50.000832,50.000828,50.000832,50.000832,50.000824,50.000816,50.000835,50.000843,50.000824,50.000824,50.000832,50.000832,50.000847,50.000824,50.000824,50.000824,50.000835,50.000835,50.000851,50.000824,50.000824,50.000832,50.000828,50.000828,50.000824,50.000832,50.000835,50.000835,50.000832,50.000847,50.000824,50.000832,50.000839,50.000839,50.000824,50.000832,50.000832,50.000832,50.000835,50.000816,50.000820,50.000824,50.000832,50.000824,50.000832,50.000835,50.000832,50.000816,50.000820,50.000839,50.000839,50.000824,50.000839,50.000820,50.000820,50.000839,50.000832,50.000835,50.000828,50.000824,50.000839,50.000839,50.000839,50.000816,50.000832,50.000824,50.000832,50.000832,50.000839,50.000824,50.000832,50.000828,50.000832,50.000828,50.000835,50.000832,50.000843,50.000839,50.000820,50.000832,50.000835,50.000824,50.000824,50.000828,50.000820,50.000820,50.000828,50.000832,50.000832,50.000828,50.000835,50.000839,50.000820,50.000832,50.000832,50.000824,50.000832,50.000832,50.000839,50.000839,50.000816,50.000828,50.000832,50.000839,50.000824,50.000824,50.000824,50.000835,50.000824,50.000832,50.000839,50.000835,50.000832,50.000828,50.000835,50.000828,50.000828,50.000824,50.000824,50.000839,50.000832,50.000824,50.000832,50.000832,50.000820,50.000851,50.000824,50.000824,50.000839,50.000824,50.000839,50.000832,50.000835,50.000820,50.000832,50.000839,50.000832,50.000832,50.000824,50.000832,50.000824,50.000832,50.000839,50.000839,50.000832,50.000816,50.000835,50.000854,50.000824,50.000816,50.000832,50.000832,50.000835,50.000816,50.000832,50.000824,50.000832,50.000832,50.000832,50.000824,50.000832,50.000824,50.000835,50.000832,50.000835,50.000832,50.000832,50.000828,50.000839,50.000824,50.000839,50.000824,50.000824,50.000839,50.000816,50.000839,50.000816,50.000832,50.000839,50.000839,50.000832,50.000824,50.000832,50.000820,50.000824,50.000835,50.000824,50.000835,50.000832,50.000824,50.000824,50.000820,50.000839,50.000816,50.000832,50.000832,50.000832,50.000824,50.000847,50.000824,50.000839]\n",
"\n",
"ax.plot(d)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"with open('impl_test_out.json') as f:\n",
" impl_measurements = json.load(f)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "dd23cf23221e4e14aaafdd58bb9416d9",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axs = plt.subplots(len(impl_measurements), figsize=(8, 20), sharex=True)\n",
"fig.tight_layout()\n",
"axs = axs.flatten()\n",
"\n",
"for (label, data), ax in zip(impl_measurements.items(), axs):\n",
" ax.set_title(label)\n",
" ax.plot(data[1:-1])\n",
" mean = np.mean(data[1:-1])\n",
" rms = np.sqrt(np.mean(np.square(data[1:-1] - mean)))\n",
" ax.text(0.2, 0.2, f'mean={mean:.3}Hz, rms={rms*1e3:.3}mHz', ha='center', va='center', transform=ax.transAxes,\n",
" bbox=dict(boxstyle=\"square\", ec=(0,0,0,0), fc=(1,1,1,0.8)))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "labenv",
"language": "python",
"name": "labenv"
},
"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.8.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}