boost/math/tools/signal_statistics.hpp

   1 //  (C) Copyright Nick Thompson 2018.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP
   7 #define BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP
   8
   9 #include <algorithm>
  10 #include <iterator>
  11 #include <boost/assert.hpp>
  12 #include <boost/math/tools/complex.hpp>
  13 #include <boost/math/tools/roots.hpp>
  14 #include <boost/math/statistics/univariate_statistics.hpp>
  15 #include <boost/config/header_deprecated.hpp>
  16
  17 BOOST_HEADER_DEPRECATED("<boost/math/statistics/signal_statistics.hpp>");
  18
  19 namespace boost::math::tools {
  20
  21 template<class ForwardIterator>
  22 auto absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
  23 {
  24     using std::abs;
  25     using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;
  26     BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");
  27
  28     std::sort(first, last,  [](RealOrComplex a, RealOrComplex b) { return abs(b) > abs(a); });
  29
  30
  31     decltype(abs(*first)) i = 1;
  32     decltype(abs(*first)) num = 0;
  33     decltype(abs(*first)) denom = 0;
  34     for (auto it = first; it != last; ++it)
  35     {
  36         decltype(abs(*first)) tmp = abs(*it);
  37         num += tmp*i;
  38         denom += tmp;
  39         ++i;
  40     }
  41
  42     // If the l1 norm is zero, all elements are zero, so every element is the same.
  43     if (denom == 0)
  44     {
  45         decltype(abs(*first)) zero = 0;
  46         return zero;
  47     }
  48     return ((2*num)/denom - i)/(i-1);
  49 }
  50
  51 template<class RandomAccessContainer>
  52 inline auto absolute_gini_coefficient(RandomAccessContainer & v)
  53 {
  54     return boost::math::tools::absolute_gini_coefficient(v.begin(), v.end());
  55 }
  56
  57 template<class ForwardIterator>
  58 auto sample_absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
  59 {
  60     size_t n = std::distance(first, last);
  61     return n*boost::math::tools::absolute_gini_coefficient(first, last)/(n-1);
  62 }
  63
  64 template<class RandomAccessContainer>
  65 inline auto sample_absolute_gini_coefficient(RandomAccessContainer & v)
  66 {
  67     return boost::math::tools::sample_absolute_gini_coefficient(v.begin(), v.end());
  68 }
  69
  70
  71 // The Hoyer sparsity measure is defined in:
  72 // https://arxiv.org/pdf/0811.4706.pdf
  73 template<class ForwardIterator>
  74 auto hoyer_sparsity(const ForwardIterator first, const ForwardIterator last)
  75 {
  76     using T = typename std::iterator_traits<ForwardIterator>::value_type;
  77     using std::abs;
  78     using std::sqrt;
  79     BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Hoyer sparsity requires at least two samples.");
  80
  81     if constexpr (std::is_unsigned<T>::value)
  82     {
  83         T l1 = 0;
  84         T l2 = 0;
  85         size_t n = 0;
  86         for (auto it = first; it != last; ++it)
  87         {
  88             l1 += *it;
  89             l2 += (*it)*(*it);
  90             n += 1;
  91         }
  92
  93         double rootn = sqrt(n);
  94         return (rootn - l1/sqrt(l2) )/ (rootn - 1);
  95     }
  96     else {
  97         decltype(abs(*first)) l1 = 0;
  98         decltype(abs(*first)) l2 = 0;
  99         // We wouldn't need to count the elements if it was a random access iterator,
 100         // but our only constraint is that it's a forward iterator.
 101         size_t n = 0;
 102         for (auto it = first; it != last; ++it)
 103         {
 104             decltype(abs(*first)) tmp = abs(*it);
 105             l1 += tmp;
 106             l2 += tmp*tmp;
 107             n += 1;
 108         }
 109         if constexpr (std::is_integral<T>::value)
 110         {
 111             double rootn = sqrt(n);
 112             return (rootn - l1/sqrt(l2) )/ (rootn - 1);
 113         }
 114         else
 115         {
 116             decltype(abs(*first)) rootn = sqrt(static_cast<decltype(abs(*first))>(n));
 117             return (rootn - l1/sqrt(l2) )/ (rootn - 1);
 118         }
 119     }
 120 }
 121
 122 template<class Container>
 123 inline auto hoyer_sparsity(Container const & v)
 124 {
 125     return boost::math::tools::hoyer_sparsity(v.cbegin(), v.cend());
 126 }
 127
 128
 129 template<class Container>
 130 auto oracle_snr(Container const & signal, Container const & noisy_signal)
 131 {
 132     using Real = typename Container::value_type;
 133     BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(),
 134                      "Signal and noisy_signal must be have the same number of elements.");
 135     if constexpr (std::is_integral<Real>::value)
 136     {
 137         double numerator = 0;
 138         double denominator = 0;
 139         for (size_t i = 0; i < signal.size(); ++i)
 140         {
 141             numerator += signal[i]*signal[i];
 142             denominator += (noisy_signal[i] - signal[i])*(noisy_signal[i] - signal[i]);
 143         }
 144         if (numerator == 0 && denominator == 0)
 145         {
 146             return std::numeric_limits<double>::quiet_NaN();
 147         }
 148         if (denominator == 0)
 149         {
 150             return std::numeric_limits<double>::infinity();
 151         }
 152         return numerator/denominator;
 153     }
 154     else if constexpr (boost::math::tools::is_complex_type<Real>::value)
 155
 156     {
 157         using std::norm;
 158         typename Real::value_type numerator = 0;
 159         typename Real::value_type denominator = 0;
 160         for (size_t i = 0; i < signal.size(); ++i)
 161         {
 162             numerator += norm(signal[i]);
 163             denominator += norm(noisy_signal[i] - signal[i]);
 164         }
 165         if (numerator == 0 && denominator == 0)
 166         {
 167             return std::numeric_limits<typename Real::value_type>::quiet_NaN();
 168         }
 169         if (denominator == 0)
 170         {
 171             return std::numeric_limits<typename Real::value_type>::infinity();
 172         }
 173
 174         return numerator/denominator;
 175     }
 176     else
 177     {
 178         Real numerator = 0;
 179         Real denominator = 0;
 180         for (size_t i = 0; i < signal.size(); ++i)
 181         {
 182             numerator += signal[i]*signal[i];
 183             denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]);
 184         }
 185         if (numerator == 0 && denominator == 0)
 186         {
 187             return std::numeric_limits<Real>::quiet_NaN();
 188         }
 189         if (denominator == 0)
 190         {
 191             return std::numeric_limits<Real>::infinity();
 192         }
 193
 194         return numerator/denominator;
 195     }
 196 }
 197
 198 template<class Container>
 199 auto mean_invariant_oracle_snr(Container const & signal, Container const & noisy_signal)
 200 {
 201     using Real = typename Container::value_type;
 202     BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(), "Signal and noisy signal must be have the same number of elements.");
 203
 204     Real mu = boost::math::tools::mean(signal);
 205     Real numerator = 0;
 206     Real denominator = 0;
 207     for (size_t i = 0; i < signal.size(); ++i)
 208     {
 209         Real tmp = signal[i] - mu;
 210         numerator += tmp*tmp;
 211         denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]);
 212     }
 213     if (numerator == 0 && denominator == 0)
 214     {
 215         return std::numeric_limits<Real>::quiet_NaN();
 216     }
 217     if (denominator == 0)
 218     {
 219         return std::numeric_limits<Real>::infinity();
 220     }
 221
 222     return numerator/denominator;
 223
 224 }
 225
 226 template<class Container>
 227 auto mean_invariant_oracle_snr_db(Container const & signal, Container const & noisy_signal)
 228 {
 229     using std::log10;
 230     return 10*log10(boost::math::tools::mean_invariant_oracle_snr(signal, noisy_signal));
 231 }
 232
 233
 234 // Follows the definition of SNR given in Mallat, A Wavelet Tour of Signal Processing, equation 11.16.
 235 template<class Container>
 236 auto oracle_snr_db(Container const & signal, Container const & noisy_signal)
 237 {
 238     using std::log10;
 239     return 10*log10(boost::math::tools::oracle_snr(signal, noisy_signal));
 240 }
 241
 242 // A good reference on the M2M4 estimator:
 243 // D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR estimation techniques for the AWGN channel," IEEE Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000.
 244 // A nice python implementation:
 245 // https://github.com/gnuradio/gnuradio/blob/master/gr-digital/examples/snr_estimators.py
 246 template<class ForwardIterator>
 247 auto m2m4_snr_estimator(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3)
 248 {
 249     BOOST_ASSERT_MSG(estimated_signal_kurtosis > 0, "The estimated signal kurtosis must be positive");
 250     BOOST_ASSERT_MSG(estimated_noise_kurtosis > 0, "The estimated noise kurtosis must be positive.");
 251     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
 252     using std::sqrt;
 253     if constexpr (std::is_floating_point<Real>::value || std::numeric_limits<Real>::max_exponent)
 254     {
 255         // If we first eliminate N, we obtain the quadratic equation:
 256         // (ka+kw-6)S^2 + 2M2(3-kw)S + kw*M2^2 - M4 = 0 =: a*S^2 + bs*N + cs = 0
 257         // If we first eliminate S, we obtain the quadratic equation:
 258         // (ka+kw-6)N^2 + 2M2(3-ka)N + ka*M2^2 - M4 = 0 =: a*N^2 + bn*N + cn = 0
 259         // I believe these equations are totally independent quadratics;
 260         // if one has a complex solution it is not necessarily the case that the other must also.
 261         // However, I can't prove that, so there is a chance that this does unnecessary work.
 262         // Future improvements: There are algorithms which can solve quadratics much more effectively than the naive implementation found here.
 263         // See: https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711#50065711
 264         auto [M1, M2, M3, M4] = boost::math::tools::first_four_moments(first, last);
 265         if (M4 == 0)
 266         {
 267             // The signal is constant. There is no noise:
 268             return std::numeric_limits<Real>::infinity();
 269         }
 270         // Change to notation in Pauluzzi, equation 41:
 271         auto kw = estimated_noise_kurtosis;
 272         auto ka = estimated_signal_kurtosis;
 273         // A common case, since it's the default:
 274         Real a = (ka+kw-6);
 275         Real bs = 2*M2*(3-kw);
 276         Real cs = kw*M2*M2 - M4;
 277         Real bn = 2*M2*(3-ka);
 278         Real cn = ka*M2*M2 - M4;
 279         auto [S0, S1] = boost::math::tools::quadratic_roots(a, bs, cs);
 280         if (S1 > 0)
 281         {
 282             auto N = M2 - S1;
 283             if (N > 0)
 284             {
 285                 return S1/N;
 286             }
 287             if (S0 > 0)
 288             {
 289                 N = M2 - S0;
 290                 if (N > 0)
 291                 {
 292                     return S0/N;
 293                 }
 294             }
 295         }
 296         auto [N0, N1] = boost::math::tools::quadratic_roots(a, bn, cn);
 297         if (N1 > 0)
 298         {
 299             auto S = M2 - N1;
 300             if (S > 0)
 301             {
 302                 return S/N1;
 303             }
 304             if (N0 > 0)
 305             {
 306                 S = M2 - N0;
 307                 if (S > 0)
 308                 {
 309                     return S/N0;
 310                 }
 311             }
 312         }
 313         // This happens distressingly often. It's a limitation of the method.
 314         return std::numeric_limits<Real>::quiet_NaN();
 315     }
 316     else
 317     {
 318         BOOST_ASSERT_MSG(false, "The M2M4 estimator has not been implemented for this type.");
 319         return std::numeric_limits<Real>::quiet_NaN();
 320     }
 321 }
 322
 323 template<class Container>
 324 inline auto m2m4_snr_estimator(Container const & noisy_signal,  typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3)
 325 {
 326     return m2m4_snr_estimator(noisy_signal.cbegin(), noisy_signal.cend(), estimated_signal_kurtosis, estimated_noise_kurtosis);
 327 }
 328
 329 template<class ForwardIterator>
 330 inline auto m2m4_snr_estimator_db(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3)
 331 {
 332     using std::log10;
 333     return 10*log10(m2m4_snr_estimator(first, last, estimated_signal_kurtosis, estimated_noise_kurtosis));
 334 }
 335
 336
 337 template<class Container>
 338 inline auto m2m4_snr_estimator_db(Container const & noisy_signal,  typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3)
 339 {
 340     using std::log10;
 341     return 10*log10(m2m4_snr_estimator(noisy_signal, estimated_signal_kurtosis, estimated_noise_kurtosis));
 342 }
 343
 344 }
 345 #endif