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