annotate maths/KLDivergence.cpp @ 309:d5014ab8b0e5

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