cannam@31: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@31: 
cannam@31: /*
cannam@31:     QM DSP Library
cannam@31: 
cannam@31:     Centre for Digital Music, Queen Mary, University of London.
cannam@31:     This file copyright 2008 QMUL
cannam@31:     All rights reserved.
cannam@31: */
cannam@31: 
cannam@31: #include "KLDivergence.h"
cannam@31: 
cannam@33: #include <cmath>
cannam@33: 
cannam@33: double KLDivergence::distanceGaussian(const vector<double> &m1,
cannam@33:                                       const vector<double> &v1,
cannam@33:                                       const vector<double> &m2,
cannam@33:                                       const vector<double> &v2)
cannam@31: {
cannam@31:     int sz = m1.size();
cannam@31: 
cannam@31:     double d = -2.0 * sz;
cannam@74:     double small = 1e-20;
cannam@31: 
cannam@31:     for (int k = 0; k < sz; ++k) {
cannam@74: 
cannam@74:         double kv1 = v1[k] + small;
cannam@74:         double kv2 = v2[k] + small;
cannam@74:         double km = (m1[k] - m2[k]) + small;
cannam@74: 
cannam@74:         d += kv1 / kv2 + kv2 / kv1;
cannam@74:         d += km * (1.0 / kv1 + 1.0 / kv2) * km;
cannam@31:     }
cannam@31: 
cannam@31:     d /= 2.0;
cannam@31: 
cannam@31:     return d;
cannam@31: }
cannam@33: 
cannam@33: double KLDivergence::distanceDistribution(const vector<double> &d1,
cannam@33:                                           const vector<double> &d2,
cannam@33:                                           bool symmetrised)
cannam@33: {
cannam@33:     int sz = d1.size();
cannam@33: 
cannam@33:     double d = 0;
cannam@33:     double small = 1e-20;
cannam@33:     
cannam@33:     for (int i = 0; i < sz; ++i) {
cannam@33:         d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
cannam@33:     }
cannam@33: 
cannam@33:     if (symmetrised) {
cannam@33:         d += distanceDistribution(d2, d1, false);
cannam@33:     }
cannam@33: 
cannam@33:     return d;
cannam@33: }
cannam@33: