# HG changeset patch
# User Chris Cannam
# Date 1444212856 -3600
# Node ID 26daede606a8dc547a3d1669641005077aa47571
# Parent  ca658c7215a97a28f3d74e444ed3d1f4a032cc05
Add L^p norms, doc, tests

diff -r ca658c7215a9 -r 26daede606a8 maths/MathUtilities.cpp
--- a/maths/MathUtilities.cpp	Wed Oct 07 10:36:09 2015 +0100
+++ b/maths/MathUtilities.cpp	Wed Oct 07 11:14:16 2015 +0100
@@ -20,6 +20,7 @@
 #include <vector>
 #include <cmath>
 
+using namespace std;
 
 double MathUtilities::mod(double x, double y)
 {
@@ -56,7 +57,7 @@
     *ANorm = a;
 }
 
-double MathUtilities::getAlphaNorm( const std::vector <double> &data, int alpha )
+double MathUtilities::getAlphaNorm( const vector <double> &data, int alpha )
 {
     int i;
     int len = data.size();
@@ -66,7 +67,6 @@
     for( i = 0; i < len; i++)
     {
 	temp = data[ i ];
-		
 	a  += ::pow( fabs(temp), double(alpha) );
     }
     a /= ( double )len;
@@ -88,9 +88,9 @@
 {
     if (len == 0) return 0;
     
-    std::vector<double> scratch;
+    vector<double> scratch;
     for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
-    std::sort(scratch.begin(), scratch.end());
+    sort(scratch.begin(), scratch.end());
 
     int middle = len/2;
     if (len % 2 == 0) {
@@ -126,7 +126,7 @@
     return retVal;
 }
 
-double MathUtilities::mean(const std::vector<double> &src,
+double MathUtilities::mean(const vector<double> &src,
                            int start,
                            int count)
 {
@@ -197,7 +197,7 @@
 	return index;
 }
 
-int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
+int MathUtilities::getMax( const vector<double> & data, double* pMax )
 {
 	int index = 0;
 	int i;
@@ -286,7 +286,7 @@
     }
 }
 
-void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
+void MathUtilities::normalise(vector<double> &data, NormaliseType type)
 {
     switch (type) {
 
@@ -317,20 +317,46 @@
     }
 }
 
-void MathUtilities::adaptiveThreshold(std::vector<double> &data)
+double MathUtilities::getLpNorm(const vector<double> &data, int p)
+{
+    double tot = 0.0;
+    for (int i = 0; i < int(data.size()); ++i) {
+        tot += abs(pow(data[i], p));
+    }
+    return pow(tot, 1.0 / p);
+}
+
+vector<double> MathUtilities::normaliseLp(const vector<double> &data,
+                                               int p,
+                                               double threshold)
+{
+    int n = int(data.size());
+    if (n == 0 || p == 0) return data;
+    double norm = getLpNorm(data, p);
+    if (norm < threshold) {
+        return vector<double>(n, 1.0 / pow(n, 1.0 / p)); // unit vector
+    }
+    vector<double> out(n);
+    for (int i = 0; i < n; ++i) {
+        out[i] = data[i] / norm;
+    }
+    return out;
+}
+    
+void MathUtilities::adaptiveThreshold(vector<double> &data)
 {
     int sz = int(data.size());
     if (sz == 0) return;
 
-    std::vector<double> smoothed(sz);
+    vector<double> smoothed(sz);
 	
     int p_pre = 8;
     int p_post = 7;
 
     for (int i = 0; i < sz; ++i) {
 
-        int first = std::max(0,      i - p_pre);
-        int last  = std::min(sz - 1, i + p_post);
+        int first = max(0,      i - p_pre);
+        int last  = min(sz - 1, i + p_post);
 
         smoothed[i] = mean(data, first, last - first + 1);
     }
diff -r ca658c7215a9 -r 26daede606a8 maths/MathUtilities.h
--- a/maths/MathUtilities.h	Wed Oct 07 10:36:09 2015 +0100
+++ b/maths/MathUtilities.h	Wed Oct 07 11:14:16 2015 +0100
@@ -73,15 +73,22 @@
      */
     static double mod( double x, double y);
 
+    /**
+     * The alpha norm is the alpha'th root of the mean alpha'th power
+     * magnitude. For example if alpha = 2 this corresponds to the RMS
+     * of the input data, and when alpha = 1 this is the mean
+     * magnitude.
+     */
     static void	  getAlphaNorm(const double *data, int len, int alpha, double* ANorm);
+
+    /**
+     * The alpha norm is the alpha'th root of the mean alpha'th power
+     * magnitude. For example if alpha = 2 this corresponds to the RMS
+     * of the input data, and when alpha = 1 this is the mean
+     * magnitude.
+     */
     static double getAlphaNorm(const std::vector <double> &data, int alpha );
 
-    static void   circShift( double* data, int length, int shift);
-
-    static int	  getMax( double* data, int length, double* max = 0 );
-    static int	  getMax( const std::vector<double> &data, double* max = 0 );
-    static int    compareInt(const void * a, const void * b);
-
     enum NormaliseType {
         NormaliseNone,
         NormaliseUnitSum,
@@ -95,11 +102,34 @@
                           NormaliseType n = NormaliseUnitMax);
 
     /**
+     * Calculate the L^p norm of a vector. Equivalent to MATLAB's
+     * norm(data, p).
+     */
+    static double getLpNorm(const std::vector<double> &data,
+                            int p);
+
+    /**
+     * Normalise a vector by dividing through by its L^p norm. If the
+     * norm is below the given threshold, the unit vector for that
+     * norm is returned. p may be 0, in which case no normalisation
+     * happens and the data is returned unchanged.
+     */
+    static std::vector<double> normaliseLp(const std::vector<double> &data,
+                                           int p,
+                                           double threshold = 1e-6);
+    
+    /**
      * Threshold the input/output vector data against a moving-mean
      * average filter.
      */
     static void adaptiveThreshold(std::vector<double> &data);
 
+    static void   circShift( double* data, int length, int shift);
+
+    static int	  getMax( double* data, int length, double* max = 0 );
+    static int	  getMax( const std::vector<double> &data, double* max = 0 );
+    static int    compareInt(const void * a, const void * b);
+
     /** 
      * Return true if x is 2^n for some integer n >= 0.
      */
diff -r ca658c7215a9 -r 26daede606a8 tests/TestMathUtilities.cpp
--- a/tests/TestMathUtilities.cpp	Wed Oct 07 10:36:09 2015 +0100
+++ b/tests/TestMathUtilities.cpp	Wed Oct 07 11:14:16 2015 +0100
@@ -3,12 +3,15 @@
 #include "maths/MathUtilities.h"
 
 #include <cmath>
+#include <iostream>
 
 #define BOOST_TEST_DYN_LINK
 #define BOOST_TEST_MAIN
 
 #include <boost/test/unit_test.hpp>
 
+using namespace std;
+
 BOOST_AUTO_TEST_SUITE(TestMathUtilities)
 
 BOOST_AUTO_TEST_CASE(round)
@@ -34,7 +37,7 @@
     BOOST_CHECK_EQUAL(MathUtilities::mean(d0, 4), 1.5);
     double d1[] = { -2.6 };
     BOOST_CHECK_EQUAL(MathUtilities::mean(d1, 1), -2.6);
-    std::vector<double> v;
+    vector<double> v;
     v.push_back(0);
     v.push_back(4);
     v.push_back(3);
@@ -148,6 +151,98 @@
     BOOST_CHECK_EQUAL(MathUtilities::gcd(37, 18), 1);
 }
 
+BOOST_AUTO_TEST_CASE(getAlphaNorm1)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    double out = MathUtilities::getAlphaNorm(in, 1);
+    double expected = 2.6250;
+    double thresh = 1e-4;
+    BOOST_CHECK_SMALL(out - expected, thresh);
+}
+
+BOOST_AUTO_TEST_CASE(getAlphaNorm2)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    double out = MathUtilities::getAlphaNorm(in, 2);
+    double expected = 3.0516;
+    double thresh = 1e-4;
+    BOOST_CHECK_SMALL(out - expected, thresh);
+}
+
+BOOST_AUTO_TEST_CASE(getL1Norm)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    double out = MathUtilities::getLpNorm(in, 1);
+    // L1 norm is the sum of magnitudes
+    double expected = 10.5;
+    double thresh = 1e-5;
+    BOOST_CHECK_SMALL(out - expected, thresh);
+}
+
+BOOST_AUTO_TEST_CASE(getL2Norm)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    double out = MathUtilities::getLpNorm(in, 2);
+    // L2 norm is the sqrt of sum of squared magnitudes
+    double expected = sqrt(37.25);
+    double thresh = 1e-5;
+    BOOST_CHECK_SMALL(out - expected, thresh);
+}
+
+BOOST_AUTO_TEST_CASE(getL3Norm)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    double out = MathUtilities::getLpNorm(in, 3);
+    double expected = 5.3875;
+    double thresh = 1e-4;
+    BOOST_CHECK_SMALL(out - expected, thresh);
+}
+
+BOOST_AUTO_TEST_CASE(normaliseL1)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    vector<double> expected { -0.095238, 0.142857, 0.285714, 0.476190 };
+    vector<double> out = MathUtilities::normaliseLp(in, 1);
+    double thresh = 1e-5;
+    for (int i = 0; i < int(out.size()); ++i) {
+        BOOST_CHECK_SMALL(out[i] - expected[i], thresh);
+    }
+    out = MathUtilities::normaliseLp({ 0, 0, 0, 0 }, 1);
+    for (int i = 0; i < int(out.size()); ++i) {
+        BOOST_CHECK_EQUAL(out[i], 0.25);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(normaliseL2)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    vector<double> expected { -0.16385, 0.24577, 0.49154, 0.81923 };
+    vector<double> out = MathUtilities::normaliseLp(in, 2);
+    double thresh = 1e-5;
+    for (int i = 0; i < int(out.size()); ++i) {
+        BOOST_CHECK_SMALL(out[i] - expected[i], thresh);
+    }
+    out = MathUtilities::normaliseLp({ 0, 0, 0, 0 }, 2);
+    for (int i = 0; i < int(out.size()); ++i) {
+        BOOST_CHECK_EQUAL(out[i], 0.5);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(normaliseL3)
+{
+    vector<double> in { -1, 1.5, 3, 5 };
+    vector<double> expected { -0.18561, 0.27842, 0.55684, 0.92807 };
+    vector<double> out = MathUtilities::normaliseLp(in, 3);
+    double thresh = 1e-5;
+    for (int i = 0; i < int(out.size()); ++i) {
+        BOOST_CHECK_SMALL(out[i] - expected[i], thresh);
+    }
+    out = MathUtilities::normaliseLp({ 0, 0, 0, 0 }, 3);
+    for (int i = 0; i < int(out.size()); ++i) {
+        BOOST_CHECK_SMALL(out[i] - 0.62996, 1e-5);
+    }
+}
+
 BOOST_AUTO_TEST_SUITE_END()