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