First version of the dynamic time warping algorithm. (b46a2f86) · Commits · Wenxi XU / Cmsis Dsp

Include/dsp/distance_functions.h

+46 −0

Original line number	Diff line number	Diff line
		@@ -36,6 +36,7 @@
		#include "dsp/statistics_functions.h"
		#include "dsp/basic_math_functions.h"
		#include "dsp/fast_math_functions.h"
		#include "dsp/matrix_functions.h"

		#ifdef __cplusplus
		extern "C"
		@@ -332,8 +333,53 @@ float32_t arm_sokalsneath_distance(const uint32_t pA, const uint32_t pB, uint3

		float32_t arm_yule_distance(const uint32_t pA, const uint32_t pB, uint32_t numberOfBools);

		typedef enum
		{
		ARM_DTW_SAKOE_CHIBA_WINDOW = 1,
		/ARM_DTW_ITAKURA_WINDOW = 2,/
		ARM_DTW_SLANTED_BAND_WINDOW = 3
		} arm_dtw_window;

		/**
		* @brief Window for dynamic time warping computation
		* @param[in] windowType Type of window
		* @param[in] windowSize Window size
		* @param[in,out] pWindow Window
		* @return Error if window type not recognized
		*
		*/
		arm_status arm_dtw_init_window_q7(const arm_dtw_window windowType,
		const int32_t windowSize,
		arm_matrix_instance_q7 *pWindow);

		/**
		* @brief Dynamic Time Warping distance
		* @param[in] pDistance Distance matrix (Query rows * Template columns)
		* @param[in] pWindow Windowing (can be NULL if no windowing used)
		* @param[out] pDTW Temporary cost buffer (same size)
		* @param[out] distance Distance
		* @return Error in case no path can be found with window constraint
		*
		*/

		arm_status arm_dtw_distance_f32(const arm_matrix_instance_f32 *pDistance,
		const arm_matrix_instance_q7 *pWindow,
		arm_matrix_instance_f32 *pDTW,
		float32_t *distance);


		/**
		* @brief Mapping between query and template
		* @param[in] pDTW Cost matrix (Query rows * Template columns)
		* @param[out] pPath Warping path in cost matrix 2*nb rows + nb columns)
		* @param[out] pathLength Length of path in number of points
		* @return none
		*
		*/

		void arm_dtw_path_f32(const arm_matrix_instance_f32 *pDTW,
		int16_t *pPath,
		uint32_t *pathLength);
		#ifdef __cplusplus
		}
		#endif

Include/dsp/utils.h

+3 −0

Original line number	Diff line number	Diff line
		@@ -40,6 +40,9 @@ extern "C"

		#define INDEX_MASK 0x0000003F

		#define MIN(x,y) ((x) < (y) ? (x) : (y))

		#define MAX(x,y) ((x) > (y) ? (x) : (y))

		#define SQ(x) ((x) * (x))

Source/DistanceFunctions/Config.cmake

+3 −0

Original line number	Diff line number	Diff line
		@@ -25,6 +25,9 @@ target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_russellrao_distance.c)
		target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_sokalmichener_distance.c)
		target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_sokalsneath_distance.c)
		target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_yule_distance.c)
		target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_dtw_distance_f32.c)
		target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_dtw_path_f32.c)
		target_sources(CMSISDSP PRIVATE DistanceFunctions/arm_dtw_init_window_q7.c)


		target_include_directories(CMSISDSP PRIVATE "DistanceFunctions")

Source/DistanceFunctions/DistanceFunctions.c

+3 −0

Original line number	Diff line number	Diff line
		@@ -49,3 +49,6 @@
		#include "arm_sokalmichener_distance.c"
		#include "arm_sokalsneath_distance.c"
		#include "arm_yule_distance.c"
		#include "arm_dtw_distance_f32.c"
		#include "arm_dtw_path_f32.c"
		#include "arm_dtw_init_window_q7.c"
		No newline at end of file

Source/DistanceFunctions/arm_dtw_distance_f32.c

0 → 100644

+153 −0

