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
|