1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_var_q31.c
4 * Description: Variance of an array of Q31 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 Q31 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 an internal 64-bit accumulator.
48 The input is represented in 1.31 format, which is then downshifted by 8 bits
49 which yields 1.23, and intermediate multiplication yields a 2.46 format.
50 The accumulator maintains full precision of the intermediate multiplication results,
51 and as a consequence has only 16 guard bits.
52 There is no saturation on intermediate additions.
53 If the accumulator overflows it wraps around and distorts the result.
54 In order to avoid overflows completely the input signal must be scaled down by
55 log2(blockSize)-8 bits, as a total of blockSize additions are performed internally.
56 After division, internal variables should be Q18.46
57 Finally, the 18.46 accumulator is right shifted by 15 bits to yield a 1.31 format value.
58 */
59 #if defined(ARM_MATH_MVEI) && !defined(ARM_MATH_AUTOVECTORIZE)
arm_var_q31(const q31_t * pSrc,uint32_t blockSize,q31_t * pResult)60 ARM_DSP_ATTRIBUTE void arm_var_q31(
61 const q31_t * pSrc,
62 uint32_t blockSize,
63 q31_t * pResult)
64 {
65 uint32_t blkCnt; /* loop counters */
66 q31x4_t vecSrc;
67 q63_t sumOfSquares = 0LL;
68 q63_t meanOfSquares, squareOfMean; /* square of mean and mean of square */
69 q63_t sum = 0LL;
70 q31_t in;
71
72 if (blockSize <= 1U) {
73 *pResult = 0;
74 return;
75 }
76
77
78 /* Compute 4 outputs at a time */
79 blkCnt = blockSize >> 2U;
80 while (blkCnt > 0U)
81 {
82 vecSrc = vldrwq_s32(pSrc);
83 /* C = (A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1]) */
84 /* Compute Sum of squares of the input samples
85 * and then store the result in a temporary variable, sumOfSquares. */
86
87 /* downscale */
88 vecSrc = vshrq(vecSrc, 8);
89 sumOfSquares = vmlaldavaq(sumOfSquares, vecSrc, vecSrc);
90 sum = vaddlvaq(sum, vecSrc);
91
92 blkCnt --;
93 pSrc += 4;
94 }
95
96
97 /*
98 * tail
99 */
100 blkCnt = blockSize & 0x3;
101 while (blkCnt > 0U)
102 {
103 /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
104 /* C = A[0] + A[1] + ... + A[blockSize-1] */
105
106 in = *pSrc++ >> 8U;
107 /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
108 sumOfSquares += ((q63_t) (in) * (in));
109 /* Compute sum and store result in a temporary variable, sum. */
110 sum += in;
111
112 /* Decrement loop counter */
113 blkCnt--;
114 }
115
116 /* Compute Mean of squares of the input samples
117 * and then store the result in a temporary variable, meanOfSquares. */
118 meanOfSquares = sumOfSquares / (q63_t) (blockSize - 1U);
119
120 /* Compute square of mean */
121 squareOfMean = sum * sum / (q63_t) (blockSize * (blockSize - 1U));
122
123 /* Compute standard deviation and then store the result to the destination */
124 *pResult = asrl(meanOfSquares - squareOfMean, 15U);
125 }
126 #else
arm_var_q31(const q31_t * pSrc,uint32_t blockSize,q31_t * pResult)127 ARM_DSP_ATTRIBUTE void arm_var_q31(
128 const q31_t * pSrc,
129 uint32_t blockSize,
130 q31_t * pResult)
131 {
132 uint32_t blkCnt; /* Loop counter */
133 q63_t sum = 0; /* Temporary result storage */
134 q63_t meanOfSquares, squareOfMean; /* Square of mean and mean of square */
135 q63_t sumOfSquares = 0; /* Sum of squares */
136 q31_t in; /* Temporary variable to store input value */
137
138 if (blockSize <= 1U)
139 {
140 *pResult = 0;
141 return;
142 }
143
144 #if defined (ARM_MATH_LOOPUNROLL)
145
146 /* Loop unrolling: Compute 4 outputs at a time */
147 blkCnt = blockSize >> 2U;
148
149 while (blkCnt > 0U)
150 {
151 /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
152 /* C = A[0] + A[1] + ... + A[blockSize-1] */
153
154 in = *pSrc++ >> 8U;
155 /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
156 sumOfSquares += ((q63_t) (in) * (in));
157 /* Compute sum and store result in a temporary variable, sum. */
158 sum += in;
159
160 in = *pSrc++ >> 8U;
161 sumOfSquares += ((q63_t) (in) * (in));
162 sum += in;
163
164 in = *pSrc++ >> 8U;
165 sumOfSquares += ((q63_t) (in) * (in));
166 sum += in;
167
168 in = *pSrc++ >> 8U;
169 sumOfSquares += ((q63_t) (in) * (in));
170 sum += in;
171
172 /* Decrement loop counter */
173 blkCnt--;
174 }
175
176 /* Loop unrolling: Compute remaining outputs */
177 blkCnt = blockSize % 0x4U;
178
179 #else
180
181 /* Initialize blkCnt with number of samples */
182 blkCnt = blockSize;
183
184 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
185
186 while (blkCnt > 0U)
187 {
188 /* C = A[0] * A[0] + A[1] * A[1] + ... + A[blockSize-1] * A[blockSize-1] */
189 /* C = A[0] + A[1] + ... + A[blockSize-1] */
190
191 in = *pSrc++ >> 8U;
192 /* Compute sum of squares and store result in a temporary variable, sumOfSquares. */
193 sumOfSquares += ((q63_t) (in) * (in));
194 /* Compute sum and store result in a temporary variable, sum. */
195 sum += in;
196
197 /* Decrement loop counter */
198 blkCnt--;
199 }
200
201 /* Compute Mean of squares and store result in a temporary variable, meanOfSquares. */
202 meanOfSquares = (sumOfSquares / (q63_t)(blockSize - 1U));
203
204 /* Compute square of mean */
205 squareOfMean = ( sum * sum / (q63_t)(blockSize * (blockSize - 1U)));
206
207 /* Compute variance and store result in destination */
208 *pResult = (meanOfSquares - squareOfMean) >> 15U;
209 }
210 #endif
211 /**
212 @} end of variance group
213 */
214