Original line number	Diff line number	Diff line

		/* ----------------------------------------------------------------------
		* Project: CMSIS DSP Library
		* Title: arm_braycurtis_distance_f32.c
		* Description: Bray-Curtis distance between two vectors
		*
		* $Date: 23 April 2021
		* $Revision: V1.9.0
		*
		* Target Processor: Cortex-M and Cortex-A cores
		* -------------------------------------------------------------------- */
		/*
		* Copyright (C) 2010-2021 ARM Limited or its affiliates. All rights reserved.
		*
		* SPDX-License-Identifier: Apache-2.0
		*
		* Licensed under the Apache License, Version 2.0 (the License); you may
		* not use this file except in compliance with the License.
		* You may obtain a copy of the License at
		*
		* www.apache.org/licenses/LICENSE-2.0
		*
		* Unless required by applicable law or agreed to in writing, software
		* distributed under the License is distributed on an AS IS BASIS, WITHOUT
		* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
		* See the License for the specific language governing permissions and
		* limitations under the License.
		*/

		#include "dsp/distance_functions.h"
		#include <limits.h>
		#include <math.h>

		#define E(MAT,R,C) \
		(((MAT)->pData + (MAT)->numCols(R) + (C)))

		#define WIN(R,C) \
		((pWindow == NULL) ? 1 : \
		((((pWindow)->pData + (pWindow)->numCols(R) + (C)))==1))

		/**
		@ingroup FloatDist
		*/

		/**
		@defgroup DTW Dynamic Time Warping Distance

		Dynamic Time Warping Distance.

		This is not really a distance since triangular inequality is
		not respected.

		The step pattern used is symmetric2.
		Future versions of this function will provide more customization options.

		*/


		/**
		@addtogroup DTW
		@{
		*/


		/**
		* @brief Dynamic Time Warping distance
		* @param[in] pDistance Distance matrix (Query rows * Template columns)
		* @param[in] pWindow Windowing matrix (can be NULL if no windowing used)
		* @param[out] pDTW Temporary cost buffer (same size)
		* @param[out] distance Distance
		* @return ARM_MATH_ARGUMENT_ERROR in case no path can be found with window constraint
		*
		* @par Windowing matrix
		*
		* The windowing matrix is used to impose some
		* constraints on the search for a path.
		* The algorithm will run faster (smaller search
		* path) but may not be able to find a solution.
		*
		* The distance matrix must be initialized only
		* where the windowing matrix is containing 1.
		* Thus, use of a window also decreases the number
		* of distances which must be computed.
		*/
		arm_status arm_dtw_distance_f32(const arm_matrix_instance_f32 *pDistance,
		const arm_matrix_instance_q7 *pWindow,
		arm_matrix_instance_f32 *pDTW,
		float32_t *distance)
		{
		const uint32_t queryLength = pDistance -> numRows;
		const uint32_t templateLength = pDistance -> numCols;
		float32_t result;

		float32_t *temp = pDTW->pData;
		for(uint32_t q=0 ; q < queryLength; q++)
		{
		for(uint32_t t= 0; t < templateLength; t++)
		{
		*temp++ = F32_MAX;
		}
		}

		pDTW->pData[0] = E(pDistance,0,0);
		for(uint32_t q = 1; q < queryLength; q++)
		{
		if (!WIN(q,0))
		{
		continue;
		}
		E(pDTW,q,0) = E(pDTW,q-1,0) + E(pDistance,q,0);
		}

		for(uint32_t t = 1; t < templateLength; t++)
		{
		if (!WIN(0,t))
		{
		continue;
		}
		E(pDTW,0,t) = E(pDTW,0,t-1) + E(pDistance,0,t);
		}


		for(uint32_t q = 1; q < queryLength; q++)
		{
		for(uint32_t t = 1; t < templateLength; t++)
		{
		if (!WIN(q,t))
		{
		continue;
		}
		E(pDTW,q,t) =
		MIN(E(pDTW,q-1,t-1) + 2.0f * E(pDistance,q,t),
		MIN(E(pDTW,q,t-1) + E(pDistance,q,t),
		E(pDTW,q-1,t) + E(pDistance,q,t)));
		}
		}

		if (E(pDTW,queryLength-1,templateLength-1) == F32_MAX)
		{
		return(ARM_MATH_ARGUMENT_ERROR);
		}

		result = E(pDTW,queryLength-1,templateLength-1);
		result = result / (queryLength + templateLength);
		*distance = result;

		return(ARM_MATH_SUCCESS);
		}

		/**
		* @} end of DTW group
		*/

Admin message