Day 8: (Sep 24) Signals

Notes for 2018-09-24. See also the Fall 2018 Calendar.

Notes from Day 7

Agenda

  1. Administrative
    • Clearances!
  2. Assignments
  3. In-class
    1. Demo 4 review.
      1. Each pair will demo and explain their puppet devices to another pair. We’ll take 5 minutes for this before starting the checkoff round.
      2. During the checkoff review, the authors will operate the device and the the opposite pair will ask questions and critique the work. Critique, not criticism.
        1. Herb Simon: “Learning results from what the student does and thinks and only from what the student does and thinks. The teacher can advance learning only by influencing what the student does to learn.”
        2. Responders, please start with simply describing what you see in the work.
        3. The authors will then pose an specific, unresolved question which would help further development.
        4. The responders will then ask neutral questions which may help reveal the work.
      3. The following prompt questions may help spark your inquiries:
        1. How would you characterize the personality of the device?
        2. What elements appear scripted? What behaviors are clearly reactive?
        3. Is the sensing logically integrated with the form?
        4. Does the experience feel responsive?
    2. Introduce demo 5 assignment and pairings.
    3. Discussion of signal processing.
      1. Sampling, resolution, and noise.
      2. Low-pass filtering.
      3. Velocity estimation.
      4. Feedback and PID control.

Lecture code samples

First-order smoothing

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
/****************************************************************/
void setup(void)
{
  Serial.begin(115200);
}

/****************************************************************/
void loop(void)
{
  static unsigned long last_update_clock = 0;

  unsigned long now = micros();
  unsigned long interval = now - last_update_clock;
  last_update_clock = now;

  const long sample_interval = 10000;   // 10 msec, 100 Hz sampling
  static long sample_timer = 0;
  
  sample_timer -= interval;
  if (sample_timer <= 0) {
    sample_timer += sample_interval;

    int raw_value = analogRead(0);                       // read the current input 
    float calibrated = ((float)raw_value) / 1024.0;      // scale to normalized units 

    static float smoothed_value = 0.0;                   // filtered value of the input
    float difference   =  calibrated - smoothed_value;   // compute the 'error'
    smoothed_value     += 0.1 * difference;              // apply a constant gain to move the smoothed value toward the reading

    Serial.println(smoothed_value);
  }
}

Ring filter

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
const unsigned LENGTH = 20;
float buffer[LENGTH];       // circular buffer of samples
unsigned position = 0;      // index of the oldest sample



// Put a new value into a circular sample buffer, overwriting the oldest sample.
void ringfilter_put(float value)
{
  if (position < LENGTH) buffer[position] = value;
  if (++position >= LENGTH) position = 0;
}


