flight/libraries/math/butterworth.c

   1 /**
   2  ******************************************************************************
   3  * @addtogroup OpenPilot Math Utilities
   4  * @{
   5  * @addtogroup Butterworth low pass filter
   6  * @{
   7  *
   8  * @file       butterworth.c
   9  * @author     The OpenPilot Team, http://www.openpilot.org Copyright (C) 2014.
  10  * @brief      Direct form two of a second order Butterworth low pass filter
  11  *
  12  * @see        The GNU Public License (GPL) Version 3
  13  *
  14  *****************************************************************************/
  15 /*
  16  * This program is free software; you can redistribute it and/or modify
  17  * it under the terms of the GNU General Public License as published by
  18  * the Free Software Foundation; either version 3 of the License, or
  19  * (at your option) any later version.
  20  *
  21  * This program is distributed in the hope that it will be useful, but
  22  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  23  * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  24  * for more details.
  25  *
  26  * You should have received a copy of the GNU General Public License along
  27  * with this program; if not, write to the Free Software Foundation, Inc.,
  28  * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  29  */
  30
  31 #include "openpilot.h"
  32 #include "math.h"
  33 #include "butterworth.h"
  34
  35 /**
  36  * Initialization function for coefficients of a second order Butterworth biquadratic filter in direct from 2.
  37  * Note that b1  = 2 * b0 and b2  = b0 is use here and in the sequel.
  38  * @param[in]  ff Cut-off frequency ratio
  39  * @param[out] filterPtr Pointer to filter coefficients
  40  * @returns Nothing
  41  */
  42 void InitButterWorthDF2Filter(const float ff, struct ButterWorthDF2Filter *filterPtr)
  43 {
  44     const float ita = 1.0f / tanf(M_PI_F * ff);
  45     const float b0  = 1.0f / (1.0f + M_SQRT2_F * ita + ita * ita);
  46     const float a1  = 2.0f * b0 * (ita * ita - 1.0f);
  47     const float a2  = -b0 * (1.0f - M_SQRT2_F * ita + ita * ita);
  48
  49     filterPtr->b0 = b0;
  50     filterPtr->a1 = a1;
  51     filterPtr->a2 = a2;
  52 }
  53
  54
  55 /**
  56  * Initialization function for intermediate values of a second order Butterworth biquadratic filter in direct from 2.
  57  * Obtained by solving a linear equation system.
  58  * @param[in]  x0 Prescribed value
  59  * @param[in]  filterPtr Pointer to filter coefficients
  60  * @param[out] wn1Ptr Pointer to first intermediate value
  61  * @param[out] wn2Ptr Pointer to second intermediate value
  62  * @returns Nothing
  63  */
  64 void InitButterWorthDF2Values(const float x0, const struct ButterWorthDF2Filter *filterPtr, float *wn1Ptr, float *wn2Ptr)
  65 {
  66     const float b0   = filterPtr->b0;
  67     const float a1   = filterPtr->a1;
  68     const float a2   = filterPtr->a2;
  69
  70     const float a11  = 2.0f + a1;
  71     const float a12  = 1.0f + a2;
  72     const float a21  = 2.0f + a1 * a1 + a2;
  73     const float a22  = 1.0f + a1 * a2;
  74     const float det  = a11 * a22 - a12 * a21;
  75     const float rhs1 = x0 / b0 - x0;
  76     const float rhs2 = x0 / b0 - x0 + a1 * x0;
  77
  78     *wn1Ptr = (a22 * rhs1 - a12 * rhs2) / det;
  79     *wn2Ptr = (-a21 * rhs1 + a11 * rhs2) / det;
  80 }
  81
  82
  83 /**
  84  * Second order Butterworth biquadratic filter in direct from 2, such that only two values wn1=w[n-1] and wn2=w[n-2] need to be stored.
  85  * Function takes care of updating the values wn1 and wn2.
  86  * @param[in]  xn New raw value
  87  * @param[in]  filterPtr Pointer to filter coefficients
  88  * @param[out] wn1Ptr Pointer to first intermediate value
  89  * @param[out] wn2Ptr Pointer to second intermediate value
  90  * @returns Filtered value
  91  */
  92 float FilterButterWorthDF2(const float xn, const struct ButterWorthDF2Filter *filterPtr, float *wn1Ptr, float *wn2Ptr)
  93 {
  94     const float wn  = xn + filterPtr->a1 * (*wn1Ptr) + filterPtr->a2 * (*wn2Ptr);
  95     const float val = filterPtr->b0 * (wn + 2.0f * (*wn1Ptr) + (*wn2Ptr));
  96
  97     *wn2Ptr = *wn1Ptr;
  98     *wn1Ptr = wn;
  99     return val;
 100 }