1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Signal Statistics</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../statistics.html" title="Chapter 6. Statistics"> 9<link rel="prev" href="bivariate_statistics.html" title="Bivariate Statistics"> 10<link rel="next" href="anderson_darling.html" title="The Anderson-Darling Test"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="bivariate_statistics.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../statistics.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="anderson_darling.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.signal_statistics"></a><a class="link" href="signal_statistics.html" title="Signal Statistics">Signal Statistics</a> 28</h2></div></div></div> 29<h4> 30<a name="math_toolkit.signal_statistics.h0"></a> 31 <span class="phrase"><a name="math_toolkit.signal_statistics.synopsis"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.synopsis">Synopsis</a> 32 </h4> 33<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">statistics</span><span class="special">/</span><span class="identifier">signal_statistics</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 34 35<span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span> <span class="special">{</span> 36 37 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 38 <span class="keyword">auto</span> <span class="identifier">absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">Container</span> <span class="special">&</span> <span class="identifier">c</span><span class="special">);</span> 39 40 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">ForwardIterator</span><span class="special">></span> 41 <span class="keyword">auto</span> <span class="identifier">absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">ForwardIterator</span> <span class="identifier">first</span><span class="special">,</span> <span class="identifier">ForwardIterator</span> <span class="identifier">last</span><span class="special">);</span> 42 43 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 44 <span class="keyword">auto</span> <span class="identifier">sample_absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">Container</span> <span class="special">&</span> <span class="identifier">c</span><span class="special">);</span> 45 46 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">ForwardIterator</span><span class="special">></span> 47 <span class="keyword">auto</span> <span class="identifier">sample_absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">ForwardIterator</span> <span class="identifier">first</span><span class="special">,</span> <span class="identifier">ForwardIterator</span> <span class="identifier">last</span><span class="special">);</span> 48 49 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 50 <span class="keyword">auto</span> <span class="identifier">hoyer_sparsity</span><span class="special">(</span><span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">c</span><span class="special">);</span> 51 52 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">ForwardIterator</span><span class="special">></span> 53 <span class="keyword">auto</span> <span class="identifier">hoyer_sparsity</span><span class="special">(</span><span class="identifier">ForwardIterator</span> <span class="identifier">first</span><span class="special">,</span> <span class="identifier">ForwardIterator</span> <span class="identifier">last</span><span class="special">);</span> 54 55 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 56 <span class="keyword">auto</span> <span class="identifier">oracle_snr</span><span class="special">(</span><span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">signal</span><span class="special">,</span> <span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">noisy_signal</span><span class="special">);</span> 57 58 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 59 <span class="keyword">auto</span> <span class="identifier">oracle_snr_db</span><span class="special">(</span><span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">signal</span><span class="special">,</span> <span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">noisy_signal</span><span class="special">);</span> 60 61 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">ForwardIterator</span><span class="special">></span> 62 <span class="keyword">auto</span> <span class="identifier">m2m4_snr_estimator</span><span class="special">(</span><span class="identifier">ForwardIterator</span> <span class="identifier">first</span><span class="special">,</span> <span class="identifier">ForwardIterator</span> <span class="identifier">last</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(*</span><span class="identifier">first</span><span class="special">)</span> <span class="identifier">estimated_signal_kurtosis</span><span class="special">=</span><span class="number">1</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(*</span><span class="identifier">first</span><span class="special">)</span> <span class="identifier">estimated_noise_kurtosis</span><span class="special">=</span><span class="number">3</span><span class="special">);</span> 63 64 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 65 <span class="keyword">auto</span> <span class="identifier">m2m4_snr_estimator</span><span class="special">(</span><span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">noisy_signal</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Container</span><span class="special">::</span><span class="identifier">value_type</span> <span class="identifier">estimated_signal_kurtosis</span><span class="special">=</span><span class="number">1</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Container</span><span class="special">::</span><span class="identifier">value_type</span> <span class="identifier">estimate_noise_kurtosis</span><span class="special">=</span><span class="number">3</span><span class="special">);</span> 66 67 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">ForwardIterator</span><span class="special">></span> 68 <span class="keyword">auto</span> <span class="identifier">m2m4_snr_estimator_db</span><span class="special">(</span><span class="identifier">ForwardIterator</span> <span class="identifier">first</span><span class="special">,</span> <span class="identifier">ForwardIterator</span> <span class="identifier">last</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(*</span><span class="identifier">first</span><span class="special">)</span> <span class="identifier">estimated_signal_kurtosis</span><span class="special">=</span><span class="number">1</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(*</span><span class="identifier">first</span><span class="special">)</span> <span class="identifier">estimated_noise_kurtosis</span><span class="special">=</span><span class="number">3</span><span class="special">);</span> 69 70 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Container</span><span class="special">></span> 71 <span class="keyword">auto</span> <span class="identifier">m2m4_snr_estimator_db</span><span class="special">(</span><span class="identifier">Container</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">noisy_signal</span><span class="special">,</span><span class="keyword">typename</span> <span class="identifier">Container</span><span class="special">::</span><span class="identifier">value_type</span> <span class="identifier">estimated_signal_kurtosis</span><span class="special">=</span><span class="number">1</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Container</span><span class="special">::</span><span class="identifier">value_type</span> <span class="identifier">estimate_noise_kurtosis</span><span class="special">=</span><span class="number">3</span><span class="special">);</span> 72 73<span class="special">}</span> 74</pre> 75<h4> 76<a name="math_toolkit.signal_statistics.h1"></a> 77 <span class="phrase"><a name="math_toolkit.signal_statistics.description"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.description">Description</a> 78 </h4> 79<p> 80 The file <code class="computeroutput"><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">statistics</span><span class="special">/</span><span class="identifier">signal_statistics</span><span class="special">.</span><span class="identifier">hpp</span></code> is a 81 set of facilities for computing quantities commonly used in signal analysis. 82 </p> 83<p> 84 Our examples use <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span></code> to 85 hold the data, but this not required. In general, you can store your data in 86 an Eigen array, and Armadillo vector, <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">array</span></code>, 87 and for many of the routines, a <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">forward_list</span></code>. 88 These routines are usable in float, double, long double, and Boost.Multiprecision 89 precision, as well as their complex extensions whenever the computation is 90 well-defined. 91 </p> 92<h4> 93<a name="math_toolkit.signal_statistics.h2"></a> 94 <span class="phrase"><a name="math_toolkit.signal_statistics.absolute_gini_coefficient"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.absolute_gini_coefficient">Absolute 95 Gini Coefficient</a> 96 </h4> 97<p> 98 The Gini coefficient, first used to measure wealth inequality, is also one 99 of the best measures of the sparsity of an expansion in a basis. A sparse expansion 100 has most of its norm concentrated in just a few coefficients, making the connection 101 with wealth inequality obvious. See <a href="https://arxiv.org/pdf/0811.4706.pdf" target="_top">Hurley 102 and Rickard</a> for details. However, for measuring sparsity, the phase 103 of the numbers is irrelevant, so we provide the <code class="computeroutput"><span class="identifier">absolute_gini_coefficient</span></code>: 104 </p> 105<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">sample_absolute_gini_coefficient</span><span class="special">;</span> 106<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">absolute_gini_coefficient</span><span class="special">;</span> 107<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">>></span> <span class="identifier">v</span><span class="special">{{</span><span class="number">0</span><span class="special">,</span><span class="number">1</span><span class="special">},</span> <span class="special">{</span><span class="number">0</span><span class="special">,</span><span class="number">0</span><span class="special">},</span> <span class="special">{</span><span class="number">0</span><span class="special">,</span><span class="number">0</span><span class="special">},</span> <span class="special">{</span><span class="number">0</span><span class="special">,</span><span class="number">0</span><span class="special">}};</span> 108<span class="keyword">double</span> <span class="identifier">abs_gini</span> <span class="special">=</span> <span class="identifier">sample_absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">v</span><span class="special">);</span> 109<span class="comment">// now abs_gini = 1; maximally unequal</span> 110 111<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">>></span> <span class="identifier">w</span><span class="special">{{</span><span class="number">0</span><span class="special">,</span><span class="number">1</span><span class="special">},</span> <span class="special">{</span><span class="number">1</span><span class="special">,</span><span class="number">0</span><span class="special">},</span> <span class="special">{</span><span class="number">0</span><span class="special">,-</span><span class="number">1</span><span class="special">},</span> <span class="special">{-</span><span class="number">1</span><span class="special">,</span><span class="number">0</span><span class="special">}};</span> 112<span class="identifier">abs_gini</span> <span class="special">=</span> <span class="identifier">absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">w</span><span class="special">);</span> 113<span class="comment">// now abs_gini = 0; every element of the vector has equal magnitude</span> 114 115<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">u</span><span class="special">{-</span><span class="number">1</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="special">-</span><span class="number">1</span><span class="special">};</span> 116<span class="identifier">abs_gini</span> <span class="special">=</span> <span class="identifier">absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">u</span><span class="special">);</span> 117<span class="comment">// now abs_gini = 0</span> 118<span class="comment">// Alternative call useful for computing over subset of the input:</span> 119<span class="identifier">abs_gini</span> <span class="special">=</span> <span class="identifier">absolute_gini_coefficient</span><span class="special">(</span><span class="identifier">u</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> <span class="identifier">u</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">1</span><span class="special">);</span> 120</pre> 121<p> 122 The sample Gini coefficient returns unity for a vector which has only one nonzero 123 coefficient. The population Gini coefficient of a vector with one non-zero 124 element is dependent on the length of the input. 125 </p> 126<p> 127 The sample Gini coefficient lacks one desirable property of the population 128 Gini coefficient, namely that "cloning" a vector has the same Gini 129 coefficient; though cloning holds to very high accuracy with the sample Gini 130 coefficient and can easily be recovered by a rescaling. 131 </p> 132<p> 133 If sorting the input data is too much expense for a sparsity measure (is it 134 going to be perfect anyway?), consider calculating the Hoyer sparsity instead. 135 </p> 136<h4> 137<a name="math_toolkit.signal_statistics.h3"></a> 138 <span class="phrase"><a name="math_toolkit.signal_statistics.hoyer_sparsity"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.hoyer_sparsity">Hoyer 139 Sparsity</a> 140 </h4> 141<p> 142 The Hoyer sparsity measures a normalized ratio of the ℓ<sup>1</sup> and ℓ<sup>2</sup> norms. 143 As the name suggests, it is used to measure the sparsity of an expansion in 144 some basis. 145 </p> 146<p> 147 The Hoyer sparsity computes (√<span class="emphasis"><em>N</em></span> - ℓ<sup>1</sup>(v)/ℓ<sup>2</sup>(v))/(√N 148 -1). For details, see <a href="http://www.jmlr.org/papers/volume5/hoyer04a/hoyer04a.pdf" target="_top">Hoyer</a> 149 as well as <a href="https://arxiv.org/pdf/0811.4706.pdf" target="_top">Hurley and Rickard</a>. 150 </p> 151<p> 152 A few special cases will serve to clarify the intended use: If <span class="emphasis"><em>v</em></span> 153 has only one nonzero coefficient, the Hoyer sparsity attains its maxima of 154 1. If the coefficients of <span class="emphasis"><em>v</em></span> all have the same magnitude, 155 then the Hoyer sparsity attains its minima of zero. If the elements of <span class="emphasis"><em>v</em></span> 156 are uniformly distributed on an interval [0, <span class="emphasis"><em>b</em></span>], then 157 the Hoyer sparsity is approximately 0.133. 158 </p> 159<p> 160 Usage: 161 </p> 162<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span> <span class="identifier">v</span><span class="special">{</span><span class="number">1</span><span class="special">,</span><span class="number">0</span><span class="special">,</span><span class="number">0</span><span class="special">};</span> 163<span class="identifier">Real</span> <span class="identifier">hs</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">hoyer_sparsity</span><span class="special">(</span><span class="identifier">v</span><span class="special">);</span> 164<span class="comment">// hs = 1</span> 165<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span> <span class="identifier">v</span><span class="special">{</span><span class="number">1</span><span class="special">,-</span><span class="number">1</span><span class="special">,</span><span class="number">1</span><span class="special">};</span> 166<span class="identifier">Real</span> <span class="identifier">hs</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">hoyer_sparsity</span><span class="special">(</span><span class="identifier">v</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> <span class="identifier">v</span><span class="special">.</span><span class="identifier">end</span><span class="special">());</span> 167<span class="comment">// hs = 0</span> 168</pre> 169<p> 170 The container must be forward iterable and the contents are not modified. Accepts 171 real, complex, and integer inputs. If the input is an integral type, the output 172 is a double precision float. 173 </p> 174<h4> 175<a name="math_toolkit.signal_statistics.h4"></a> 176 <span class="phrase"><a name="math_toolkit.signal_statistics.oracle_signal_to_noise_ratio"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.oracle_signal_to_noise_ratio">Oracle 177 Signal-to-noise ratio</a> 178 </h4> 179<p> 180 The function <code class="computeroutput"><span class="identifier">oracle_snr</span></code> computes 181 the ratio ‖ <span class="emphasis"><em>s</em></span> ‖<sub>2</sub><sup>2</sup> / ‖ <span class="emphasis"><em>s</em></span> 182 - <span class="emphasis"><em>x</em></span> ‖<sub>2</sub><sup>2</sup>, where <span class="emphasis"><em>s</em></span> is signal 183 and <span class="emphasis"><em>x</em></span> is a noisy signal. The function <code class="computeroutput"><span class="identifier">oracle_snr_db</span></code> 184 computes 10<code class="computeroutput"><span class="identifier">log</span></code><sub>10</sub>(‖ 185 <span class="emphasis"><em>s</em></span> ‖<sup>2</sup> / ‖ <span class="emphasis"><em>s</em></span> - <span class="emphasis"><em>x</em></span> 186 ‖<sup>2</sup>). The functions are so named because in general, one does not know 187 how to decompose a real signal <span class="emphasis"><em>x</em></span> into <span class="emphasis"><em>s</em></span> 188 + <span class="emphasis"><em>w</em></span> and as such <span class="emphasis"><em>s</em></span> is regarded as 189 oracle information. Hence this function is mainly useful for unit testing other 190 SNR estimators. 191 </p> 192<p> 193 Usage: 194 </p> 195<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">signal</span><span class="special">(</span><span class="number">500</span><span class="special">,</span> <span class="number">3.2</span><span class="special">);</span> 196<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">noisy_signal</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 197<span class="comment">// fill 'noisy_signal' signal + noise</span> 198<span class="keyword">double</span> <span class="identifier">snr_db</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">oracle_snr_db</span><span class="special">(</span><span class="identifier">signal</span><span class="special">,</span> <span class="identifier">noisy_signal</span><span class="special">);</span> 199<span class="keyword">double</span> <span class="identifier">snr</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">oracle_snr</span><span class="special">(</span><span class="identifier">signal</span><span class="special">,</span> <span class="identifier">noisy_signal</span><span class="special">);</span> 200</pre> 201<p> 202 The input can be real, complex, or integral. Integral inputs produce double 203 precision floating point outputs. The input data is not modified and must satisfy 204 the requirements of a <code class="computeroutput"><span class="identifier">RandomAccessContainer</span></code>. 205 </p> 206<h4> 207<a name="math_toolkit.signal_statistics.h5"></a> 208 <span class="phrase"><a name="math_toolkit.signal_statistics.m_sub_2_m_sub_4_snr_estimation"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.m_sub_2_m_sub_4_snr_estimation"><span class="emphasis"><em>M</em></span><sub>2</sub><span class="emphasis"><em>M</em></span><sub>4</sub> SNR 209 Estimation</a> 210 </h4> 211<p> 212 Estimates the SNR of a noisy signal via the <span class="emphasis"><em>M</em></span><sub>2</sub><span class="emphasis"><em>M</em></span><sub>4</sub> method. 213 See <a href="https://doi.org/10.1109/26.871393" target="_top">Pauluzzi and N.C. Beaulieu</a> 214 and <a href="https://doi.org/10.1109/ISIT.1994.394869" target="_top">Matzner and Englberger</a> 215 for details. 216 </p> 217<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">noisy_signal</span><span class="special">(</span><span class="number">512</span><span class="special">);</span> 218<span class="comment">// fill noisy_signal with data contaminated by Gaussian white noise:</span> 219<span class="keyword">double</span> <span class="identifier">est_snr_db</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">m2m4_snr_estimator_db</span><span class="special">(</span><span class="identifier">noisy_signal</span><span class="special">);</span> 220</pre> 221<p> 222 The <span class="emphasis"><em>M</em></span><sub>2</sub><span class="emphasis"><em>M</em></span><sub>4</sub> SNR estimator is an "in-service" 223 estimator, meaning that the estimate is made using the noisy, data-bearing 224 signal, and does not require a background estimate. This estimator has been 225 found to be work best between roughly -3 and 15db, tending to overestimate 226 the noise below -3db, and underestimate the noise above 15db. See <a href="https://www.mdpi.com/2078-2489/8/3/75/pdf" target="_top">Xue 227 et al</a> for details. 228 </p> 229<p> 230 The <span class="emphasis"><em>M</em></span><sub>2</sub><span class="emphasis"><em>M</em></span><sub>4</sub> SNR estimator, by default, 231 assumes that the kurtosis of the signal is 1 and the kurtosis of the noise 232 is 3, the latter corresponding to Gaussian noise. These parameters, however, 233 can be overridden: 234 </p> 235<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">noisy_signal</span><span class="special">(</span><span class="number">512</span><span class="special">);</span> 236<span class="comment">// fill noisy_signal with the data:</span> 237<span class="keyword">double</span> <span class="identifier">signal_kurtosis</span> <span class="special">=</span> <span class="number">1.5</span><span class="special">;</span> 238<span class="comment">// Noise is assumed to follow Laplace distribution, which has kurtosis of 6:</span> 239<span class="keyword">double</span> <span class="identifier">noise_kurtosis</span> <span class="special">=</span> <span class="number">6</span><span class="special">;</span> 240<span class="keyword">double</span> <span class="identifier">est_snr</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">statistics</span><span class="special">::</span><span class="identifier">m2m4_snr_estimator_db</span><span class="special">(</span><span class="identifier">noisy_signal</span><span class="special">,</span> <span class="identifier">signal_kurtosis</span><span class="special">,</span> <span class="identifier">noise_kurtosis</span><span class="special">);</span> 241</pre> 242<p> 243 Now, technically the method is a "blind SNR estimator", meaning that 244 the no <span class="emphasis"><em>a-priori</em></span> information about the signal is required 245 to use the method. However, the performance of the method is <span class="emphasis"><em>vastly</em></span> 246 better if you can come up with a better estimate of the signal and noise kurtosis. 247 How can we do this? Suppose we know that the SNR is much greater than 1. Then 248 we can estimate the signal kurtosis simply by using the noisy signal kurtosis. 249 If the SNR is much less than one, this method breaks down as the noisy signal 250 kurtosis will tend to the noise kurtosis-though in this limit we have an excellent 251 estimator of the noise kurtosis! In addition, if you have a model of what your 252 signal should look like, you can precompute the signal kurtosis. For example, 253 sinusoids have a kurtosis of 1.5. See <a href="http://www.jcomputers.us/vol8/jcp0808-21.pdf" target="_top">here</a> 254 for a study which uses estimates of this sort to improve the performance of 255 the <span class="emphasis"><em>M</em></span><sub>2</sub><span class="emphasis"><em>M</em></span><sub>4</sub> estimator. 256 </p> 257<p> 258 <span class="emphasis"><em>Nota bene</em></span>: The traditional definition of SNR is <span class="emphasis"><em>not</em></span> 259 mean invariant. By this we mean that if a constant is added to every sample 260 of a signal, the SNR is changed. For example, adding DC bias to a signal changes 261 its SNR. For most use cases, this is really not what you intend; for example 262 a signal consisting of zeros plus Gaussian noise has an SNR of zero, whereas 263 a signal with a constant DC bias and random Gaussian noise might have a very 264 large SNR. 265 </p> 266<p> 267 The <span class="emphasis"><em>M</em></span><sub>2</sub><span class="emphasis"><em>M</em></span><sub>4</sub> SNR estimator is computed 268 from mean-invariant quantities, and hence it should really be compared to the 269 mean-invariant SNR. 270 </p> 271<p> 272 <span class="emphasis"><em>Nota bene</em></span>: This computation requires the solution of a 273 system of quadratic equations involving the noise kurtosis, the signal kurtosis, 274 and the second and fourth moments of the data. There is no guarantee that a 275 solution of this system exists for all value of these parameters, in fact nonexistence 276 can easily be demonstrated for certain data. If there is no solution to the 277 system, then failure is communicated by returning NaNs. This happens distressingly 278 often; if a user is aware of any blind SNR estimators which do not suffer from 279 this drawback, please open a github ticket and let us know. 280 </p> 281<p> 282 The author has not managed to fully characterize the conditions under which 283 a real solution with <span class="emphasis"><em>S > 0</em></span> and <span class="emphasis"><em>N >0</em></span> 284 exists. However, a very intuitive example demonstrates why nonexistence can 285 occur. Suppose the signal and noise kurtosis are equal. Then the method has 286 no way to distinguish between the signal and the noise, and the solution is 287 non-unique. 288 </p> 289<h4> 290<a name="math_toolkit.signal_statistics.h6"></a> 291 <span class="phrase"><a name="math_toolkit.signal_statistics.references"></a></span><a class="link" href="signal_statistics.html#math_toolkit.signal_statistics.references">References</a> 292 </h4> 293<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 294<li class="listitem"> 295 Mallat, Stephane. <span class="emphasis"><em>A wavelet tour of signal processing: the sparse 296 way.</em></span> Academic press, 2008. 297 </li> 298<li class="listitem"> 299 Hurley, Niall, and Scott Rickard. <span class="emphasis"><em>Comparing measures of sparsity.</em></span> 300 IEEE Transactions on Information Theory 55.10 (2009): 4723-4741. 301 </li> 302<li class="listitem"> 303 Jensen, Arne, and Anders la Cour-Harbo. <span class="emphasis"><em>Ripples in mathematics: 304 the discrete wavelet transform.</em></span> Springer Science & Business 305 Media, 2001. 306 </li> 307<li class="listitem"> 308 D. R. Pauluzzi and N. C. Beaulieu, <span class="emphasis"><em>A comparison of SNR estimation 309 techniques for the AWGN channel,</em></span> IEEE Trans. Communications, 310 Vol. 48, No. 10, pp. 1681-1691, 2000. 311 </li> 312<li class="listitem"> 313 Hoyer, Patrik O. <span class="emphasis"><em>Non-negative matrix factorization with sparseness 314 constraints.</em></span>, Journal of machine learning research 5.Nov (2004): 315 1457-1469. 316 </li> 317</ul></div> 318</div> 319<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 320<td align="left"></td> 321<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 322 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 323 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 324 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 325 Daryle Walker and Xiaogang Zhang<p> 326 Distributed under the Boost Software License, Version 1.0. (See accompanying 327 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 328 </p> 329</div></td> 330</tr></table> 331<hr> 332<div class="spirit-nav"> 333<a accesskey="p" href="bivariate_statistics.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../statistics.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="anderson_darling.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 334</div> 335</body> 336</html> 337
