Saturday, October 3, 2015

Signals, Symbolic Assignments and the Plotter


Symbolic assignments are mostly used to adjust readings coming from connected devices or to define so called virtual devices. In this post, Rafael Cobo will get you familiar with the extended symbolic assignments he just implemented for the plotter in Artisan v0.9.8.

The plotter is a mechanism that allows to show additional curves defined by symbolic formulas based on data taken from the current profile. To get started we will first dive into the theory of profile representation, before we introduce the new plotter features and motivate them through examples.

Signals


In Artisan, profiles are mathematical representation of signals. Profiles are finite data structures arranged in one dimensional arrays, whose elements are indexed by a subscript. Suppose that the temperatures acquired in a profile are listed as follows:

Y = 300, 304, 307, 308, 308, … 

One can denote all the values in the array using only one symbol, for example Y, but with different indices:

Y[1], Y[2], Y[3], Y[4], Y[5], …

Index order: 1, 2, 3, 4, 5, …

Observe that the each index denotes the position of the value in the array. For example,

Y[1] = 300 (the first number)
Y[3] = 307 (the third value)


Profiles


If we sampled a signal at equally spaced points in time, the result would be a sequence of numbers that could be represented as a function of an index variable. However, a computer with a busy operating system does not guarantee a precise sampling rate. Thus, Artisan records both the sensor data and the time data.

Profiles in Artisan are composed of two signals, sensor information and time information. When recording a roast, sensor readings are taken approximately at every sample interval.  That is, at every sample interval, data is brought from the device to the computer and the time is recorded. For example, with a sampling rate of three seconds we would get

Y = 300, 304, 307, 308, 308, …
t = 0.0, 3.01, 6.02, 9.02, 12.03, 15.03, …  

The values of both signals Y and t can again be represented by indices. For example:

Y[2] = 304t[2] = 3.01

This means that at index 2, the temperature Y=304 was acquired when time t=3.01. If several sensor readings were taken simultaneously, all the signals could be accessible by indices.


Relative Indices


Indices could also be counted from right to left instead of from the left to right. For example, suppose we were looking at the end of a profile

Y = …, 400, 402, 404, 405
t = …, 600.1, 603.1, 606.2, 609.2

with 405 the last Y value and 609.2 the last t value.

Index order: …, 4, 3, 2, 1

So we have
Y[1] = 405, Y[2] = 404, ... 
t[1] = 609.2, t[2] = 606.2, ...

The left or the right end of a signal could be used as absolute points of reference as they cannot be moved. They mark the first and the last point of the signal respectively. But could a middle point be used as a reference instead? Yes. The index reference to other points would be relative to the present point. That is what Artisan uses in the formulas. Other values can be accessed with a + or sign and a number relative to the present point. For example, supposed we were at Y=360, and t=546.2 in the following signal.

Y = …, 355, 357, 360, 362, 365 …
t = …, 540.1, 543.1, 546.2, 549.2, 562.3, …

Index order:  …, -2, -1, 0, +1, +2 , …

0 : present time, +1 : future index, -1 : past index

Y[0]  = 360   and   t[0]  = 546.2
Y[+1] = 362   and   t[+1] = 549.2
Y[-1] = 357   and   t[-1] = 543.1

In Artisan, you access values relative to your present location with the format [<signnumber>] like [-3] (no space in between). For example:

Y     : present value, equivalent to Y[+0] or Y[-0]
Y[-1] : previous value
Y[-2] : second previous value
Y[+1] : next future value


Formula Evaluation


Artisan evaluates a formula for every index of the signal. For example:

A = 1, 2, 2, 3 
B = 2, 3, 4, 1
A+B = 1+2, 2+3, 2+4, 3+1 = 3, 5, 6, 4
A*B = 1*2, 2*3, 2*4, 3*1 = 2, 6, 8, 3


Time Shift


The output of a shifted signal is a signal moved in time.  ”“ shifts the signal to the right while “+” shifts the signal to the left by the number that follows.  For example:

