tool/mbed/mbed-sdk/libraries/dsp/cmsis_dsp/FilteringFunctions/arm_biquad_cascade_df2T_f32.c

   1 /* ----------------------------------------------------------------------
   2 * Copyright (C) 2010-2013 ARM Limited. All rights reserved.
   3 *
   4 * $Date:        17. January 2013
   5 *
   6 * Project:          CMSIS DSP Library
   7 * Title:            arm_biquad_cascade_df2T_f32.c
   8 *
   9 * Description:  Processing function for the floating-point transposed
  10 *               direct form II Biquad cascade filter.
  11 *
  12 * Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
  13 *
  14 * Redistribution and use in source and binary forms, with or without
  15 * modification, are permitted provided that the following conditions
  16 * are met:
  17 *   - Redistributions of source code must retain the above copyright
  18 *     notice, this list of conditions and the following disclaimer.
  19 *   - Redistributions in binary form must reproduce the above copyright
  20 *     notice, this list of conditions and the following disclaimer in
  21 *     the documentation and/or other materials provided with the
  22 *     distribution.
  23 *   - Neither the name of ARM LIMITED nor the names of its contributors
  24 *     may be used to endorse or promote products derived from this
  25 *     software without specific prior written permission.
  26 *
  27 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  28 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  29 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  30 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  31 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  32 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  33 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  34 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  35 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  36 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  37 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  38 * POSSIBILITY OF SUCH DAMAGE.
  39 * -------------------------------------------------------------------- */
  40
  41 #include "arm_math.h"
  42
  43 /**
  44 * @ingroup groupFilters
  45 */
  46
  47 /**
  48 * @defgroup BiquadCascadeDF2T Biquad Cascade IIR Filters Using a Direct Form II Transposed Structure
  49 *
  50 * This set of functions implements arbitrary order recursive (IIR) filters using a transposed direct form II structure.
  51 * The filters are implemented as a cascade of second order Biquad sections.
  52 * These functions provide a slight memory savings as compared to the direct form I Biquad filter functions.
  53 * Only floating-point data is supported.
  54 *
  55 * This function operate on blocks of input and output data and each call to the function
  56 * processes <code>blockSize</code> samples through the filter.
  57 * <code>pSrc</code> points to the array of input data and
  58 * <code>pDst</code> points to the array of output data.
  59 * Both arrays contain <code>blockSize</code> values.
  60 *
  61 * \par Algorithm
  62 * Each Biquad stage implements a second order filter using the difference equation:
  63 * <pre>
  64 *    y[n] = b0 * x[n] + d1
  65 *    d1 = b1 * x[n] + a1 * y[n] + d2
  66 *    d2 = b2 * x[n] + a2 * y[n]
  67 * </pre>
  68 * where d1 and d2 represent the two state values.
  69 *
  70 * \par
  71 * A Biquad filter using a transposed Direct Form II structure is shown below.
  72 * \image html BiquadDF2Transposed.gif "Single transposed Direct Form II Biquad"
  73 * Coefficients <code>b0, b1, and b2 </code> multiply the input signal <code>x[n]</code> and are referred to as the feedforward coefficients.
  74 * Coefficients <code>a1</code> and <code>a2</code> multiply the output signal <code>y[n]</code> and are referred to as the feedback coefficients.
  75 * Pay careful attention to the sign of the feedback coefficients.
  76 * Some design tools flip the sign of the feedback coefficients:
  77 * <pre>
  78 *    y[n] = b0 * x[n] + d1;
  79 *    d1 = b1 * x[n] - a1 * y[n] + d2;
  80 *    d2 = b2 * x[n] - a2 * y[n];
  81 * </pre>
  82 * In this case the feedback coefficients <code>a1</code> and <code>a2</code> must be negated when used with the CMSIS DSP Library.
  83 *
  84 * \par
  85 * Higher order filters are realized as a cascade of second order sections.
  86 * <code>numStages</code> refers to the number of second order stages used.
  87 * For example, an 8th order filter would be realized with <code>numStages=4</code> second order stages.
  88 * A 9th order filter would be realized with <code>numStages=5</code> second order stages with the
  89 * coefficients for one of the stages configured as a first order filter (<code>b2=0</code> and <code>a2=0</code>).
  90 *
  91 * \par
  92 * <code>pState</code> points to the state variable array.
  93 * Each Biquad stage has 2 state variables <code>d1</code> and <code>d2</code>.
  94 * The state variables are arranged in the <code>pState</code> array as:
  95 * <pre>
  96 *     {d11, d12, d21, d22, ...}
  97 * </pre>
  98 * where <code>d1x</code> refers to the state variables for the first Biquad and
  99 * <code>d2x</code> refers to the state variables for the second Biquad.
 100 * The state array has a total length of <code>2*numStages</code> values.
 101 * The state variables are updated after each block of data is processed; the coefficients are untouched.
 102 *
 103 * \par
 104 * The CMSIS library contains Biquad filters in both Direct Form I and transposed Direct Form II.
 105 * The advantage of the Direct Form I structure is that it is numerically more robust for fixed-point data types.
 106 * That is why the Direct Form I structure supports Q15 and Q31 data types.
 107 * The transposed Direct Form II structure, on the other hand, requires a wide dynamic range for the state variables <code>d1</code> and <code>d2</code>.
 108 * Because of this, the CMSIS library only has a floating-point version of the Direct Form II Biquad.
 109 * The advantage of the Direct Form II Biquad is that it requires half the number of state variables, 2 rather than 4, per Biquad stage.
 110 *
 111 * \par Instance Structure
 112 * The coefficients and state variables for a filter are stored together in an instance data structure.
 113 * A separate instance structure must be defined for each filter.
 114 * Coefficient arrays may be shared among several instances while state variable arrays cannot be shared.
 115 *
 116 * \par Init Functions
 117 * There is also an associated initialization function.
 118 * The initialization function performs following operations:
 119 * - Sets the values of the internal structure fields.
 120 * - Zeros out the values in the state buffer.
 121 * To do this manually without calling the init function, assign the follow subfields of the instance structure:
 122 * numStages, pCoeffs, pState. Also set all of the values in pState to zero.
 123 *
 124 * \par
 125 * Use of the initialization function is optional.
 126 * However, if the initialization function is used, then the instance structure cannot be placed into a const data section.
 127 * To place an instance structure into a const data section, the instance structure must be manually initialized.
 128 * Set the values in the state buffer to zeros before static initialization.
 129 * For example, to statically initialize the instance structure use
 130 * <pre>
 131 *     arm_biquad_cascade_df2T_instance_f32 S1 = {numStages, pState, pCoeffs};
 132 * </pre>
 133 * where <code>numStages</code> is the number of Biquad stages in the filter; <code>pState</code> is the address of the state buffer.
 134 * <code>pCoeffs</code> is the address of the coefficient buffer;
 135 *
 136 */
 137
 138 /**
 139 * @addtogroup BiquadCascadeDF2T
 140 * @{
 141 */
 142
 143 /**
 144 * @brief Processing function for the floating-point transposed direct form II Biquad cascade filter.
 145 * @param[in]  *S        points to an instance of the filter data structure.
 146 * @param[in]  *pSrc     points to the block of input data.
 147 * @param[out] *pDst     points to the block of output data
 148 * @param[in]  blockSize number of samples to process.
 149 * @return none.
 150 */
 151
 152
 153 LOW_OPTIMIZATION_ENTER
 154 void arm_biquad_cascade_df2T_f32(
 155 const arm_biquad_cascade_df2T_instance_f32 * S,
 156 float32_t * pSrc,
 157 float32_t * pDst,
 158 uint32_t blockSize)
 159 {
 160
 161    float32_t *pIn = pSrc;                         /*  source pointer            */
 162    float32_t *pOut = pDst;                        /*  destination pointer       */
 163    float32_t *pState = S->pState;                 /*  State pointer             */
 164    float32_t *pCoeffs = S->pCoeffs;               /*  coefficient pointer       */
 165    float32_t acc1;                                /*  accumulator               */
 166    float32_t b0, b1, b2, a1, a2;                  /*  Filter coefficients       */
 167    float32_t Xn1;                                 /*  temporary input           */
 168    float32_t d1, d2;                              /*  state variables           */
 169    uint32_t sample, stage = S->numStages;         /*  loop counters             */
 170
 171 #ifndef ARM_MATH_CM0_FAMILY_FAMILY
 172
 173    float32_t Xn2, Xn3, Xn4;                       /*  Input State variables     */
 174    float32_t acc2, acc3, acc4;                            /*  accumulator               */
 175
 176
 177    float32_t p0, p1, p2, p3, p4, A1;
 178
 179    /* Run the below code for Cortex-M4 and Cortex-M3 */
 180    do
 181    {
 182       /* Reading the coefficients */
 183       b0 = *pCoeffs++;
 184       b1 = *pCoeffs++;
 185       b2 = *pCoeffs++;
 186       a1 = *pCoeffs++;
 187       a2 = *pCoeffs++;
 188
 189
 190       /*Reading the state values */
 191       d1 = pState[0];
 192       d2 = pState[1];
 193
 194       /* Apply loop unrolling and compute 4 output values simultaneously. */
 195       sample = blockSize >> 2u;
 196
 197       /* First part of the processing with loop unrolling.  Compute 4 outputs at a time.
 198    ** a second loop below computes the remaining 1 to 3 samples. */
 199       while(sample > 0u) {
 200
 201          /* y[n] = b0 * x[n] + d1 */
 202          /* d1 = b1 * x[n] + a1 * y[n] + d2 */
 203          /* d2 = b2 * x[n] + a2 * y[n] */
 204
 205          /* Read the four inputs */
 206          Xn1 = pIn[0];
 207          Xn2 = pIn[1];
 208          Xn3 = pIn[2];
 209          Xn4 = pIn[3];
 210          pIn += 4;
 211
 212          p0 = b0 * Xn1;
 213          p1 = b1 * Xn1;
 214          acc1 = p0 + d1;
 215          p0 = b0 * Xn2;
 216          p3 = a1 * acc1;
 217          p2 = b2 * Xn1;
 218          A1 = p1 + p3;
 219          p4 = a2 * acc1;
 220          d1 = A1 + d2;
 221          d2 = p2 + p4;
 222
 223          p1 = b1 * Xn2;
 224          acc2 = p0 + d1;
 225          p0 = b0 * Xn3;
 226          p3 = a1 * acc2;
 227          p2 = b2 * Xn2;
 228          A1 = p1 + p3;
 229          p4 = a2 * acc2;
 230          d1 = A1 + d2;
 231          d2 = p2 + p4;
 232
 233          p1 = b1 * Xn3;
 234          acc3 = p0 + d1;
 235          p0 = b0 * Xn4;
 236          p3 = a1 * acc3;
 237          p2 = b2 * Xn3;
 238          A1 = p1 + p3;
 239          p4 = a2 * acc3;
 240          d1 = A1 + d2;
 241          d2 = p2 + p4;
 242
 243          acc4 = p0 + d1;
 244          p1 = b1 * Xn4;
 245          p3 = a1 * acc4;
 246          p2 = b2 * Xn4;
 247          A1 = p1 + p3;
 248          p4 = a2 * acc4;
 249          d1 = A1 + d2;
 250          d2 = p2 + p4;
 251
 252          pOut[0] = acc1;
 253          pOut[1] = acc2;
 254          pOut[2] = acc3;
 255          pOut[3] = acc4;
 256                                  pOut += 4;
 257
 258          sample--;
 259       }
 260
 261       sample = blockSize & 0x3u;
 262       while(sample > 0u) {
 263          Xn1 = *pIn++;
 264
 265          p0 = b0 * Xn1;
 266          p1 = b1 * Xn1;
 267          acc1 = p0 + d1;
 268          p3 = a1 * acc1;
 269          p2 = b2 * Xn1;
 270          A1 = p1 + p3;
 271          p4 = a2 * acc1;
 272          d1 = A1 + d2;
 273          d2 = p2 + p4;
 274
 275          *pOut++ = acc1;
 276
 277          sample--;
 278       }
 279
 280       /* Store the updated state variables back into the state array */
 281       *pState++ = d1;
 282       *pState++ = d2;
 283
 284       /* The current stage input is given as the output to the next stage */
 285       pIn = pDst;
 286
 287       /*Reset the output working pointer */
 288       pOut = pDst;
 289
 290       /* decrement the loop counter */
 291       stage--;
 292
 293    } while(stage > 0u);
 294
 295 #else
 296
 297    /* Run the below code for Cortex-M0 */
 298
 299    do
 300    {
 301       /* Reading the coefficients */
 302       b0 = *pCoeffs++;
 303       b1 = *pCoeffs++;
 304       b2 = *pCoeffs++;
 305       a1 = *pCoeffs++;
 306       a2 = *pCoeffs++;
 307
 308       /*Reading the state values */
 309       d1 = pState[0];
 310       d2 = pState[1];
 311
 312
 313       sample = blockSize;
 314
 315       while(sample > 0u)
 316       {
 317          /* Read the input */
 318          Xn1 = *pIn++;
 319
 320          /* y[n] = b0 * x[n] + d1 */
 321          acc1 = (b0 * Xn1) + d1;
 322
 323          /* Store the result in the accumulator in the destination buffer. */
 324          *pOut++ = acc1;
 325
 326          /* Every time after the output is computed state should be updated. */
 327          /* d1 = b1 * x[n] + a1 * y[n] + d2 */
 328          d1 = ((b1 * Xn1) + (a1 * acc1)) + d2;
 329
 330          /* d2 = b2 * x[n] + a2 * y[n] */
 331          d2 = (b2 * Xn1) + (a2 * acc1);
 332
 333          /* decrement the loop counter */
 334          sample--;
 335       }
 336
 337       /* Store the updated state variables back into the state array */
 338       *pState++ = d1;
 339       *pState++ = d2;
 340
 341       /* The current stage input is given as the output to the next stage */
 342       pIn = pDst;
 343
 344       /*Reset the output working pointer */
 345       pOut = pDst;
 346
 347       /* decrement the loop counter */
 348       stage--;
 349
 350    } while(stage > 0u);
 351
 352 #endif /*  #ifndef ARM_MATH_CM0_FAMILY         */
 353
 354 }
 355 LOW_OPTIMIZATION_EXIT
 356
 357 /**
 358    * @} end of BiquadCascadeDF2T group
 359    */