but I am not sure about the whole idea behind it though it resembles the derivative idea of Rob.
Also, how to determine the dt (time between measurements) given that I need the fastest response possible. Should I eliminate the "wait" statement and make dt just = millis() - last_millis() ?
Kp, Ki and Kd are tuning constants and that, I suspect, is more a matter of trial and errors.
How could I apply it to my problem (slow decay curve)?
robtillaart:
Yes, please post the code, you want to share, I'll comment (maybe others too
Here it is (I wanted to triple check it and test it before posting):
/* Breath Sensor and Moving average
By Luca Brigatti
Adapted from an algorithm by robtillaart
Calculates the derivative (difference) between two subsequent sensor readings
and calculates a moving average of the derivative signal
Optimized for a running average of 2^n elements ( 2, 4, 8 , 16 etc.)
Sensor used for testing is breath sensor by Modern Devices:
http://shop.moderndevice.com/products/wind-sensor
*/
const boolean DEBUG = true ; // true, to print data on the serial port. false to do something else
const int MAVSIZE = 16 ; // Elements (size) of the moving average
// Signal variables
int oldval ;
int newval ;
int derivative ; // will hold oldval - newval
// Moving Average variables
int MavValues [MAVSIZE] ; // Holds the last MAVSIZE values of the moving average
int mavtot = 0 ; // total of the moving average (before the division)
int mav = 0 ; // Moving average = mavtot / MAVSIZE
int mavindex = 0 ; // next value to use for the moving average
int sizemask ; // Mask to use for the index counter
void setup() {
Serial.begin(9600); // <--- Low speed: easier to collect data from the serial port, High Speed: smoother curve and tighter moving average
sizemask = MAVSIZE -1 ; // If MAVSIZE is B10000 (16), sizemask is B01111 (15)
setupmav (MAVSIZE) ; // Initialize the Moving Average with the first MAVSIZE values
if (DEBUG) Serial.println("Square breathing 0.1");
}
void loop() {
newval = analogRead(A2); // measure the wind
derivative = (newval - oldval); // first derivative of the signal
updatemav (MAVSIZE) ; // Updates the moving average, a new mav value is assigned.
if (DEBUG) serialdump() ;
else ; // Do something with mav , the moving average of the derivative or newval, the actual reading
oldval = newval ; // Remember previous value for next derivative
}
void serialdump () {
// OUTPUT CAN BE COPIED TO EXCEL
// Put here what you want to print
Serial.print(millis());
Serial.print(";");
Serial.print(oldval);
Serial.print(";");
Serial.print(newval);
Serial.print(";");
Serial.print(derivative);
Serial.print(";");
Serial.print(mav);
Serial.print(";");
// SQUARE OUTPUT (% as example)
if (derivative > 0 && newval > 40) Serial.println(100);
else Serial.println(0);
}
void setupmav(int mavsize) {
mav = 0 ;
mavtot = 0 ;
oldval = analogRead(A2) ; // Starting point, change with Pin actually used
for (int i = 0 ; i < mavsize; i++ ) {
newval = analogRead(A2); // measure the wind, change with pin used
derivative = (newval - oldval); // first derivative of the signal
MavValues [i] = derivative ; // Remember this value
mavtot += derivative; // Update the running total
oldval = newval ; // New Value is now old value
}
mav = mavtot / mavsize ; // First moving average
}
void updatemav (int mavsize) {
mavtot += derivative ; // Add latest value
mavtot -= MavValues [mavindex] ; // Subtract oldest value in the array
MavValues [mavindex] = derivative ; // Newest value in array, replaces oldest value
mav = mavtot / mavsize ; // New moving average
mavindex++ ; // Point to the next element in the array, now the oldest element.
mavindex &= sizemask ; // Circle back to 0 if reached size of the moving average (e.g. 10000 & 01111 (16 & 15) = 00000 (0) )
}
My algorithm is not as abstracted and nearly as useful as your library but it is optimized for a situation like the one at hand:
Putting the initial loading of the moving average array in setup and the main moving average calculation in the loop function eliminates the need for an "if" statement in the main loop to keep the count on how many values are in the array.
Also, this algorithm uses only integer math.
Finally
mavindex++ ; // Point to the next element in the array, now the oldest element.
mavindex &= sizemask ; // Circle back to 0 if reached size of the moving average (e.g. 10000 & 01111 (16 & 15) = 00000 (0) )
Replaces your modulus and my
if (mavindex=16) mavindex = 15 ;
As I feel it may be faster than both (but I really don't know).
dhenry:
A fast moving average would be something like this:
mav = (1-alpha)*mav + alpha * new_data;
where 0<alpha<1 as a weight for new observations. The smaller alpha is, the longer memory mav has of historical data.
If you pick alpha to be 1/2^n, the math gets considerably easier as you can simply shift. For example, for alpha = 1/16:
mav = mav + (new_data - mav) / 16.
essentially, 1 subtraction, 1 shift, 1 increment.
That's very elegant.
It eliminates the need for an array to keep track of the last "n" numbers (and the ++ and &= operations on the counter) and eliminates the initializing function altogether.
You can easily test if this works (depending on the datatypes you can have certain side effects !)
faster? think yes you need to measure 10000 times or so to see diffs.
You can easily test if this works (depending on the datatypes you can have certain side effects !)
faster? think yes you need to measure 10000 times or so to see diffs.
So would this in effect be acting as a low pass filter on the signal? If so how is the 16 related to the 'corner frequency of the effective filtering action, or is that more a function of how often the statement is executed per time unit?
retrolefty:
So would this in effect be acting as a low pass filter on the signal? If so how is the 16 related to the 'corner frequency of the effective filtering action, or is that more a function of how often the statement is executed per time unit?
Ah... I yield to dhenry and robtillaart (not that I would not know the answer....)
Before you get too excited about this approach, it has limitations:
It requires that the new_data to be greater than 1/16th. Otherwise, the moving average will degenerate to 0. One way to avoid this is to use fixed point math. For example, use the last two four (binary) digits as "decimal points":
new_data=analogRead(MY_CHANNEL) << 4; //read the data
mav += ((signed short) (new_data - mav)) >> 4; //compute moving average
You just need to recognize that mav is 16 times bigger.
You can also multiply the reading by 100 (for example) to obtain two decimal points.
The key is to avoid degeneration with small numbers.
dhenry:
There are two issues: signed vs. unsigned, and integer math.
Which one are you referring to?
The value I am averaging is a signed integer which varies roughly between -5 and +5 (but in theory could be more) and the 16 periods moving average (in my tests) is very often 0 with some peaks between -2 and +2.
Yes, it will have issues with long memory: the smoothed signal would accurate reproduce the peaks / valleys.
Two solutions:
use shorter memory: rather than 1/16th, use 1/8th, or 1/4th. The data generated this way will be more volatile.
scale up observations: rather than having data series in -5 - +5 range, scale them up by 10 (or 16, using left shifts) for example. The resulting moving average will be correspondingly scaled up.
Just wanted to report on my results, for whomever may be interested.
I applied a simple "decision" to the 16 period moving average:
IF moving average > 0 I am blowing on the sensor
IF moving average < 0 I am NOT blowing on the sensor
with the corollary that:
IF moving average == 0 I am doing whatever I was doing before // No need to actually state
This is the result in graphic form:
The blue line is the sensor data, the red line is the moving average (multiplied by 100 for ease of plotting) and the green line is the result of the decision: 100 I am blowing, 0 I am not blowing.
It's getting accurate, though there are some false positives and false negatives now and then and, in this very example, there are 86 ms of delay from the beginning of the decay (when I stop blowing) to the response of the moving average (when the software realized that I stopped blowing).
Unfortunately, there isn't a whole lot one can do about it: the process we implemented here is an integration process. The longer the memory, the better it is at eliminating noises, and the longer the delay is.
You will just have to be careful in selecting the memory.
Yes, I'm coming more or less to the same conclusion.
Since the signal is quite clean to begin with, the next thing will be to see if an 8 period average gives a shorter delay with an acceptable low number of false positives/negatives.
Interesting (and very frustrating) the portion of sketch that sends the MIDI signal to the serial would not work unless i put a delay (10) ; In the main loop.
But it will work without the delay, if I Serial.print() ; the data on the serial monitor for debugging.