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
|
cannam@493
|
20 using std::vector;
|
cannam@493
|
21
|
c@258
|
22 double KLDivergence::distanceGaussian(const vector<double> &m1,
|
c@258
|
23 const vector<double> &v1,
|
c@258
|
24 const vector<double> &m2,
|
c@258
|
25 const vector<double> &v2)
|
c@256
|
26 {
|
c@256
|
27 int sz = m1.size();
|
c@256
|
28
|
c@256
|
29 double d = -2.0 * sz;
|
c@299
|
30 double small = 1e-20;
|
c@256
|
31
|
c@256
|
32 for (int k = 0; k < sz; ++k) {
|
c@299
|
33
|
c@299
|
34 double kv1 = v1[k] + small;
|
c@299
|
35 double kv2 = v2[k] + small;
|
c@299
|
36 double km = (m1[k] - m2[k]) + small;
|
c@299
|
37
|
c@299
|
38 d += kv1 / kv2 + kv2 / kv1;
|
c@299
|
39 d += km * (1.0 / kv1 + 1.0 / kv2) * km;
|
c@256
|
40 }
|
c@256
|
41
|
c@256
|
42 d /= 2.0;
|
c@256
|
43
|
c@256
|
44 return d;
|
c@256
|
45 }
|
c@258
|
46
|
c@258
|
47 double KLDivergence::distanceDistribution(const vector<double> &d1,
|
c@258
|
48 const vector<double> &d2,
|
c@258
|
49 bool symmetrised)
|
c@258
|
50 {
|
c@258
|
51 int sz = d1.size();
|
c@258
|
52
|
c@258
|
53 double d = 0;
|
c@258
|
54 double small = 1e-20;
|
c@258
|
55
|
c@258
|
56 for (int i = 0; i < sz; ++i) {
|
c@258
|
57 d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
|
c@258
|
58 }
|
c@258
|
59
|
c@258
|
60 if (symmetrised) {
|
c@258
|
61 d += distanceDistribution(d2, d1, false);
|
c@258
|
62 }
|
c@258
|
63
|
c@258
|
64 return d;
|
c@258
|
65 }
|
c@258
|
66
|