1
2 /* ----------------------------------------------------------------------
3 * Project: CMSIS DSP Library
4 * Title: arm_canberra_distance_f32.c
5 * Description: Canberra distance between two vectors
6 *
7 * $Date: 23 April 2021
8 * $Revision: V1.9.0
9 *
10 * Target Processor: Cortex-M and Cortex-A cores
11 * -------------------------------------------------------------------- */
12 /*
13 * Copyright (C) 2010-2021 ARM Limited or its affiliates. All rights reserved.
14 *
15 * SPDX-License-Identifier: Apache-2.0
16 *
17 * Licensed under the Apache License, Version 2.0 (the License); you may
18 * not use this file except in compliance with the License.
19 * You may obtain a copy of the License at
20 *
21 * www.apache.org/licenses/LICENSE-2.0
22 *
23 * Unless required by applicable law or agreed to in writing, software
24 * distributed under the License is distributed on an AS IS BASIS, WITHOUT
25 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
26 * See the License for the specific language governing permissions and
27 * limitations under the License.
28 */
29
30 #include "dsp/distance_functions.h"
31 #include <limits.h>
32 #include <math.h>
33
34
35 /**
36 @addtogroup Canberra
37 @{
38 */
39
40
41 /**
42 * @brief Canberra distance between two vectors
43 *
44 * This function may divide by zero when samples pA[i] and pB[i] are both zero.
45 * The result of the computation will be correct. So the division per zero may be
46 * ignored.
47 *
48 * @param[in] pA First vector
49 * @param[in] pB Second vector
50 * @param[in] blockSize vector length
51 * @return distance
52 *
53 */
54
55 #if defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE)
56
57 #include "arm_helium_utils.h"
58 #include "arm_vec_math.h"
59
arm_canberra_distance_f32(const float32_t * pA,const float32_t * pB,uint32_t blockSize)60 float32_t arm_canberra_distance_f32(const float32_t *pA,const float32_t *pB, uint32_t blockSize)
61 {
62 float32_t accum = 0.0f;
63 uint32_t blkCnt;
64 f32x4_t a, b, c, accumV;
65
66 accumV = vdupq_n_f32(0.0f);
67
68 blkCnt = blockSize >> 2;
69 while (blkCnt > 0) {
70 a = vld1q(pA);
71 b = vld1q(pB);
72
73 c = vabdq(a, b);
74
75 a = vabsq(a);
76 b = vabsq(b);
77 a = vaddq(a, b);
78
79 /*
80 * May divide by zero when a and b have both the same lane at zero.
81 */
82 a = vrecip_medprec_f32(a);
83
84 /*
85 * Force result of a division by 0 to 0. It the behavior of the
86 * sklearn canberra function.
87 */
88 a = vdupq_m_n_f32(a, 0.0f, vcmpeqq(a, 0.0f));
89 c = vmulq(c, a);
90 accumV = vaddq(accumV, c);
91
92 pA += 4;
93 pB += 4;
94 blkCnt--;
95 }
96
97 blkCnt = blockSize & 3;
98 if (blkCnt > 0U) {
99 mve_pred16_t p0 = vctp32q(blkCnt);
100
101 a = vldrwq_z_f32(pA, p0);
102 b = vldrwq_z_f32(pB, p0);
103
104 c = vabdq(a, b);
105
106 a = vabsq(a);
107 b = vabsq(b);
108 a = vaddq(a, b);
109
110 /*
111 * May divide by zero when a and b have both the same lane at zero.
112 */
113 a = vrecip_medprec_f32(a);
114
115 /*
116 * Force result of a division by 0 to 0. It the behavior of the
117 * sklearn canberra function.
118 */
119 a = vdupq_m_n_f32(a, 0.0f, vcmpeqq(a, 0.0f));
120 c = vmulq(c, a);
121 accumV = vaddq_m(accumV, accumV, c, p0);
122 }
123
124 accum = vecAddAcrossF32Mve(accumV);
125
126 return (accum);
127 }
128
129 #else
130 #if defined(ARM_MATH_NEON)
131
132 #include "NEMath.h"
133
arm_canberra_distance_f32(const float32_t * pA,const float32_t * pB,uint32_t blockSize)134 float32_t arm_canberra_distance_f32(const float32_t *pA,const float32_t *pB, uint32_t blockSize)
135 {
136 float32_t accum=0.0f, tmpA, tmpB,diff,sum;
137 uint32_t blkCnt;
138 float32x4_t a,b,c,accumV;
139 float32x2_t accumV2;
140 uint32x4_t isZeroV;
141 float32x4_t zeroV = vdupq_n_f32(0.0f);
142
143 accumV = vdupq_n_f32(0.0f);
144
145 blkCnt = blockSize >> 2;
146 while(blkCnt > 0)
147 {
148 a = vld1q_f32(pA);
149 b = vld1q_f32(pB);
150
151 c = vabdq_f32(a,b);
152
153 a = vabsq_f32(a);
154 b = vabsq_f32(b);
155 a = vaddq_f32(a,b);
156 isZeroV = vceqq_f32(a,zeroV);
157
158 /*
159 * May divide by zero when a and b have both the same lane at zero.
160 */
161 a = vinvq_f32(a);
162
163 /*
164 * Force result of a division by 0 to 0. It the behavior of the
165 * sklearn canberra function.
166 */
167 a = vreinterpretq_f32_s32(vbicq_s32(vreinterpretq_s32_f32(a),vreinterpretq_s32_u32(isZeroV)));
168 c = vmulq_f32(c,a);
169 accumV = vaddq_f32(accumV,c);
170
171 pA += 4;
172 pB += 4;
173 blkCnt --;
174 }
175 accumV2 = vpadd_f32(vget_low_f32(accumV),vget_high_f32(accumV));
176 accum = vget_lane_f32(accumV2, 0) + vget_lane_f32(accumV2, 1);
177
178
179 blkCnt = blockSize & 3;
180 while(blkCnt > 0)
181 {
182 tmpA = *pA++;
183 tmpB = *pB++;
184
185 diff = fabsf(tmpA - tmpB);
186 sum = fabsf(tmpA) + fabsf(tmpB);
187 if ((tmpA != 0.0f) || (tmpB != 0.0f))
188 {
189 accum += (diff / sum);
190 }
191 blkCnt --;
192 }
193 return(accum);
194 }
195
196 #else
arm_canberra_distance_f32(const float32_t * pA,const float32_t * pB,uint32_t blockSize)197 float32_t arm_canberra_distance_f32(const float32_t *pA,const float32_t *pB, uint32_t blockSize)
198 {
199 float32_t accum=0.0f, tmpA, tmpB,diff,sum;
200
201 while(blockSize > 0)
202 {
203 tmpA = *pA++;
204 tmpB = *pB++;
205
206 diff = fabsf(tmpA - tmpB);
207 sum = fabsf(tmpA) + fabsf(tmpB);
208 if ((tmpA != 0.0f) || (tmpB != 0.0f))
209 {
210 accum += (diff / sum);
211 }
212 blockSize --;
213 }
214 return(accum);
215 }
216 #endif
217 #endif /* defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE) */
218
219
220 /**
221 * @} end of Canberra group
222 */
223
