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_lms_q31.c
simon 0:1014af42efd9 9 *
simon 0:1014af42efd9 10 * Description: Processing function for the Q31 LMS filter.
simon 0:1014af42efd9 11 *
simon 0:1014af42efd9 12 * Target Processor: Cortex-M4/Cortex-M3
simon 0:1014af42efd9 13 *
simon 0:1014af42efd9 14 * Version 1.0.3 2010/11/29
simon 0:1014af42efd9 15 * Re-organized the CMSIS folders and updated documentation.
simon 0:1014af42efd9 16 *
simon 0:1014af42efd9 17 * Version 1.0.2 2010/11/11
simon 0:1014af42efd9 18 * Documentation updated.
simon 0:1014af42efd9 19 *
simon 0:1014af42efd9 20 * Version 1.0.1 2010/10/05
simon 0:1014af42efd9 21 * Production release and review comments incorporated.
simon 0:1014af42efd9 22 *
simon 0:1014af42efd9 23 * Version 1.0.0 2010/09/20
simon 0:1014af42efd9 24 * Production release and review comments incorporated
simon 0:1014af42efd9 25 *
simon 0:1014af42efd9 26 * Version 0.0.7 2010/06/10
simon 0:1014af42efd9 27 * Misra-C changes done
simon 0:1014af42efd9 28 * -------------------------------------------------------------------- */
simon 0:1014af42efd9 29
simon 0:1014af42efd9 30 #include "arm_math.h"
simon 0:1014af42efd9 31 /**
simon 0:1014af42efd9 32 * @ingroup groupFilters
simon 0:1014af42efd9 33 */
simon 0:1014af42efd9 34
simon 0:1014af42efd9 35 /**
simon 0:1014af42efd9 36 * @addtogroup LMS
simon 0:1014af42efd9 37 * @{
simon 0:1014af42efd9 38 */
simon 0:1014af42efd9 39
simon 0:1014af42efd9 40 /**
simon 0:1014af42efd9 41 * @brief Processing function for Q31 LMS filter.
simon 0:1014af42efd9 42 * @param[in] *S points to an instance of the Q15 LMS filter structure.
simon 0:1014af42efd9 43 * @param[in] *pSrc points to the block of input data.
simon 0:1014af42efd9 44 * @param[in] *pRef points to the block of reference data.
simon 0:1014af42efd9 45 * @param[out] *pOut points to the block of output data.
simon 0:1014af42efd9 46 * @param[out] *pErr points to the block of error data.
simon 0:1014af42efd9 47 * @param[in] blockSize number of samples to process.
simon 0:1014af42efd9 48 * @return none.
simon 0:1014af42efd9 49 *
simon 0:1014af42efd9 50 * \par Scaling and Overflow Behavior:
simon 0:1014af42efd9 51 * The function is implemented using an internal 64-bit accumulator.
simon 0:1014af42efd9 52 * The accumulator has a 2.62 format and maintains full precision of the intermediate multiplication results but provides only a single guard bit.
simon 0:1014af42efd9 53 * Thus, if the accumulator result overflows it wraps around rather than clips.
simon 0:1014af42efd9 54 * In order to avoid overflows completely the input signal must be scaled down by log2(numTaps) bits.
simon 0:1014af42efd9 55 * The reference signal should not be scaled down.
simon 0:1014af42efd9 56 * After all multiply-accumulates are performed, the 2.62 accumulator is truncated to 1.32 format and then saturated to 1.31 format.
simon 0:1014af42efd9 57 * The output signal and error signal are in 1.31 format.
simon 0:1014af42efd9 58 *
simon 0:1014af42efd9 59 * \par
simon 0:1014af42efd9 60 * In this filter, filter coefficients are updated for each sample and the updation of filter cofficients are saturted.
simon 0:1014af42efd9 61 */
simon 0:1014af42efd9 62
simon 0:1014af42efd9 63 void arm_lms_q31(
simon 0:1014af42efd9 64 const arm_lms_instance_q31 * S,
simon 0:1014af42efd9 65 q31_t * pSrc,
simon 0:1014af42efd9 66 q31_t * pRef,
simon 0:1014af42efd9 67 q31_t * pOut,
simon 0:1014af42efd9 68 q31_t * pErr,
simon 0:1014af42efd9 69 uint32_t blockSize)
simon 0:1014af42efd9 70 {
simon 0:1014af42efd9 71 q31_t *pState = S->pState; /* State pointer */
simon 0:1014af42efd9 72 uint32_t numTaps = S->numTaps; /* Number of filter coefficients in the filter */
simon 0:1014af42efd9 73 q31_t *pCoeffs = S->pCoeffs; /* Coefficient pointer */
simon 0:1014af42efd9 74 q31_t *pStateCurnt; /* Points to the current sample of the state */
simon 0:1014af42efd9 75 q31_t mu = S->mu; /* Adaptive factor */
simon 0:1014af42efd9 76 q31_t *px; /* Temporary pointer for state */
simon 0:1014af42efd9 77 q31_t *pb; /* Temporary pointer for coefficient buffer */
simon 0:1014af42efd9 78 uint32_t tapCnt, blkCnt; /* Loop counters */
simon 0:1014af42efd9 79 q63_t acc; /* Accumulator */
simon 0:1014af42efd9 80 q31_t e = 0; /* error of data sample */
simon 0:1014af42efd9 81 q31_t alpha; /* Intermediate constant for taps update */
simon 0:1014af42efd9 82 uint8_t shift = (uint8_t) (32u - (S->postShift + 1u)); /* Shift to be applied to the output */
simon 0:1014af42efd9 83 q31_t coef; /* Temporary variable for coef */
simon 0:1014af42efd9 84
simon 0:1014af42efd9 85 /* S->pState points to buffer which contains previous frame (numTaps - 1) samples */
simon 0:1014af42efd9 86 /* pStateCurnt points to the location where the new input data should be written */
simon 0:1014af42efd9 87 pStateCurnt = &(S->pState[(numTaps - 1u)]);
simon 0:1014af42efd9 88
simon 0:1014af42efd9 89 /* Initializing blkCnt with blockSize */
simon 0:1014af42efd9 90 blkCnt = blockSize;
simon 0:1014af42efd9 91
simon 0:1014af42efd9 92 while(blkCnt > 0u)
simon 0:1014af42efd9 93 {
simon 0:1014af42efd9 94 /* Copy the new input sample into the state buffer */
simon 0:1014af42efd9 95 *pStateCurnt++ = *pSrc++;
simon 0:1014af42efd9 96
simon 0:1014af42efd9 97 /* Initialize state pointer */
simon 0:1014af42efd9 98 px = pState;
simon 0:1014af42efd9 99
simon 0:1014af42efd9 100 /* Initialize coefficient pointer */
simon 0:1014af42efd9 101 pb = pCoeffs;
simon 0:1014af42efd9 102
simon 0:1014af42efd9 103 /* Set the accumulator to zero */
simon 0:1014af42efd9 104 acc = 0;
simon 0:1014af42efd9 105
simon 0:1014af42efd9 106 /* Loop unrolling. Process 4 taps at a time. */
simon 0:1014af42efd9 107 tapCnt = numTaps >> 2;
simon 0:1014af42efd9 108
simon 0:1014af42efd9 109 while(tapCnt > 0u)
simon 0:1014af42efd9 110 {
simon 0:1014af42efd9 111 /* Perform the multiply-accumulate */
simon 0:1014af42efd9 112 /* acc += b[N] * x[n-N] */
simon 0:1014af42efd9 113 acc += ((q63_t) (*px++)) * (*pb++);
simon 0:1014af42efd9 114
simon 0:1014af42efd9 115 /* acc += b[N-1] * x[n-N-1] */
simon 0:1014af42efd9 116 acc += ((q63_t) (*px++)) * (*pb++);
simon 0:1014af42efd9 117
simon 0:1014af42efd9 118 /* acc += b[N-2] * x[n-N-2] */
simon 0:1014af42efd9 119 acc += ((q63_t) (*px++)) * (*pb++);
simon 0:1014af42efd9 120
simon 0:1014af42efd9 121 /* acc += b[N-3] * x[n-N-3] */
simon 0:1014af42efd9 122 acc += ((q63_t) (*px++)) * (*pb++);
simon 0:1014af42efd9 123
simon 0:1014af42efd9 124 /* Decrement the loop counter */
simon 0:1014af42efd9 125 tapCnt--;
simon 0:1014af42efd9 126 }
simon 0:1014af42efd9 127
simon 0:1014af42efd9 128 /* If the filter length is not a multiple of 4, compute the remaining filter taps */
simon 0:1014af42efd9 129 tapCnt = numTaps % 0x4u;
simon 0:1014af42efd9 130
simon 0:1014af42efd9 131 while(tapCnt > 0u)
simon 0:1014af42efd9 132 {
simon 0:1014af42efd9 133 /* Perform the multiply-accumulate */
simon 0:1014af42efd9 134 acc += ((q63_t) (*px++)) * (*pb++);
simon 0:1014af42efd9 135
simon 0:1014af42efd9 136 /* Decrement the loop counter */
simon 0:1014af42efd9 137 tapCnt--;
simon 0:1014af42efd9 138 }
simon 0:1014af42efd9 139
simon 0:1014af42efd9 140 /* Converting the result to 1.31 format */
simon 0:1014af42efd9 141 /* Store the result from accumulator into the destination buffer. */
simon 0:1014af42efd9 142 acc = (q31_t) (acc >> shift);
simon 0:1014af42efd9 143
simon 0:1014af42efd9 144 *pOut++ = (q31_t) acc;
simon 0:1014af42efd9 145
simon 0:1014af42efd9 146 /* Compute and store error */
simon 0:1014af42efd9 147 e = *pRef++ - (q31_t) acc;
simon 0:1014af42efd9 148
simon 0:1014af42efd9 149 *pErr++ = (q31_t) e;
simon 0:1014af42efd9 150
simon 0:1014af42efd9 151 /* Compute alpha i.e. intermediate constant for taps update */
simon 0:1014af42efd9 152 alpha = (q31_t) (((q63_t) e * mu) >> 31);
simon 0:1014af42efd9 153
simon 0:1014af42efd9 154 /* Initialize state pointer */
simon 0:1014af42efd9 155 /* Advance state pointer by 1 for the next sample */
simon 0:1014af42efd9 156 px = pState++;
simon 0:1014af42efd9 157
simon 0:1014af42efd9 158 /* Initialize coefficient pointer */
simon 0:1014af42efd9 159 pb = pCoeffs;
simon 0:1014af42efd9 160
simon 0:1014af42efd9 161 /* Loop unrolling. Process 4 taps at a time. */
simon 0:1014af42efd9 162 tapCnt = numTaps >> 2;
simon 0:1014af42efd9 163
simon 0:1014af42efd9 164 /* Update filter coefficients */
simon 0:1014af42efd9 165 while(tapCnt > 0u)
simon 0:1014af42efd9 166 {
simon 0:1014af42efd9 167 /* coef is in 2.30 format */
simon 0:1014af42efd9 168 coef = (q31_t) (((q63_t) alpha * (*px++)) >> (32));
simon 0:1014af42efd9 169 /* get coef in 1.31 format by left shifting */
simon 0:1014af42efd9 170 *pb = clip_q63_to_q31((q63_t) *pb + (coef << 1u));
simon 0:1014af42efd9 171 /* update coefficient buffer to next coefficient */
simon 0:1014af42efd9 172 pb++;
simon 0:1014af42efd9 173
simon 0:1014af42efd9 174 coef = (q31_t) (((q63_t) alpha * (*px++)) >> (32));
simon 0:1014af42efd9 175 *pb = clip_q63_to_q31((q63_t) *pb + (coef << 1u));
simon 0:1014af42efd9 176 pb++;
simon 0:1014af42efd9 177
simon 0:1014af42efd9 178 coef = (q31_t) (((q63_t) alpha * (*px++)) >> (32));
simon 0:1014af42efd9 179 *pb = clip_q63_to_q31((q63_t) *pb + (coef << 1u));
simon 0:1014af42efd9 180 pb++;
simon 0:1014af42efd9 181
simon 0:1014af42efd9 182 coef = (q31_t) (((q63_t) alpha * (*px++)) >> (32));
simon 0:1014af42efd9 183 *pb = clip_q63_to_q31((q63_t) *pb + (coef << 1u));
simon 0:1014af42efd9 184 pb++;
simon 0:1014af42efd9 185
simon 0:1014af42efd9 186 /* Decrement the loop counter */
simon 0:1014af42efd9 187 tapCnt--;
simon 0:1014af42efd9 188 }
simon 0:1014af42efd9 189
simon 0:1014af42efd9 190 /* If the filter length is not a multiple of 4, compute the remaining filter taps */
simon 0:1014af42efd9 191 tapCnt = numTaps % 0x4u;
simon 0:1014af42efd9 192
simon 0:1014af42efd9 193 while(tapCnt > 0u)
simon 0:1014af42efd9 194 {
simon 0:1014af42efd9 195 /* Perform the multiply-accumulate */
simon 0:1014af42efd9 196 coef = (q31_t) (((q63_t) alpha * (*px++)) >> (32));
simon 0:1014af42efd9 197 *pb = clip_q63_to_q31((q63_t) *pb + (coef << 1u));
simon 0:1014af42efd9 198 pb++;
simon 0:1014af42efd9 199
simon 0:1014af42efd9 200 /* Decrement the loop counter */
simon 0:1014af42efd9 201 tapCnt--;
simon 0:1014af42efd9 202 }
simon 0:1014af42efd9 203
simon 0:1014af42efd9 204 /* Decrement the loop counter */
simon 0:1014af42efd9 205 blkCnt--;
simon 0:1014af42efd9 206 }
simon 0:1014af42efd9 207
simon 0:1014af42efd9 208 /* Processing is complete. Now copy the last numTaps - 1 samples to the
simon 0:1014af42efd9 209 satrt of the state buffer. This prepares the state buffer for the
simon 0:1014af42efd9 210 next function call. */
simon 0:1014af42efd9 211
simon 0:1014af42efd9 212 /* Points to the start of the pState buffer */
simon 0:1014af42efd9 213 pStateCurnt = S->pState;
simon 0:1014af42efd9 214
simon 0:1014af42efd9 215 /* Loop unrolling for (numTaps - 1u) samples copy */
simon 0:1014af42efd9 216 tapCnt = (numTaps - 1u) >> 2u;
simon 0:1014af42efd9 217
simon 0:1014af42efd9 218 /* copy data */
simon 0:1014af42efd9 219 while(tapCnt > 0u)
simon 0:1014af42efd9 220 {
simon 0:1014af42efd9 221 *pStateCurnt++ = *pState++;
simon 0:1014af42efd9 222 *pStateCurnt++ = *pState++;
simon 0:1014af42efd9 223 *pStateCurnt++ = *pState++;
simon 0:1014af42efd9 224 *pStateCurnt++ = *pState++;
simon 0:1014af42efd9 225
simon 0:1014af42efd9 226 /* Decrement the loop counter */
simon 0:1014af42efd9 227 tapCnt--;
simon 0:1014af42efd9 228 }
simon 0:1014af42efd9 229
simon 0:1014af42efd9 230 /* Calculate remaining number of copies */
simon 0:1014af42efd9 231 tapCnt = (numTaps - 1u) % 0x4u;
simon 0:1014af42efd9 232
simon 0:1014af42efd9 233 /* Copy the remaining q31_t data */
simon 0:1014af42efd9 234 while(tapCnt > 0u)
simon 0:1014af42efd9 235 {
simon 0:1014af42efd9 236 *pStateCurnt++ = *pState++;
simon 0:1014af42efd9 237
simon 0:1014af42efd9 238 /* Decrement the loop counter */
simon 0:1014af42efd9 239 tapCnt--;
simon 0:1014af42efd9 240 }
simon 0:1014af42efd9 241
simon 0:1014af42efd9 242 }
simon 0:1014af42efd9 243
simon 0:1014af42efd9 244 /**
simon 0:1014af42efd9 245 * @} end of LMS group
simon 0:1014af42efd9 246 */