libs/qm-dsp/maths/KLDivergence.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP Library
   5
   6     Centre for Digital Music, Queen Mary, University of London.
   7     This file copyright 2008 QMUL
   8
   9     This program is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU General Public License as
  11     published by the Free Software Foundation; either version 2 of the
  12     License, or (at your option) any later version.  See the file
  13     COPYING included with this distribution for more information.
  14 */
  15
  16 #include "KLDivergence.h"
  17
  18 #include <cmath>
  19
  20 double KLDivergence::distanceGaussian(const vector<double> &m1,
  21                                       const vector<double> &v1,
  22                                       const vector<double> &m2,
  23                                       const vector<double> &v2)
  24 {
  25     int sz = m1.size();
  26
  27     double d = -2.0 * sz;
  28     double small = 1e-20;
  29
  30     for (int k = 0; k < sz; ++k) {
  31
  32         double kv1 = v1[k] + small;
  33         double kv2 = v2[k] + small;
  34         double km = (m1[k] - m2[k]) + small;
  35
  36         d += kv1 / kv2 + kv2 / kv1;
  37         d += km * (1.0 / kv1 + 1.0 / kv2) * km;
  38     }
  39
  40     d /= 2.0;
  41
  42     return d;
  43 }
  44
  45 double KLDivergence::distanceDistribution(const vector<double> &d1,
  46                                           const vector<double> &d2,
  47                                           bool symmetrised)
  48 {
  49     int sz = d1.size();
  50
  51     double d = 0;
  52     double small = 1e-20;
  53
  54     for (int i = 0; i < sz; ++i) {
  55         d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
  56     }
  57
  58     if (symmetrised) {
  59         d += distanceDistribution(d2, d1, false);
  60     }
  61
  62     return d;
  63 }
  64