1  /* ----------------------------------------------------------------------
  2   * Project:      CMSIS DSP Library
  3   * Title:        arm_std_q15.c
  4   * Description:  Standard deviation of an array of Q15 vector
  5   *
  6   * $Date:        23 April 2021
  7   * $Revision:    V1.9.0
  8   *
  9   * Target Processor: Cortex-M and Cortex-A cores
 10   * -------------------------------------------------------------------- */
 11  /*
 12   * Copyright (C) 2010-2021 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 "dsp/statistics_functions.h"
 30  
 31  /**
 32    @ingroup groupStats
 33   */
 34  
 35  /**
 36    @addtogroup STD
 37    @{
 38   */
 39  
 40  /**
 41    @brief         Standard deviation of the elements of a Q15 vector.
 42    @param[in]     pSrc       points to the input vector
 43    @param[in]     blockSize  number of samples in input vector
 44    @param[out]    pResult    standard deviation value returned here
 45    @return        none
 46  
 47    @par           Scaling and Overflow Behavior
 48                     The function is implemented using a 64-bit internal accumulator.
 49                     The input is represented in 1.15 format.
 50                     Intermediate multiplication yields a 2.30 format, and this
 51                     result is added without saturation to a 64-bit accumulator in 34.30 format.
 52                     With 33 guard bits in the accumulator, there is no risk of overflow, and the
 53                     full precision of the intermediate multiplication is preserved.
 54                     Finally, the 34.30 result is truncated to 34.15 format by discarding the lower
 55                     15 bits, and then saturated to yield a result in 1.15 format.
 56   */
 57  #if defined(ARM_MATH_MVEI) && !defined(ARM_MATH_AUTOVECTORIZE)
 58  void arm_std_q15(
 59    const q15_t * pSrc,
 60          uint32_t blockSize,
 61          q15_t * pResult)
 62  {
 63      q15_t var=0;
 64  
 65      arm_var_q15(pSrc, blockSize, &var);
 66      arm_sqrt_q15(var,pResult);
 67  }
 68  #else
 69  void arm_std_q15(
 70    const q15_t * pSrc,
 71          uint32_t blockSize,
 72          q15_t * pResult)
 73  {
 74          uint32_t blkCnt;                               /* Loop counter */
 75          q31_t sum = 0;                                 /* Accumulator */
 76          q31_t meanOfSquares, squareOfMean;             /* Square of mean and mean of square */
 77          q63_t sumOfSquares = 0;                        /* Sum of squares */
 78          q15_t in;                                      /* Temporary variable to store input value */
 79  
 80  #if defined (ARM_MATH_LOOPUNROLL) && defined (ARM_MATH_DSP)
 81    q31_t in32;                                    /* Temporary variable to store input value */
 82  #endif
 83  
 84    if (blockSize <= 1U)
 85    {
 86      *pResult = 0;
 87      return;
 88    }
 89  
 90  #if defined (ARM_MATH_LOOPUNROLL)
 91  
 92    /* Loop unrolling: Compute 4 outputs at a time */
 93    blkCnt = blockSize >> 2U;
 94  
 95    while (blkCnt > 0U)
 96    {
 97      /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
 98      /* C = A[0] + A[1] + ... + A[blockSize-1] */
 99  
100      /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
101      /* Compute sum and store result in a temporary variable, sum. */
102  #if defined (ARM_MATH_DSP)
103      in32 = read_q15x2_ia (&pSrc);
104      sumOfSquares = __SMLALD(in32, in32, sumOfSquares);
105      sum += ((in32 << 16U) >> 16U);
106      sum +=  (in32 >> 16U);
107  
108      in32 = read_q15x2_ia (&pSrc);
109      sumOfSquares = __SMLALD(in32, in32, sumOfSquares);
110      sum += ((in32 << 16U) >> 16U);
111      sum +=  (in32 >> 16U);
112  #else
113      in = *pSrc++;
114      sumOfSquares += (in * in);
115      sum += in;
116  
117      in = *pSrc++;
118      sumOfSquares += (in * in);
119      sum += in;
120  
121      in = *pSrc++;
122      sumOfSquares += (in * in);
123      sum += in;
124  
125      in = *pSrc++;
126      sumOfSquares += (in * in);
127      sum += in;
128  #endif /* #if defined (ARM_MATH_DSP) */
129  
130      /* Decrement loop counter */
131      blkCnt--;
132    }
133  
134    /* Loop unrolling: Compute remaining outputs */
135    blkCnt = blockSize % 0x4U;
136  
137  #else
138  
139    /* Initialize blkCnt with number of samples */
140    blkCnt = blockSize;
141  
142  #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
143  
144    while (blkCnt > 0U)
145    {
146      /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
147      /* C = A[0] + A[1] + ... + A[blockSize-1] */
148  
149      in = *pSrc++;
150      /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
151      sumOfSquares += (in * in);
152      /* Compute sum and store result in a temporary variable, sum. */
153      sum += in;
154  
155      /* Decrement loop counter */
156      blkCnt--;
157    }
158  
159    /* Compute Mean of squares and store result in a temporary variable, meanOfSquares. */
160    meanOfSquares = (q31_t) (sumOfSquares / (q63_t)(blockSize - 1U));
161  
162    /* Compute square of mean */
163    squareOfMean = (q31_t) ((q63_t) sum * sum / (q63_t)(blockSize * (blockSize - 1U)));
164  
165    /* mean of squares minus the square of mean. */
166    /* Compute standard deviation and store result in destination */
167    arm_sqrt_q15(__SSAT((meanOfSquares - squareOfMean) >> 15U, 16U), pResult);
168  }
169  #endif /* defined(ARM_MATH_MVEI) */
170  
171  /**
172    @} end of STD group
173   */