|
| 1 | +#include <math.h> |
| 2 | +#include "fast_kalman.h" |
| 3 | + |
| 4 | +fastKalman_t filter[3]; |
| 5 | +filterTypdef filterType; |
| 6 | + |
| 7 | +void initFilter(fastKalman_t *filter, float q, float r, float p, float intialValue) { |
| 8 | + filter->q = q * 0.001f; //add multiplier to make tuning easier |
| 9 | + filter->r = r; //add multiplier to make tuning easier |
| 10 | + filter->p = p; //add multiplier to make tuning easier |
| 11 | + filter->x = intialValue; //set intial value, can be zero if unknown |
| 12 | + filter->lastX = intialValue; //set intial value, can be zero if unknown |
| 13 | + filter->k = 0.0f; //kalman gain, |
| 14 | + filter->gyroDfkfDataPtr = 0; |
| 15 | +} |
| 16 | + |
| 17 | +void fastKalmanInit(float q, float r, float p, float intialValue, filterTypdef type) |
| 18 | +{ |
| 19 | + filterType = type; |
| 20 | + initFilter(&filter[0], q, r, p, intialValue); |
| 21 | + initFilter(&filter[1], q, r, p, intialValue); |
| 22 | + initFilter(&filter[2], q, r, p, intialValue); |
| 23 | +} |
| 24 | + |
| 25 | + |
| 26 | +float noiseEstimate(float data[], uint32_t size) |
| 27 | +{ |
| 28 | + uint32_t i; |
| 29 | + float sum = 0.0, sumOfSquares = 0.0, stdDev; |
| 30 | + for (int i = size+1; i >= 0; i--) |
| 31 | + { |
| 32 | + sum += data[i]; |
| 33 | + sumOfSquares += powf(data[i], 2); |
| 34 | + } |
| 35 | + |
| 36 | + arm_sqrt_f32((sumOfSquares-powf(sum, 2))/(float)size, &stdDev); |
| 37 | + return (stdDev * 0.003f); |
| 38 | +} |
| 39 | + |
| 40 | +float distanceEstimate(float data[], uint32_t size) |
| 41 | +{ |
| 42 | + float largest = -10000.0f; |
| 43 | + float smallest = 10000.0f; |
| 44 | + float range; |
| 45 | + |
| 46 | + //find the range |
| 47 | + for(int x=size;x>=0;x--) |
| 48 | + { |
| 49 | + if(array[x]>largest) |
| 50 | + largest = array[x]; |
| 51 | + |
| 52 | + if(array[x]<smallest) |
| 53 | + smallest = array[x]; |
| 54 | + } |
| 55 | + |
| 56 | + range = largest - smallest; |
| 57 | + |
| 58 | + return (range * 0.001f); |
| 59 | +} |
| 60 | + |
| 61 | +float _fastKalmanUpdate(fastKalman_t *axisFilter) |
| 62 | +{ |
| 63 | + //project the state ahead using acceleration |
| 64 | + axisFilter->x += (axisFilter->x - axisFilter->lastX); |
| 65 | + |
| 66 | + //update last state |
| 67 | + axisFilter->lastX = axisFilter->x; |
| 68 | + |
| 69 | + //prediction update |
| 70 | + axisFilter->p = axisFilter->p + axisFilter->q; |
| 71 | + |
| 72 | + //measurement update |
| 73 | + axisFilter->k = axisFilter->p / (axisFilter->p + axisFilter->r); |
| 74 | + axisFilter->x += axisFilter->k * (input - axisFilter->x); |
| 75 | + axisFilter->p = (1.0f - axisFilter->k) * axisFilter->p; |
| 76 | + return axisFilter->x; |
| 77 | +} |
| 78 | + |
| 79 | +float fastKalmanUpdate(filterAxisTypeDef axis, float input) |
| 80 | +{ |
| 81 | + if (filterType == STD_DEV_ESTIMATION) { |
| 82 | + filter[axis]->r = noiseEstimate(filter[axis]->gyroDfkfData, 6); |
| 83 | + } else if (filterType == DISTANCE_ESTIMATION) { |
| 84 | + filter[axis]->r = distanceEstimate(filter[axis]->gyroDfkfData, 6); |
| 85 | + } |
| 86 | + filter[axis]->gyroDfkfData[filter[axis]->gyroDfkfDataPtr] = input; |
| 87 | + _fastKalmanUpdate(&filter[axis]); |
| 88 | + filter[axis]->gyroDfkfDataPtr++; |
| 89 | + if (filter[axis]->gyroDfkfDataPtr > 5) |
| 90 | + { |
| 91 | + filter[axis]->gyroDfkfDataPtr=0; |
| 92 | + } |
| 93 | + return filter[axis]->x; |
| 94 | +} |
0 commit comments