CMSIS DSP Library from CMSIS 2.0. See http://www.onarm.com/cmsis/ for full details

Dependents:   K22F_DSP_Matrix_least_square BNO055-ELEC3810 1BNO055 ECE4180Project--Slave2 ... more

Committer:
simon
Date:
Thu Mar 10 15:07:50 2011 +0000
Revision:
0:1014af42efd9

        

Who changed what in which revision?

UserRevisionLine numberNew contents of line
simon 0:1014af42efd9 1 /* ----------------------------------------------------------------------
simon 0:1014af42efd9 2 * Copyright (C) 2010 ARM Limited. All rights reserved.
simon 0:1014af42efd9 3 *
simon 0:1014af42efd9 4 * $Date: 29. November 2010
simon 0:1014af42efd9 5 * $Revision: V1.0.3
simon 0:1014af42efd9 6 *
simon 0:1014af42efd9 7 * Project: CMSIS DSP Library
simon 0:1014af42efd9 8 * Title: arm_biquad_cascade_df1_fast_q31.c
simon 0:1014af42efd9 9 *
simon 0:1014af42efd9 10 * Description: Processing function for the
simon 0:1014af42efd9 11 * Q31 Fast Biquad cascade DirectFormI(DF1) filter.
simon 0:1014af42efd9 12 *
simon 0:1014af42efd9 13 * Target Processor: Cortex-M4/Cortex-M3
simon 0:1014af42efd9 14 *
simon 0:1014af42efd9 15 * Version 1.0.3 2010/11/29
simon 0:1014af42efd9 16 * Re-organized the CMSIS folders and updated documentation.
simon 0:1014af42efd9 17 *
simon 0:1014af42efd9 18 * Version 1.0.2 2010/11/11
simon 0:1014af42efd9 19 * Documentation updated.
simon 0:1014af42efd9 20 *
simon 0:1014af42efd9 21 * Version 1.0.1 2010/10/05
simon 0:1014af42efd9 22 * Production release and review comments incorporated.
simon 0:1014af42efd9 23 *
simon 0:1014af42efd9 24 * Version 1.0.0 2010/09/20
simon 0:1014af42efd9 25 * Production release and review comments incorporated.
simon 0:1014af42efd9 26 *
simon 0:1014af42efd9 27 * Version 0.0.9 2010/08/27
simon 0:1014af42efd9 28 * Initial version
simon 0:1014af42efd9 29 *
simon 0:1014af42efd9 30 * -------------------------------------------------------------------- */
simon 0:1014af42efd9 31
simon 0:1014af42efd9 32 #include "arm_math.h"
simon 0:1014af42efd9 33
simon 0:1014af42efd9 34 /**
simon 0:1014af42efd9 35 * @ingroup groupFilters
simon 0:1014af42efd9 36 */
simon 0:1014af42efd9 37
simon 0:1014af42efd9 38 /**
simon 0:1014af42efd9 39 * @addtogroup BiquadCascadeDF1
simon 0:1014af42efd9 40 * @{
simon 0:1014af42efd9 41 */
simon 0:1014af42efd9 42
simon 0:1014af42efd9 43 /**
simon 0:1014af42efd9 44 * @details
simon 0:1014af42efd9 45 *
simon 0:1014af42efd9 46 * @param[in] *S points to an instance of the Q31 Biquad cascade structure.
simon 0:1014af42efd9 47 * @param[in] *pSrc points to the block of input data.
simon 0:1014af42efd9 48 * @param[out] *pDst points to the block of output data.
simon 0:1014af42efd9 49 * @param[in] blockSize number of samples to process per call.
simon 0:1014af42efd9 50 * @return none.
simon 0:1014af42efd9 51 *
simon 0:1014af42efd9 52 * <b>Scaling and Overflow Behavior:</b>
simon 0:1014af42efd9 53 * \par
simon 0:1014af42efd9 54 * This function is optimized for speed at the expense of fixed-point precision and overflow protection.
simon 0:1014af42efd9 55 * The result of each 1.31 x 1.31 multiplication is truncated to 2.30 format.
simon 0:1014af42efd9 56 * These intermediate results are added to a 2.30 accumulator.
simon 0:1014af42efd9 57 * Finally, the accumulator is saturated and converted to a 1.31 result.
simon 0:1014af42efd9 58 * The fast version has the same overflow behavior as the standard version and provides less precision since it discards the low 32 bits of each multiplication result.
simon 0:1014af42efd9 59 * In order to avoid overflows completely the input signal must be scaled down by two bits and lie in the range [-0.25 +0.25). Use the intialization function
simon 0:1014af42efd9 60 * arm_biquad_cascade_df1_init_q31() to initialize filter structure.
simon 0:1014af42efd9 61 *
simon 0:1014af42efd9 62 * \par
simon 0:1014af42efd9 63 * Refer to the function <code>arm_biquad_cascade_df1_q31()</code> for a slower implementation of this function which uses 64-bit accumulation to provide higher precision. Both the slow and the fast versions use the same instance structure.
simon 0:1014af42efd9 64 * Use the function <code>arm_biquad_cascade_df1_init_q31()</code> to initialize the filter structure.
simon 0:1014af42efd9 65 */
simon 0:1014af42efd9 66
simon 0:1014af42efd9 67 void arm_biquad_cascade_df1_fast_q31(
simon 0:1014af42efd9 68 const arm_biquad_casd_df1_inst_q31 * S,
simon 0:1014af42efd9 69 q31_t * pSrc,
simon 0:1014af42efd9 70 q31_t * pDst,
simon 0:1014af42efd9 71 uint32_t blockSize)
simon 0:1014af42efd9 72 {
simon 0:1014af42efd9 73 q31_t *pIn = pSrc; /* input pointer initialization */
simon 0:1014af42efd9 74 q31_t *pOut = pDst; /* output pointer initialization */
simon 0:1014af42efd9 75 q31_t *pState = S->pState; /* pState pointer initialization */
simon 0:1014af42efd9 76 q31_t *pCoeffs = S->pCoeffs; /* coeff pointer initialization */
simon 0:1014af42efd9 77 q31_t acc; /* accumulator */
simon 0:1014af42efd9 78 q31_t Xn1, Xn2, Yn1, Yn2; /* Filter state variables */
simon 0:1014af42efd9 79 q31_t b0, b1, b2, a1, a2; /* Filter coefficients */
simon 0:1014af42efd9 80 q31_t Xn; /* temporary input */
simon 0:1014af42efd9 81 int32_t shift = (int32_t) S->postShift + 1; /* Shift to be applied to the output */
simon 0:1014af42efd9 82 uint32_t sample, stage = S->numStages; /* loop counters */
simon 0:1014af42efd9 83
simon 0:1014af42efd9 84
simon 0:1014af42efd9 85 do
simon 0:1014af42efd9 86 {
simon 0:1014af42efd9 87 /* Reading the coefficients */
simon 0:1014af42efd9 88 b0 = *pCoeffs++;
simon 0:1014af42efd9 89 b1 = *pCoeffs++;
simon 0:1014af42efd9 90 b2 = *pCoeffs++;
simon 0:1014af42efd9 91 a1 = *pCoeffs++;
simon 0:1014af42efd9 92 a2 = *pCoeffs++;
simon 0:1014af42efd9 93
simon 0:1014af42efd9 94 /* Reading the state values */
simon 0:1014af42efd9 95 Xn1 = pState[0];
simon 0:1014af42efd9 96 Xn2 = pState[1];
simon 0:1014af42efd9 97 Yn1 = pState[2];
simon 0:1014af42efd9 98 Yn2 = pState[3];
simon 0:1014af42efd9 99
simon 0:1014af42efd9 100 /* Apply loop unrolling and compute 4 output values simultaneously. */
simon 0:1014af42efd9 101 /* The variables acc ... acc3 hold output values that are being computed:
simon 0:1014af42efd9 102 *
simon 0:1014af42efd9 103 * acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]
simon 0:1014af42efd9 104 */
simon 0:1014af42efd9 105
simon 0:1014af42efd9 106 sample = blockSize >> 2u;
simon 0:1014af42efd9 107
simon 0:1014af42efd9 108 /* First part of the processing with loop unrolling. Compute 4 outputs at a time.
simon 0:1014af42efd9 109 ** a second loop below computes the remaining 1 to 3 samples. */
simon 0:1014af42efd9 110 while(sample > 0u)
simon 0:1014af42efd9 111 {
simon 0:1014af42efd9 112 /* Read the input */
simon 0:1014af42efd9 113 Xn = *pIn++;
simon 0:1014af42efd9 114
simon 0:1014af42efd9 115 /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
simon 0:1014af42efd9 116 /* acc = b0 * x[n] */
simon 0:1014af42efd9 117 acc = (q31_t) (((q63_t) b0 * Xn) >> 32);
simon 0:1014af42efd9 118 /* acc += b1 * x[n-1] */
simon 0:1014af42efd9 119 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn1))) >> 32);
simon 0:1014af42efd9 120 /* acc += b[2] * x[n-2] */
simon 0:1014af42efd9 121 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn2))) >> 32);
simon 0:1014af42efd9 122 /* acc += a1 * y[n-1] */
simon 0:1014af42efd9 123 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn1))) >> 32);
simon 0:1014af42efd9 124 /* acc += a2 * y[n-2] */
simon 0:1014af42efd9 125 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn2))) >> 32);
simon 0:1014af42efd9 126
simon 0:1014af42efd9 127 /* The result is converted to 1.31 , Yn2 variable is reused */
simon 0:1014af42efd9 128 Yn2 = acc << shift;
simon 0:1014af42efd9 129
simon 0:1014af42efd9 130 /* Store the output in the destination buffer. */
simon 0:1014af42efd9 131 *pOut++ = Yn2;
simon 0:1014af42efd9 132
simon 0:1014af42efd9 133 /* Read the second input */
simon 0:1014af42efd9 134 Xn2 = *pIn++;
simon 0:1014af42efd9 135
simon 0:1014af42efd9 136 /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
simon 0:1014af42efd9 137 /* acc = b0 * x[n] */
simon 0:1014af42efd9 138 acc = (q31_t) (((q63_t) b0 * (Xn2)) >> 32);
simon 0:1014af42efd9 139 /* acc += b1 * x[n-1] */
simon 0:1014af42efd9 140 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn))) >> 32);
simon 0:1014af42efd9 141 /* acc += b[2] * x[n-2] */
simon 0:1014af42efd9 142 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn1))) >> 32);
simon 0:1014af42efd9 143 /* acc += a1 * y[n-1] */
simon 0:1014af42efd9 144 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn2))) >> 32);
simon 0:1014af42efd9 145 /* acc += a2 * y[n-2] */
simon 0:1014af42efd9 146 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn1))) >> 32);
simon 0:1014af42efd9 147
simon 0:1014af42efd9 148 /* The result is converted to 1.31, Yn1 variable is reused */
simon 0:1014af42efd9 149 Yn1 = acc << shift;
simon 0:1014af42efd9 150
simon 0:1014af42efd9 151 /* Store the output in the destination buffer. */
simon 0:1014af42efd9 152 *pOut++ = Yn1;
simon 0:1014af42efd9 153
simon 0:1014af42efd9 154 /* Read the third input */
simon 0:1014af42efd9 155 Xn1 = *pIn++;
simon 0:1014af42efd9 156
simon 0:1014af42efd9 157 /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
simon 0:1014af42efd9 158 /* acc = b0 * x[n] */
simon 0:1014af42efd9 159 acc = (q31_t) (((q63_t) b0 * (Xn1)) >> 32);
simon 0:1014af42efd9 160 /* acc += b1 * x[n-1] */
simon 0:1014af42efd9 161 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn2))) >> 32);
simon 0:1014af42efd9 162 /* acc += b[2] * x[n-2] */
simon 0:1014af42efd9 163 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn))) >> 32);
simon 0:1014af42efd9 164 /* acc += a1 * y[n-1] */
simon 0:1014af42efd9 165 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn1))) >> 32);
simon 0:1014af42efd9 166 /* acc += a2 * y[n-2] */
simon 0:1014af42efd9 167 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn2))) >> 32);
simon 0:1014af42efd9 168
simon 0:1014af42efd9 169 /* The result is converted to 1.31, Yn2 variable is reused */
simon 0:1014af42efd9 170 Yn2 = acc << shift;
simon 0:1014af42efd9 171
simon 0:1014af42efd9 172 /* Store the output in the destination buffer. */
simon 0:1014af42efd9 173 *pOut++ = Yn2;
simon 0:1014af42efd9 174
simon 0:1014af42efd9 175 /* Read the forth input */
simon 0:1014af42efd9 176 Xn = *pIn++;
simon 0:1014af42efd9 177
simon 0:1014af42efd9 178 /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
simon 0:1014af42efd9 179 /* acc = b0 * x[n] */
simon 0:1014af42efd9 180 acc = (q31_t) (((q63_t) b0 * (Xn)) >> 32);
simon 0:1014af42efd9 181 /* acc += b1 * x[n-1] */
simon 0:1014af42efd9 182 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn1))) >> 32);
simon 0:1014af42efd9 183 /* acc += b[2] * x[n-2] */
simon 0:1014af42efd9 184 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn2))) >> 32);
simon 0:1014af42efd9 185 /* acc += a1 * y[n-1] */
simon 0:1014af42efd9 186 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn2))) >> 32);
simon 0:1014af42efd9 187 /* acc += a2 * y[n-2] */
simon 0:1014af42efd9 188 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn1))) >> 32);
simon 0:1014af42efd9 189
simon 0:1014af42efd9 190 /* The result is converted to 1.31, Yn1 variable is reused */
simon 0:1014af42efd9 191 Yn1 = acc << shift;
simon 0:1014af42efd9 192
simon 0:1014af42efd9 193 /* Every time after the output is computed state should be updated. */
simon 0:1014af42efd9 194 /* The states should be updated as: */
simon 0:1014af42efd9 195 /* Xn2 = Xn1 */
simon 0:1014af42efd9 196 /* Xn1 = Xn */
simon 0:1014af42efd9 197 /* Yn2 = Yn1 */
simon 0:1014af42efd9 198 /* Yn1 = acc */
simon 0:1014af42efd9 199 Xn2 = Xn1;
simon 0:1014af42efd9 200 Xn1 = Xn;
simon 0:1014af42efd9 201
simon 0:1014af42efd9 202 /* Store the output in the destination buffer. */
simon 0:1014af42efd9 203 *pOut++ = Yn1;
simon 0:1014af42efd9 204
simon 0:1014af42efd9 205 /* decrement the loop counter */
simon 0:1014af42efd9 206 sample--;
simon 0:1014af42efd9 207 }
simon 0:1014af42efd9 208
simon 0:1014af42efd9 209 /* If the blockSize is not a multiple of 4, compute any remaining output samples here.
simon 0:1014af42efd9 210 ** No loop unrolling is used. */
simon 0:1014af42efd9 211 sample = (blockSize & 0x3u);
simon 0:1014af42efd9 212
simon 0:1014af42efd9 213 while(sample > 0u)
simon 0:1014af42efd9 214 {
simon 0:1014af42efd9 215 /* Read the input */
simon 0:1014af42efd9 216 Xn = *pIn++;
simon 0:1014af42efd9 217
simon 0:1014af42efd9 218 /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
simon 0:1014af42efd9 219 /* acc = b0 * x[n] */
simon 0:1014af42efd9 220 acc = (q31_t) (((q63_t) b0 * (Xn)) >> 32);
simon 0:1014af42efd9 221 /* acc += b1 * x[n-1] */
simon 0:1014af42efd9 222 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn1))) >> 32);
simon 0:1014af42efd9 223 /* acc += b[2] * x[n-2] */
simon 0:1014af42efd9 224 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn2))) >> 32);
simon 0:1014af42efd9 225 /* acc += a1 * y[n-1] */
simon 0:1014af42efd9 226 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn1))) >> 32);
simon 0:1014af42efd9 227 /* acc += a2 * y[n-2] */
simon 0:1014af42efd9 228 acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn2))) >> 32);
simon 0:1014af42efd9 229 /* The result is converted to 1.31 */
simon 0:1014af42efd9 230 acc = acc << shift;
simon 0:1014af42efd9 231
simon 0:1014af42efd9 232 /* Every time after the output is computed state should be updated. */
simon 0:1014af42efd9 233 /* The states should be updated as: */
simon 0:1014af42efd9 234 /* Xn2 = Xn1 */
simon 0:1014af42efd9 235 /* Xn1 = Xn */
simon 0:1014af42efd9 236 /* Yn2 = Yn1 */
simon 0:1014af42efd9 237 /* Yn1 = acc */
simon 0:1014af42efd9 238 Xn2 = Xn1;
simon 0:1014af42efd9 239 Xn1 = Xn;
simon 0:1014af42efd9 240 Yn2 = Yn1;
simon 0:1014af42efd9 241 Yn1 = acc;
simon 0:1014af42efd9 242
simon 0:1014af42efd9 243 /* Store the output in the destination buffer. */
simon 0:1014af42efd9 244 *pOut++ = acc;
simon 0:1014af42efd9 245
simon 0:1014af42efd9 246 /* decrement the loop counter */
simon 0:1014af42efd9 247 sample--;
simon 0:1014af42efd9 248 }
simon 0:1014af42efd9 249
simon 0:1014af42efd9 250 /* The first stage goes from the input buffer to the output buffer. */
simon 0:1014af42efd9 251 /* Subsequent stages occur in-place in the output buffer */
simon 0:1014af42efd9 252 pIn = pDst;
simon 0:1014af42efd9 253
simon 0:1014af42efd9 254 /* Reset to destination pointer */
simon 0:1014af42efd9 255 pOut = pDst;
simon 0:1014af42efd9 256
simon 0:1014af42efd9 257 /* Store the updated state variables back into the pState array */
simon 0:1014af42efd9 258 *pState++ = Xn1;
simon 0:1014af42efd9 259 *pState++ = Xn2;
simon 0:1014af42efd9 260 *pState++ = Yn1;
simon 0:1014af42efd9 261 *pState++ = Yn2;
simon 0:1014af42efd9 262
simon 0:1014af42efd9 263 } while(--stage);
simon 0:1014af42efd9 264 }
simon 0:1014af42efd9 265
simon 0:1014af42efd9 266 /**
simon 0:1014af42efd9 267 * @} end of BiquadCascadeDF1 group
simon 0:1014af42efd9 268 */