A     = 1, 2, 3, 4, 5, 6   (original signal)
A[-1] = 1, 1, 2, 3, 4, 5   (shifts one to the right)
A[-2] = 1, 1, 1, 2, 3, 4   (shifts two to the right)
A[+1] = 2, 3, 4, 5, 6, 6   (shifts one to the left)
A[+2] = 3, 4, 5, 6, 6, 6   (shifts 2 to the left)

Operations on Symbols and Time Shifts


When applying operations on formulas with shifting, the operation process is kept the same but the shifted signal uses different values. 

A     = [1,2,2,3,4]
A[-1] = [1,1,2,2,3] 
so we have

(A - A[-1])/2 
   = ([1-1, 2-1, 2-2, 3-2, 4-3])/2 
   = ([0, 1, 0, 1, 1])/2
   = [0, 0.5, 0, 0.5, 0.5]

Lets take a closer look at [-1] to see what it does. The [-1] part picks the value of the previous index (left) and puts it in the current (present) index of the new signal through all indexes. For example, the value 3 in A is replaced by the value on its left, which is 2 (index #4). It used the previous index (-1) value. If you look at the whole signal, it can be seen that the signal has moved to the right (shifted). This is also interpreted as delaying the signal since the new signal with [-1] is a copy of the original signal A but delayed in time (moved horizontally). When plotting in Artisan, the shape of A[-1] will appear in time after the shape of A. Thus, moving the signal to the right is delaying it as it would appear later. And moving the signal to the left is advancing it as it would appear earlier.


Symbolic Variables


Devices Variables


All device variables start with a Y plus a number (like Y1). The number following Y are not indices.   Indices are always enclosed in brackets.  Each device in Artisan is conceptually assumed to create two signals (as most temperature meters support dual probes). There is one device considered the main device and a list of extra devices, each contributing two signals. For example, Y1 and Y2 would come from the same (main) device. Y3 and Y4 would come from the another additional device.

Y1 : ET (environmental temperature)
Y2 : BT (bean temperature)
Y3 : Extra device #1 thermocouple A
Y4 : Extra device #1 thermocouple B
Y5 : Extra device #2 thermocouple A
...


Time Access


The variable t holds the sample times. During a roast, CHARGE sets the relative time to zero. Thus, the times shown after CHARGE are relative to CHARGE. The variable t is affected by this. The variable t holds the relative time w.r.t. CHARGE if is registered or the begin of the recording. Thus a formula P(t) = t will always represent a straight line through the origin (0,0) of the graph.


Background Curves


The symbolic variables to access the background curves are as follows.
B1 = ET Background
B2 = BT background
B3 = Extra 1 background A
B4 = Extra 1 background B
...

Plotter Results


The plotter in Artisan supports 9 plotter definitions generating up to 9 extra plots.

P1,..,P9

The variables P1,..,P9 represent the results from plot #1,..,#9. You can perform calculations in a later plot on variables of an earlier plot. That way, the plot variables P1,..,P9 allow the cascading or intermediate results. For example, plot #3 can refer to the results of plot 1 using the variable P1.

F1,..,F9 

F1 refers to the previous result of the actual formula to realize a feedback loop. This is useful in filter designs. Similarly, F2 refers to the second previous result etc.


Plotter Operators and Commands


# operator

# used as a first character in a plot window, will comment out the formula such that it will not be evaluated. As plotter formulas are kept persistent, this commenting feature can be used to keep currently unused formulas.

$ operator

$ used as a first character in a plot window,  will hide the curve on the screen as will do the # operator, but in contrast will still evaluate the rest formula. This is useful for calculating intermediate results to be accessed using the plot variables Pn variable and thus reused in later plots.

annotate(<text>, <time>, <temperature>, <fontsize>)

The annotate command allows to place an annotation label at a certain time/temperature coordinate. For example, annotate(Hello, 02:00, 300, 20) places the label "Hello" in size 20pt at minute 2 and 300 degrees.


Putting Theory into Practice


Rate of Change


A derivative of a function is the rate of change of a function, often called rate of rise (RoR) for the corresponding (mostly positive) signal in coffee roasting. In Artisan, the rate of change of a profile is mathematically speaking the derivative of the signal Y with respect to the signal t.

First Derivative (Delta-Y/Delta-t):  dY/dt = (Y-Y[-1])/(t-t[-1])

This first derivative shows how fast the curve is changing. Thus, it is an indicator of the physical speed of the roast. When the derivative is positive, the curve is sloping upwards.

The following formulas calculates the rate of change of BT:
(Y2-Y2[-1])/(t-t[-1])       (degrees per second)
((Y2-Y2[-1])/(t-t[-1]))*60  (degrees per minute)

60 is the multiplication factor that converts to degrees per minutes from degrees per second.  Artisan uses degrees per minute as unit per default. The rate of change of any variable can be computed by substituting the variable Y2 to any other variable in the formula above.

If we put the formula ((Y2-Y2[-1])/(t-t[-1]))*60 into the plotter P1 field


and press the "Plot" button, we get a new black line drawn on the bottom of our profile.



Hm. Not what we expected. Instead we expected a new line drawn along the dotted blue line that represents the DeltaBT (also representing the RoR of the BT signal). The reason for the lower signal is that the plotter plots on the temperature scale shown on the left and not on the RoR scale shown on the right. This can be easy fixed by computing a corresponding scaling factor k to k=300/40=7.5 (range of the temperature scale divided by the range of the delta scale) and adding it to our formula. With this update formula ((Y2-Y2[-1])/(t-t[-1]))*60*7.5 we get the expected result, thus slightly more noisy. If the minimum of either the temperature or RoR axis is set to anything different from 0, also an offset o has to be added. Two variables k and o are predefined to compute to the correct scaling factors needed based on both y-axis limits. Therefore, we can also use the formula ((Y2-Y2[-1])/(t-t[-1]))*60*k+o to get the very same result with the advantage that we do not have to change the formula on changes of the axis limits.






A smother signal results from increasing the offset between the two readings defining the rate-of-change function. Increasing that offset by two as in  ((Y2-Y2[-3])/(t-t[-3]))*60*k we get a way smoother signal, thus slightly shifted to the right.






With the formula ((Y2[+2]-Y2[-2])/(t[+2]-t[-2]))*60*k we can compute an even smoother rate-of-change signal without a significant shift as seen in the following.





Filtering and Smoothing Signals


Once you have a rate of change in a plot, say in P1, you can smooth it further by making it an input to another plot. For example, let's load our profile again and start the plotter with our initial raw and noisy rate-of-change formula.

P1: ((Y2-Y2[-1])/(t-t[-1]))*60*k

We want to smooth it in plot P3 with a running average filter. Averaging is often applied if data fluctuates too much and must be smoothed out to allow for a proper interpretation. Below we use a “sliding window” of 6 samples. That is, every output requires 6 inputs.  The higher the number of  inputs the smoother.

P3: (P1+P1[-1]+P1[-2]+P1[-3]+P1[-4]+P1[-5])/6

The length of the averaging interval has an effect on the output too.  For example, P4 is an alternative to P3 with a larger time shift:

P4:  (P1 + P1[-9])/2

To formulate an IIR filters, one can use the variable feedback F(1-9) as in the definition of P5:

P5: F1+(P1-F1)/(15/(t-t[-1]))

The above IIR low-pass filter follows the general format

Output = 
    feedback + (Signal – Feedback)/(SmoothingFactor/time-delta)

The higher the smoothing factor (15 in P5 above), the smoother the output signal becomes and the shift in time increases. An averaging filter shifts in time by (length-1)/2. This can be corrected by making its output the input of another plotter window and shifting it back. For example, P3 can be shifted by the following factor sf as in P6:

sf = (6-1)/2  = 2.5 
P6 = P3[+3]  
(we take the offset to 3 as we can shift only by whole integers)

The picture below shows the final plotter configuration as discussed above together with a zoomed version of the resulting graph that also shows the labels placed via the annotation feature.



The resulting graph shows the plots P1(black), P3 (green), P4 (blue), P5 (red) and P6 (pink).





The $ operator could now be used to hide intermediate curves. This method of computation is called cascading, where the output of one filter becomes the input of another. For example, we could put a $ character as the first symbol in P1 and P3 to hide those, since we want only to see the filtered and time corrected version P6.

Another option to smooth a signal would be to apply a FIR filter. So to filter a signal in P1, one could use formulas like the following ones:


  • FIR filter {7,24,34,24,7}/96

(7*P1[-2]+24*P1[-1]+34*P1+24*P1[+1]+7*P1[+3])/96


  • FIR filter {1,6,12,14,12,6,1}/52

(P1[-3]+6*P1[-2]+12*P1[-1]+14*P1+12*P1[+1] 
   +6*P1[+2]+1*P1[+3])/52


Notes

  • the curves drawn by the plotter are only temporal, thus closing the plotter dialog removes the plots too
  • the first two plots, P1 and P2, can be send to the background profile by pressing the "Background" button. This replaces the background ET and BT curves, and thus survive the closing of the plotter
  • the first two plots, P1 and P2, can be turned also into an additional extra Virtual Device entry by pressing the "ET/BT" button. That way, the generated data can be saved along with the current profile. If no profile is loaded, the P1 and P2 plots are used to populate the ET and BT curves of a new profile that can be stored as any regular profile to the file system.


Symbolic Assignments on Logging


The new math formula mechanism is also available during logging to be used as symbolic formulas in device definitions (menu Config >> Devices)
  • time-shifting to the future is obviously not available here
  • feedback via the F variables is not available
  • an additional variable x holds the last reading of the corresponding device channel

Deviation from Background


One application is to add two LCDs via a new extra Virtual Device that computes the temperature deltas of BT and ET w.r.t. the current background curve and its temporal alignment via two symbolic assignments Y1-B1 and Y2-B2. Here we deactivate the drawing of the curves and only activate the LCDs.


Note that the variables Y1 an Y2 are assigned to the raw values retrieved in the actual sampling interval. So any computations that has to be applied to those raw values has to be applied here too.


Rate of Change for Extra Device Curves


Artisan provides a mechanism for the rate-of-change (Rate-of-Rise, RoR) computations only for the main temperature channels as DeltaET and DeltaBT (menu Tools >> Extras).

However, a virtual device with a symbolic assignment formula making use of some time shifting can be used to define such signals. Assuming a first extra device that provides data in variable Y3, we can define a second virtual device with the following definitions

(Y3-Y3[-1])/(t-t[-1])*60
(Y3-Y3[-1])/(t-t[-1])*60*k + o

Note that the first definition does not scale to the temperature range using this k factor such that this signal holds the correct values (in C/min) to be displayed on an LCDs only. The second channel holds the scaled values (in C/min*k), with offset o, to be drawn correctly w.r.t. the temperature axis.




Calculating the rate-of-change will result potentially in a quite noisy signal dramatically amplifying any noise of the underlying signal. One way to improve on this is to increase the delta span, that is the range over which we compute the rate-of-changes. Above it is just a span of one taking exactly two readings in a distance of one sampling interval. The formula below takes a span of 3.

(Y3-Y3[-3])/(t-t[-3])*60
(Y3-Y3[-3])/(t-t[-3])*60*k + o

Of course we can also integrate some smoothing in the delta computation like in the following.

(((Y3-Y3[-1])/(t-t[-1])+(Y3[-1]-Y3[-2])/(t[-1]-t[-2])
   +(Y3[-2]-Y3[-3])/(t[-2]-t[-3]))*60)/3
(((Y3-Y3[-1])/(t-t[-1])+(Y3[-1]-Y3[-2])/(t[-1]-t[-2])
   +(Y3[-2]-Y3[-3])/(t[-2]-t[-3]))*60*k)/3 + o

That smoothing and the extended delta span idea can of course be combined.

2 comments:

  1. Would love to see some 'plug and play' formulas for those of us who get lost in the math. Like "plug in this formula to see a curve that represents [fill in the blank]".

    ReplyDelete