src/main/flight/lowpass.c

   1 /*
   2  * This file is part of Cleanflight.
   3  *
   4  * Cleanflight is free software: you can redistribute it and/or modify
   5  * it under the terms of the GNU General Public License as published by
   6  * the Free Software Foundation, either version 3 of the License, or
   7  * (at your option) any later version.
   8  *
   9  * Cleanflight is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with Cleanflight.  If not, see <http://www.gnu.org/licenses/>.
  16  */
  17
  18 #include <stdbool.h>
  19 #include <stdint.h>
  20 #include <stdlib.h>
  21 #include <math.h>
  22
  23 #include "flight/lowpass.h"
  24
  25 //#define DEBUG_LOWPASS
  26
  27 void generateLowpassCoeffs2(int16_t freq, lowpass_t *filter)
  28 {
  29     float fixedScaler;
  30     int i;
  31
  32     // generates coefficients for a 2nd-order butterworth low-pass filter
  33     float freqf = (float)freq*0.001f;
  34     float omega = tanf((float)M_PI*freqf/2.0f);
  35     float scaling = 1.0f / (omega*omega + 1.4142136f*omega + 1.0f);
  36
  37
  38 #if defined(UNIT_TEST) && defined(DEBUG_LOWPASS)
  39     printf("lowpass cutoff: %f, omega: %f\n", freqf, omega);
  40 #endif
  41
  42     filter->bf[0] = scaling * omega*omega;
  43     filter->bf[1] = 2.0f * filter->bf[0];
  44     filter->bf[2] = filter->bf[0];
  45
  46     filter->af[0] = 1.0f;
  47     filter->af[1] = scaling * (2.0f * omega*omega - 2.0f);
  48     filter->af[2] = scaling * (omega*omega - 1.4142136f * omega + 1.0f);
  49
  50     // Scale for fixed-point
  51     filter->input_bias = 1500; // Typical servo range is 1500 +/- 500
  52     filter->input_shift = 16;
  53     filter->coeff_shift = 24;
  54     fixedScaler = (float)(1ULL << filter->coeff_shift);
  55     for (i = 0; i < LOWPASS_NUM_COEF; i++) {
  56         filter->a[i] = LPF_ROUND(filter->af[i] * fixedScaler);
  57         filter->b[i] = LPF_ROUND(filter->bf[i] * fixedScaler);
  58 #if defined(UNIT_TEST) && defined(DEBUG_LOWPASS)
  59         printf("(%d) bf: %f af: %f b: %ld a: %ld\n", i,
  60                 filter->bf[i], filter->af[i], filter->b[i], filter->a[i]);
  61 #endif
  62     }
  63
  64     filter->freq = freq;
  65 }
  66
  67 int32_t lowpassFixed(lowpass_t *filter, int32_t in, int16_t freq)
  68 {
  69     int16_t coefIdx;
  70     int64_t out;
  71     int32_t in_s;
  72
  73     // Check to see if cutoff frequency changed
  74     if (freq != filter->freq) {
  75         filter->init = false;
  76     }
  77
  78     // Initialize if needed
  79     if (!filter->init) {
  80         generateLowpassCoeffs2(freq, filter);
  81         for (coefIdx = 0; coefIdx < LOWPASS_NUM_COEF; coefIdx++) {
  82             filter->x[coefIdx] = (in - filter->input_bias) << filter->input_shift;
  83             filter->y[coefIdx] = (in - filter->input_bias) << filter->input_shift;
  84         }
  85         filter->init = true;
  86     }
  87
  88     // Unbias input and scale
  89     in_s = (in - filter->input_bias) << filter->input_shift;
  90
  91     // Delays
  92     for (coefIdx = LOWPASS_NUM_COEF-1; coefIdx > 0; coefIdx--) {
  93         filter->x[coefIdx] = filter->x[coefIdx-1];
  94         filter->y[coefIdx] = filter->y[coefIdx-1];
  95     }
  96     filter->x[0] = in_s;
  97
  98     // Accumulate result
  99     out = filter->x[0] * filter->b[0];
 100     for (coefIdx = 1; coefIdx < LOWPASS_NUM_COEF; coefIdx++) {
 101         out -= filter->y[coefIdx] * filter->a[coefIdx];
 102         out += filter->x[coefIdx] * filter->b[coefIdx];
 103     }
 104
 105     // Scale output by coefficient shift
 106     out >>= filter->coeff_shift;
 107     filter->y[0] = (int32_t)out;
 108
 109     // Scale output by input shift and round
 110     out = (out + (1 << (filter->input_shift-1))) >> filter->input_shift;
 111
 112     // Reapply bias
 113     out += filter->input_bias;
 114
 115     return (int32_t)out;
 116 }
 117