// Calculate the first derivative as the finite difference between the newest
// and oldest values.
float ringfilter_deriv(float delta_t)
{
  float oldest = buffer[position];
  float newest = (position < LENGTH-1) ? buffer[position+1] : buffer[0];
  return (newest - oldest) / ((LENGTH-1) * delta_t);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
/****************************************************************/
// The following was generated using ringfilter_gen.py

// Ring filter coefficients for quadratic.
static const float coeff_dt  = 0.010000; // sampling interval
static const float coeff_end = 0.000000; // sample end time
static const int coeff_order = 3;
static const int coeff_length = 20;
static float coeff[3][20] = {
  {   0.09935,    0.05519,    0.01753,   -0.01364,   -0.03831,   -0.05649,   -0.06818,   -0.07338, 
     -0.07208,   -0.06429,   -0.05000,   -0.02922,   -0.00195,    0.03182,    0.07208,    0.11883, 
      0.17208,    0.23182,    0.29805,    0.37078 },
  {   4.74026,    2.94258,    1.36136,   -0.00342,   -1.15174,   -2.08362,   -2.79904,   -3.29802, 
     -3.58054,   -3.64662,   -3.49624,   -3.12941,   -2.54614,   -1.74641,   -0.73023,    0.50239, 
      1.95147,    3.61700,    5.49897,    7.59740 },
  {  32.46753,   22.21463,   13.10093,    5.12645,   -1.70882,   -7.40488,  -11.96172,  -15.37936, 
    -17.65778,  -18.79699,  -18.79699,  -17.65778,  -15.37936,  -11.96172,   -7.40488,   -1.70882, 
      5.12645,   13.10093,   22.21463,   32.46753 }
};

/****************************************************************/
// Multiply the filter coefficient matrix by the circular sample buffer to
// produce a result vector.
void ringfilter_fit_quadratic(float result[3])
{
  unsigned index = position;

  // iterate over the coefficient rows
  for (int i = 0; i < coeff_order; i++) {    
    result[i] = 0.0; // clear accumulator

    // Iterate over the samples and the coefficient rows.  The index cycles
    // around the circular buffer once per row.
    for (int j = 0; j < coeff_length; j++) {
      result[i] += coeff[i][j] * buffer[index];
      if (++index >= LENGTH) index = 0;
    }
  }
}
/****************************************************************/

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
/****************************************************************/
void setup(void)
{
  Serial.begin(115200);
}

/****************************************************************/
void loop(void)
{
  static unsigned long last_update_clock = 0;

  unsigned long now = micros();
  unsigned long interval = now - last_update_clock;
  last_update_clock = now;

  const long sample_interval = 10000;   // 10 msec, 100 Hz sampling
  static long sample_timer = 0;
  
  sample_timer -= interval;
  if (sample_timer <= 0) {
    sample_timer += sample_interval;

    int raw_value = analogRead(0);
    float calibrated = ((float)raw_value) / 1024.0;

    ringfilter_put(calibrated);

    float velocity = ringfilter_deriv(0.01);
    
    float quadratic_params[3];
    ringfilter_fit_quadratic(quadratic_params);

    Serial.print("v: "); Serial.print(velocity);
    Serial.print(" q: "); Serial.print(quadratic_params[0]);
    Serial.print(" "); Serial.print(quadratic_params[1]);
    Serial.print(" "); Serial.println(quadratic_params[2]);
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#!/usr/bin/env python

# Generate a coefficient table for a ringfilter. This solves an
# overconstrained linear system using a pseudoinverse based on SVD.  This
# represents a least-squares parametric fit to uniformly sampled data.

# E.g. a quadratic:
#
#    y = a0 + a1*t + a2*t*t
#
# can be sampled in time to be represented as a matrix F multiplied by
# the quadratic coefficients:
#
#    y = F x
#
# where
#    y = []               a Nx1 matrix vector of measurements over time
#    F = [1 t t**2]       a Nx3 matrix of polynomial terms for a set of time values
#    x = [a0 a1 a2]       a 3x1 matrix of the unknown coefficients
#
# The solution for x uses the pseudoinverse F+:
#
#   x = F+ y, where F+ is the peudoinverse

import sys, argparse
import numpy as np
import scipy.linalg

################################################################
# Generate a fitting matrix which can compute the weight terms (a0, a1, a2) for
# a quadratic polynomial by left-multiplying with a vector of uniformly sampled
# y values.  In the ring buffer, the values are ordered from oldest to newest,
# with the newest by default considered to be at t=0.0.  The polynomial basis
# terms are evaluated over the negative time range corresponding to the ring
# buffer data, so the default is a causal filter estimating the parameters for
# the quadratic with t=0 at the present moment.  The optional 'endt' parameter
# can change the assumed endpoint to t=endt, e.g. to make an acausal filter with
# the quadratic t=0 centered in the buffer data.

# The quadratic function has the following linear form:
#  y(t) = [1 t t**2] * [a0 a1 a2]'

def quadratic_fit_matrix(dt = 0.005, N = 50, endt = 0.0):
    t0 = np.ones(N)
    t1 = np.linspace(endt - (N-1)*dt, endt, N, endpoint=True)
    t2 = t1*t1
    F = np.vstack((t0,t1,t2)).transpose()
    Finv = scipy.linalg.pinv2(F)
    return Finv

################################################################
# Write out a C declaration for a fitting matrix in a
# form convenient for use with the ringfilter

def write_fitting_matrix_in_C(stream, F, name="coeff_e", precision='float'):
    order = F.shape[0]
    N     = F.shape[1]
    stream.write("static const int %s_order = %d;\nstatic const int %s_length = %d;\n" % (name, order, name, N))
    stream.write("static %s %s[%d][%d] = {\n" % (precision, name, order, N))    
    for i in range(order):
        stream.write("  {")
        for j,v in enumerate(F[i][:]):
            stream.write("%10.5f" % v)
            if j < N-1:
                stream.write(", ")
                if j%8 == 7: stream.write("\n   ");
        if i < order-1: stream.write(" },\n")
        else: stream.write(" }\n")
    stream.write("};\n")
        
################################################################
if __name__ == "__main__":
    parser = argparse.ArgumentParser(description = """Generate fitting matrices for use with ringfilter.""")
    parser.add_argument('--dt', default=0.01, type=float, help = 'Sampling interval in seconds.')
    parser.add_argument('-N', default=20, type=int, help = 'Number of samples.')
    parser.add_argument('--end', default=0.0, type=float, help = 'Time reference for end sample.')
    parser.add_argument('--name', default='coeff', type=str, help = 'C name of fitting matrix.')
    parser.add_argument( '--out', type=argparse.FileType('w'), default=None, help='Path of C output file.')
    args = parser.parse_args()

    F = quadratic_fit_matrix(dt = args.dt, N = args.N, endt = args.end)
    stream = sys.stdout if args.out is None else args.out
    stream.write("// Ring filter coefficients for quadratic.\n")
    stream.write("static const float %s_dt  = %f; // sampling interval\n" % (args.name, args.dt))
    stream.write("static const float %s_end = %f; // sample end time\n" % (args.name, args.end))
    write_fitting_matrix_in_C(stream, F, args.name)