1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_var_q15.c
4 * Description: Variance of an array of Q15 type
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 variance
37 @{
38 */
39
40 /**
41 @brief Variance 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 variance value returned here
45
46 @par Scaling and Overflow Behavior
47 The function is implemented using a 64-bit internal accumulator.
48 The input is represented in 1.15 format.
49 Intermediate multiplication yields a 2.30 format, and this
50 result is added without saturation to a 64-bit accumulator in 34.30 format.
51 With 33 guard bits in the accumulator, there is no risk of overflow, and the
52 full precision of the intermediate multiplication is preserved.
53 Finally, the 34.30 result is truncated to 34.15 format by discarding the lower
54 15 bits, and then saturated to yield a result in 1.15 format.
55 */
56 #if defined(ARM_MATH_MVEI) && !defined(ARM_MATH_AUTOVECTORIZE)
arm_var_q15(const q15_t * pSrc,uint32_t blockSize,q15_t * pResult)57 ARM_DSP_ATTRIBUTE void arm_var_q15(
58 const q15_t * pSrc,
59 uint32_t blockSize,
60 q15_t * pResult)
61 {
62 uint32_t blkCnt; /* loop counters */
63 q15x8_t vecSrc;
64 q63_t sumOfSquares = 0LL;
65 q63_t meanOfSquares, squareOfMean; /* square of mean and mean of square */
66 q63_t sum = 0LL;
67 q15_t in;
68
69 if (blockSize <= 1U) {
70 *pResult = 0;
71 return;
72 }
73
74
75 blkCnt = blockSize >> 3;
76 while (blkCnt > 0U)
77 {
78 vecSrc = vldrhq_s16(pSrc);
79 /* C = (A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1]) */
80 /* Compute Sum of squares of the input samples
81 * and then store the result in a temporary variable, sumOfSquares. */
82
83 sumOfSquares = vmlaldavaq_s16(sumOfSquares, vecSrc, vecSrc);
84 sum = vaddvaq_s16(sum, vecSrc);
85
86 blkCnt --;
87 pSrc += 8;
88 }
89
90 /* Tail */
91 blkCnt = blockSize & 7;
92 while (blkCnt > 0U)
93 {
94 /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
95 /* C = A[0] + A[1] + ... + A[blockSize-1] */
96
97 in = *pSrc++;
98 /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
99 #if defined (ARM_MATH_DSP)
100 sumOfSquares = __SMLALD(in, in, sumOfSquares);
101 #else
102 sumOfSquares += (in * in);
103 #endif /* #if defined (ARM_MATH_DSP) */
104 /* Compute sum and store result in a temporary variable, sum. */
105 sum += in;
106
107 /* Decrement loop counter */
108 blkCnt--;
109 }
110
111 /* Compute Mean of squares of the input samples
112 * and then store the result in a temporary variable, meanOfSquares. */
113 meanOfSquares = arm_div_int64_to_int32(sumOfSquares, (blockSize - 1U));
114
115 /* Compute square of mean */
116 squareOfMean = arm_div_int64_to_int32((q63_t)sum * sum, (q31_t)(blockSize * (blockSize - 1U)));
117
118 /* mean of the squares minus the square of the mean. */
119 *pResult = (meanOfSquares - squareOfMean) >> 15;
120 }
121 #else
arm_var_q15(const q15_t * pSrc,uint32_t blockSize,q15_t * pResult)122 ARM_DSP_ATTRIBUTE void arm_var_q15(
123 const q15_t * pSrc,
124 uint32_t blockSize,
125 q15_t * pResult)
126 {
127 uint32_t blkCnt; /* Loop counter */
128 q31_t sum = 0; /* Accumulator */
129 q31_t meanOfSquares, squareOfMean; /* Square of mean and mean of square */
130 q63_t sumOfSquares = 0; /* Sum of squares */
131 q15_t in; /* Temporary variable to store input value */
132
133 #if defined (ARM_MATH_LOOPUNROLL) && defined (ARM_MATH_DSP)
134 q31_t in32; /* Temporary variable to store input value */
135 #endif
136
137 if (blockSize <= 1U)
138 {
139 *pResult = 0;
140 return;
141 }
142
143 #if defined (ARM_MATH_LOOPUNROLL)
144
145 /* Loop unrolling: Compute 4 outputs at a time */
146 blkCnt = blockSize >> 2U;
147
148 while (blkCnt > 0U)
149 {
150 /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
151 /* C = A[0] + A[1] + ... + A[blockSize-1] */
152
153 /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
154 /* Compute sum and store result in a temporary variable, sum. */
155 #if defined (ARM_MATH_DSP)
156 in32 = read_q15x2_ia (&pSrc);
157 sumOfSquares = __SMLALD(in32, in32, sumOfSquares);
158 sum += ((in32 << 16U) >> 16U);
159 sum += (in32 >> 16U);
160
161 in32 = read_q15x2_ia (&pSrc);
162 sumOfSquares = __SMLALD(in32, in32, sumOfSquares);
163 sum += ((in32 << 16U) >> 16U);
164 sum += (in32 >> 16U);
165 #else
166 in = *pSrc++;
167 sumOfSquares += (in * in);
168 sum += in;
169
170 in = *pSrc++;
171 sumOfSquares += (in * in);
172 sum += in;
173
174 in = *pSrc++;
175 sumOfSquares += (in * in);
176 sum += in;
177
178 in = *pSrc++;
179 sumOfSquares += (in * in);
180 sum += in;
181 #endif /* #if defined (ARM_MATH_DSP) */
182
183 /* Decrement loop counter */
184 blkCnt--;
185 }
186
187 /* Loop unrolling: Compute remaining outputs */
188 blkCnt = blockSize % 0x4U;
189
190 #else
191
192 /* Initialize blkCnt with number of samples */
193 blkCnt = blockSize;
194
195 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
196
197 while (blkCnt > 0U)
198 {
199 /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
200 /* C = A[0] + A[1] + ... + A[blockSize-1] */
201
202 in = *pSrc++;
203 /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
204 #if defined (ARM_MATH_DSP)
205 sumOfSquares = __SMLALD(in, in, sumOfSquares);
206 #else
207 sumOfSquares += (in * in);
208 #endif /* #if defined (ARM_MATH_DSP) */
209 /* Compute sum and store result in a temporary variable, sum. */
210 sum += in;
211
212 /* Decrement loop counter */
213 blkCnt--;
214 }
215
216 /* Compute Mean of squares and store result in a temporary variable, meanOfSquares. */
217 meanOfSquares = (q31_t) (sumOfSquares / (q63_t)(blockSize - 1U));
218
219 /* Compute square of mean */
220 squareOfMean = (q31_t) ((q63_t) sum * sum / (q63_t)(blockSize * (blockSize - 1U)));
221
222 /* mean of squares minus the square of mean. */
223 *pResult = (meanOfSquares - squareOfMean) >> 15U;
224 }
225 #endif /* defined(ARM_MATH_MVEI) */
226
227 /**
228 @} end of variance group
229 */
230
