1
2 /* ----------------------------------------------------------------------
3 * Project: CMSIS DSP Library
4 * Title: arm_braycurtis_distance_f32.c
5 * Description: Bray-Curtis 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 "dsp/matrix_utils.h"
32 #include <limits.h>
33 #include <math.h>
34
35 #define E(MAT,R,C) \
36 (*((MAT)->pData + (MAT)->numCols*(R) + (C)))
37
38 #define WIN(R,C) \
39 ((pWindow == NULL) ? 1 : \
40 ((*((pWindow)->pData + (pWindow)->numCols*(R) + (C)))==1))
41
42 /**
43 @ingroup FloatDist
44 */
45
46 /**
47 @defgroup DTW Dynamic Time Warping Distance
48
49 Dynamic Time Warping Distance.
50
51 This is not really a distance since triangular inequality is
52 not respected.
53
54 The step pattern used is symmetric2.
55 Future versions of this function will provide more customization options.
56
57 */
58
59
60 /**
61 @addtogroup DTW
62 @{
63 */
64
65
66 /**
67 * @brief Dynamic Time Warping distance
68 * @param[in] pDistance Distance matrix (Query rows * Template columns)
69 * @param[in] pWindow Windowing matrix (can be NULL if no windowing used)
70 * @param[out] pDTW Temporary cost buffer (same size)
71 * @param[out] distance Distance
72 * @return ARM_MATH_ARGUMENT_ERROR in case no path can be found with window constraint
73 *
74 * @par Windowing matrix
75 *
76 * The windowing matrix is used to impose some
77 * constraints on the search for a path.
78 * The algorithm will run faster (smaller search
79 * path) but may not be able to find a solution.
80 *
81 * The distance matrix must be initialized only
82 * where the windowing matrix is containing 1.
83 * Thus, use of a window also decreases the number
84 * of distances which must be computed.
85 */
arm_dtw_distance_f32(const arm_matrix_instance_f32 * pDistance,const arm_matrix_instance_q7 * pWindow,arm_matrix_instance_f32 * pDTW,float32_t * distance)86 arm_status arm_dtw_distance_f32(const arm_matrix_instance_f32 *pDistance,
87 const arm_matrix_instance_q7 *pWindow,
88 arm_matrix_instance_f32 *pDTW,
89 float32_t *distance)
90 {
91 const uint32_t queryLength = pDistance -> numRows;
92 const uint32_t templateLength = pDistance -> numCols;
93 float32_t result;
94
95 float32_t *temp = pDTW->pData;
96 for(uint32_t q=0 ; q < queryLength; q++)
97 {
98 for(uint32_t t= 0; t < templateLength; t++)
99 {
100 *temp++ = F32_MAX;
101 }
102 }
103
104 pDTW->pData[0] = E(pDistance,0,0);
105 for(uint32_t q = 1; q < queryLength; q++)
106 {
107 if (!WIN(q,0))
108 {
109 continue;
110 }
111 E(pDTW,q,0) = E(pDTW,q-1,0) + E(pDistance,q,0);
112 }
113
114 for(uint32_t t = 1; t < templateLength; t++)
115 {
116 if (!WIN(0,t))
117 {
118 continue;
119 }
120 E(pDTW,0,t) = E(pDTW,0,t-1) + E(pDistance,0,t);
121 }
122
123
124 for(uint32_t q = 1; q < queryLength; q++)
125 {
126 for(uint32_t t = 1; t < templateLength; t++)
127 {
128 if (!WIN(q,t))
129 {
130 continue;
131 }
132 E(pDTW,q,t) =
133 MIN(E(pDTW,q-1,t-1) + 2.0f * E(pDistance,q,t),
134 MIN(E(pDTW,q,t-1) + E(pDistance,q,t),
135 E(pDTW,q-1,t) + E(pDistance,q,t)));
136 }
137 }
138
139 if (E(pDTW,queryLength-1,templateLength-1) == F32_MAX)
140 {
141 return(ARM_MATH_ARGUMENT_ERROR);
142 }
143
144 result = E(pDTW,queryLength-1,templateLength-1);
145 result = result / (queryLength + templateLength);
146 *distance = result;
147
148 return(ARM_MATH_SUCCESS);
149 }
150
151 /**
152 * @} end of DTW group
153 */
154
155