c@256: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@256: 
c@256: /*
c@256:     QM DSP Library
c@256: 
c@256:     Centre for Digital Music, Queen Mary, University of London.
c@256:     This file copyright 2008 QMUL
c@309: 
c@309:     This program is free software; you can redistribute it and/or
c@309:     modify it under the terms of the GNU General Public License as
c@309:     published by the Free Software Foundation; either version 2 of the
c@309:     License, or (at your option) any later version.  See the file
c@309:     COPYING included with this distribution for more information.
c@256: */
c@256: 
c@256: #include "KLDivergence.h"
c@256: 
c@258: #include <cmath>
c@258: 
cannam@493: using std::vector;
cannam@493: 
c@258: double KLDivergence::distanceGaussian(const vector<double> &m1,
c@258:                                       const vector<double> &v1,
c@258:                                       const vector<double> &m2,
c@258:                                       const vector<double> &v2)
c@256: {
c@256:     int sz = m1.size();
c@256: 
c@256:     double d = -2.0 * sz;
c@299:     double small = 1e-20;
c@256: 
c@256:     for (int k = 0; k < sz; ++k) {
c@299: 
c@299:         double kv1 = v1[k] + small;
c@299:         double kv2 = v2[k] + small;
c@299:         double km = (m1[k] - m2[k]) + small;
c@299: 
c@299:         d += kv1 / kv2 + kv2 / kv1;
c@299:         d += km * (1.0 / kv1 + 1.0 / kv2) * km;
c@256:     }
c@256: 
c@256:     d /= 2.0;
c@256: 
c@256:     return d;
c@256: }
c@258: 
c@258: double KLDivergence::distanceDistribution(const vector<double> &d1,
c@258:                                           const vector<double> &d2,
c@258:                                           bool symmetrised)
c@258: {
c@258:     int sz = d1.size();
c@258: 
c@258:     double d = 0;
c@258:     double small = 1e-20;
c@258:     
c@258:     for (int i = 0; i < sz; ++i) {
c@258:         d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
c@258:     }
c@258: 
c@258:     if (symmetrised) {
c@258:         d += distanceDistribution(d2, d1, false);
c@258:     }
c@258: 
c@258:     return d;
c@258: }
c@258: