arm_kullback_leibler_f64.c 1.9 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. /* ----------------------------------------------------------------------
  2. * Project: CMSIS DSP Library
  3. * Title: arm_logsumexp_f64.c
  4. * Description: LogSumExp
  5. *
  6. *
  7. * Target Processor: Cortex-M and Cortex-A cores
  8. * -------------------------------------------------------------------- */
  9. /*
  10. * Copyright (C) 2010-2020 ARM Limited or its affiliates. All rights reserved.
  11. *
  12. * SPDX-License-Identifier: Apache-2.0
  13. *
  14. * Licensed under the Apache License, Version 2.0 (the License); you may
  15. * not use this file except in compliance with the License.
  16. * You may obtain a copy of the License at
  17. *
  18. * www.apache.org/licenses/LICENSE-2.0
  19. *
  20. * Unless required by applicable law or agreed to in writing, software
  21. * distributed under the License is distributed on an AS IS BASIS, WITHOUT
  22. * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  23. * See the License for the specific language governing permissions and
  24. * limitations under the License.
  25. */
  26. #include "arm_math.h"
  27. #include <limits.h>
  28. #include <math.h>
  29. /**
  30. * @addtogroup groupStats
  31. * @{
  32. */
  33. /**
  34. * @brief Kullback-Leibler
  35. *
  36. * @param[in] *pSrcA points to an array of input values for probaility distribution A.
  37. * @param[in] *pSrcB points to an array of input values for probaility distribution B.
  38. * @param[in] blockSize number of samples in the input array.
  39. * @return Kullback-Leibler divergence D(A || B)
  40. *
  41. */
  42. float64_t arm_kullback_leibler_f64(const float64_t * pSrcA, const float64_t * pSrcB, uint32_t blockSize)
  43. {
  44. const float64_t *pInA, *pInB;
  45. uint32_t blkCnt;
  46. float64_t accum, pA,pB;
  47. pInA = pSrcA;
  48. pInB = pSrcB;
  49. blkCnt = blockSize;
  50. accum = 0.0f;
  51. while(blkCnt > 0)
  52. {
  53. pA = *pInA++;
  54. pB = *pInB++;
  55. accum += pA * log(pB / pA);
  56. blkCnt--;
  57. }
  58. return(-accum);
  59. }
  60. /**
  61. * @} end of groupStats group
  62. */