lib/main/CMSIS/DSP/Source/StatisticsFunctions/arm_var_q31.c

   1 /* ----------------------------------------------------------------------
   2  * Project:      CMSIS DSP Library
   3  * Title:        arm_var_q31.c
   4  * Description:  Variance of an array of Q31 type
   5  *
   6  * $Date:        27. January 2017
   7  * $Revision:    V.1.5.1
   8  *
   9  * Target Processor: Cortex-M cores
  10  * -------------------------------------------------------------------- */
  11 /*
  12  * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
  13  *
  14  * SPDX-License-Identifier: Apache-2.0
  15  *
  16  * Licensed under the Apache License, Version 2.0 (the License); you may
  17  * not use this file except in compliance with the License.
  18  * You may obtain a copy of the License at
  19  *
  20  * www.apache.org/licenses/LICENSE-2.0
  21  *
  22  * Unless required by applicable law or agreed to in writing, software
  23  * distributed under the License is distributed on an AS IS BASIS, WITHOUT
  24  * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  25  * See the License for the specific language governing permissions and
  26  * limitations under the License.
  27  */
  28
  29 #include "arm_math.h"
  30
  31 /**
  32  * @ingroup groupStats
  33  */
  34
  35 /**
  36  * @addtogroup variance
  37  * @{
  38  */
  39
  40 /**
  41  * @brief Variance of the elements of a Q31 vector.
  42  * @param[in]       *pSrc points to the input vector
  43  * @param[in]       blockSize length of the input vector
  44  * @param[out]      *pResult variance value returned here
  45  * @return none.
  46  * @details
  47  * <b>Scaling and Overflow Behavior:</b>
  48  *
  49  *\par
  50  * The function is implemented using an internal 64-bit accumulator.
  51  * The input is represented in 1.31 format, which is then downshifted by 8 bits
  52  * which yields 1.23, and intermediate multiplication yields a 2.46 format.
  53  * The accumulator maintains full precision of the intermediate multiplication results,
  54  * but provides only a 16 guard bits.
  55  * There is no saturation on intermediate additions.
  56  * If the accumulator overflows it wraps around and distorts the result.
  57  * In order to avoid overflows completely the input signal must be scaled down by
  58  * log2(blockSize)-8 bits, as a total of blockSize additions are performed internally.
  59  * After division, internal variables should be Q18.46
  60  * Finally, the 18.46 accumulator is right shifted by 15 bits to yield a 1.31 format value.
  61  *
  62  */
  63
  64 void arm_var_q31(
  65   q31_t * pSrc,
  66   uint32_t blockSize,
  67   q31_t * pResult)
  68 {
  69   q63_t sum = 0;                                 /* Accumulator */
  70   q63_t meanOfSquares, squareOfMean;             /* square of mean and mean of square */
  71   q31_t in;                                      /* input value */
  72   uint32_t blkCnt;                               /* loop counter */
  73   q63_t sumOfSquares = 0;                        /* Accumulator */
  74
  75   if (blockSize == 1U)
  76   {
  77     *pResult = 0;
  78     return;
  79   }
  80
  81 #if defined (ARM_MATH_DSP)
  82   /* Run the below code for Cortex-M4 and Cortex-M3 */
  83
  84   /*loop Unrolling */
  85   blkCnt = blockSize >> 2U;
  86
  87   /* First part of the processing with loop unrolling.  Compute 4 outputs at a time.
  88    ** a second loop below computes the remaining 1 to 3 samples. */
  89   while (blkCnt > 0U)
  90   {
  91     /* C = (A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1])  */
  92     /* Compute Sum of squares of the input samples
  93      * and then store the result in a temporary variable, sum. */
  94     in = *pSrc++ >> 8U;
  95     sum += in;
  96     sumOfSquares += ((q63_t) (in) * (in));
  97     in = *pSrc++ >> 8U;
  98     sum += in;
  99     sumOfSquares += ((q63_t) (in) * (in));
 100     in = *pSrc++ >> 8U;
 101     sum += in;
 102     sumOfSquares += ((q63_t) (in) * (in));
 103     in = *pSrc++ >> 8U;
 104     sum += in;
 105     sumOfSquares += ((q63_t) (in) * (in));
 106
 107     /* Decrement the loop counter */
 108     blkCnt--;
 109   }
 110
 111   /* If the blockSize is not a multiple of 4, compute any remaining output samples here.
 112    ** No loop unrolling is used. */
 113   blkCnt = blockSize % 0x4U;
 114
 115   while (blkCnt > 0U)
 116   {
 117     /* C = (A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1]) */
 118     /* Compute Sum of squares of the input samples
 119      * and then store the result in a temporary variable, sum. */
 120     in = *pSrc++ >> 8U;
 121     sum += in;
 122     sumOfSquares += ((q63_t) (in) * (in));
 123
 124     /* Decrement the loop counter */
 125     blkCnt--;
 126   }
 127
 128   /* Compute Mean of squares of the input samples
 129    * and then store the result in a temporary variable, meanOfSquares. */
 130   meanOfSquares = sumOfSquares / (q63_t)(blockSize - 1U);
 131
 132 #else
 133   /* Run the below code for Cortex-M0 */
 134
 135   /* Loop over blockSize number of values */
 136   blkCnt = blockSize;
 137
 138   while (blkCnt > 0U)
 139   {
 140     /* C = (A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1]) */
 141     /* Compute Sum of squares of the input samples
 142      * and then store the result in a temporary variable, sumOfSquares. */
 143     in = *pSrc++ >> 8U;
 144     sumOfSquares += ((q63_t) (in) * (in));
 145
 146     /* C = (A[0] + A[1] + A[2] + ... + A[blockSize-1]) */
 147     /* Compute sum of all input values and then store the result in a temporary variable, sum. */
 148     sum += in;
 149
 150     /* Decrement the loop counter */
 151     blkCnt--;
 152   }
 153
 154   /* Compute Mean of squares of the input samples
 155    * and then store the result in a temporary variable, meanOfSquares. */
 156   meanOfSquares = sumOfSquares / (q63_t)(blockSize - 1U);
 157
 158 #endif /* #if defined (ARM_MATH_DSP) */
 159
 160   /* Compute square of mean */
 161   squareOfMean = sum * sum / (q63_t)(blockSize * (blockSize - 1U));
 162
 163   /* Compute standard deviation and then store the result to the destination */
 164   *pResult = (meanOfSquares - squareOfMean) >> 15U;
 165 }
 166
 167 /**
 168  * @} end of variance group
 169  */