qm-dsp  1.8
KLDivergence.cpp
Go to the documentation of this file.
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
00002 
00003 /*
00004     QM DSP Library
00005 
00006     Centre for Digital Music, Queen Mary, University of London.
00007     This file copyright 2008 QMUL
00008 
00009     This program is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU General Public License as
00011     published by the Free Software Foundation; either version 2 of the
00012     License, or (at your option) any later version.  See the file
00013     COPYING included with this distribution for more information.
00014 */
00015 
00016 #include "KLDivergence.h"
00017 
00018 #include <cmath>
00019 
00020 double KLDivergence::distanceGaussian(const vector<double> &m1,
00021                                       const vector<double> &v1,
00022                                       const vector<double> &m2,
00023                                       const vector<double> &v2)
00024 {
00025     int sz = m1.size();
00026 
00027     double d = -2.0 * sz;
00028     double small = 1e-20;
00029 
00030     for (int k = 0; k < sz; ++k) {
00031 
00032         double kv1 = v1[k] + small;
00033         double kv2 = v2[k] + small;
00034         double km = (m1[k] - m2[k]) + small;
00035 
00036         d += kv1 / kv2 + kv2 / kv1;
00037         d += km * (1.0 / kv1 + 1.0 / kv2) * km;
00038     }
00039 
00040     d /= 2.0;
00041 
00042     return d;
00043 }
00044 
00045 double KLDivergence::distanceDistribution(const vector<double> &d1,
00046                                           const vector<double> &d2,
00047                                           bool symmetrised)
00048 {
00049     int sz = d1.size();
00050 
00051     double d = 0;
00052     double small = 1e-20;
00053     
00054     for (int i = 0; i < sz; ++i) {
00055         d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
00056     }
00057 
00058     if (symmetrised) {
00059         d += distanceDistribution(d2, d1, false);
00060     }
00061 
00062     return d;
00063 }
00064