src/stats.c

   1 /* libSoX effect: stats   (c) 2009 robs@users.sourceforge.net
   2  *
   3  * This library is free software; you can redistribute it and/or modify it
   4  * under the terms of the GNU Lesser General Public License as published by
   5  * the Free Software Foundation; either version 2.1 of the License, or (at
   6  * your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful, but
   9  * WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
  11  * General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public License
  14  * along with this library; if not, write to the Free Software Foundation,
  15  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  16  */
  17
  18 #include "sox_i.h"
  19 #include <ctype.h>
  20 #include <string.h>
  21
  22 typedef struct {
  23   int       scale_bits, hex_bits;
  24   double    time_constant, scale;
  25
  26   double    last, sigma_x, sigma_x2, avg_sigma_x2, min_sigma_x2, max_sigma_x2;
  27   double    min, max, mult, min_run, min_runs, max_run, max_runs;
  28   off_t     num_samples, tc_samples, min_count, max_count;
  29   uint32_t  mask;
  30 } priv_t;
  31
  32 static int getopts(sox_effect_t * effp, int argc, char **argv)
  33 {
  34   priv_t * p = (priv_t *)effp->priv;
  35   int c;
  36   lsx_getopt_t optstate;
  37   lsx_getopt_init(argc, argv, "+x:b:w:s:", NULL, lsx_getopt_flag_none, 1, &optstate);
  38
  39   p->time_constant = .05;
  40   p->scale = 1;
  41   while ((c = lsx_getopt(&optstate)) != -1) switch (c) {
  42     GETOPT_NUMERIC(optstate, 'x', hex_bits      ,  2 , 32)
  43     GETOPT_NUMERIC(optstate, 'b', scale_bits    ,  2 , 32)
  44     GETOPT_NUMERIC(optstate, 'w', time_constant ,  .01 , 10)
  45     GETOPT_NUMERIC(optstate, 's', scale         ,  -99, 99)
  46     default: lsx_fail("invalid option `-%c'", optstate.opt); return lsx_usage(effp);
  47   }
  48   if (p->hex_bits)
  49     p->scale_bits = p->hex_bits;
  50   return optstate.ind != argc? lsx_usage(effp) : SOX_SUCCESS;
  51 }
  52
  53 static int start(sox_effect_t * effp)
  54 {
  55   priv_t * p = (priv_t *)effp->priv;
  56
  57   p->last = 0;
  58   p->mult = exp((-1 / p->time_constant / effp->in_signal.rate));
  59   p->tc_samples = 5 * p->time_constant * effp->in_signal.rate + .5;
  60   p->sigma_x = p->sigma_x2 = p->avg_sigma_x2 = p->max_sigma_x2 = 0;
  61   p->min = p->min_sigma_x2 = 2;
  62   p->max = -p->min;
  63   p->num_samples = 0;
  64   p->mask = 0;
  65   return SOX_SUCCESS;
  66 }
  67
  68 static int flow(sox_effect_t * effp, const sox_sample_t * ibuf,
  69     sox_sample_t * obuf, size_t * ilen, size_t * olen)
  70 {
  71   priv_t * p = (priv_t *)effp->priv;
  72   size_t len = *ilen = *olen = min(*ilen, *olen);
  73   memcpy(obuf, ibuf, len * sizeof(*obuf));
  74
  75   for (; len--; ++ibuf, ++p->num_samples) {
  76     double d = SOX_SAMPLE_TO_FLOAT_64BIT(*ibuf,);
  77
  78     if (d < p->min)
  79       p->min = d, p->min_count = 1, p->min_run = 1, p->min_runs = 0;
  80     else if (d == p->min) {
  81       ++p->min_count;
  82       p->min_run = d == p->last? p->min_run + 1 : 1;
  83     }
  84     else if (p->last == p->min)
  85       p->min_runs += sqr(p->min_run);
  86
  87     if (d > p->max)
  88       p->max = d, p->max_count = 1, p->max_run = 1, p->max_runs = 0;
  89     else if (d == p->max) {
  90       ++p->max_count;
  91       p->max_run = d == p->last? p->max_run + 1 : 1;
  92     }
  93     else if (p->last == p->max)
  94       p->max_runs += sqr(p->max_run);
  95
  96     p->sigma_x += d;
  97     p->sigma_x2 += sqr(d);
  98     p->avg_sigma_x2 = p->avg_sigma_x2 * p->mult + (1 - p->mult) * sqr(d);
  99
 100     if (p->num_samples >= p->tc_samples) {
 101       if (p->avg_sigma_x2 > p->max_sigma_x2)
 102         p->max_sigma_x2 = p->avg_sigma_x2;
 103       if (p->avg_sigma_x2 < p->min_sigma_x2)
 104         p->min_sigma_x2 = p->avg_sigma_x2;
 105     }
 106     p->last = d;
 107     p->mask |= *ibuf;
 108   }
 109   return SOX_SUCCESS;
 110 }
 111
 112 static int drain(sox_effect_t * effp, sox_sample_t * obuf, size_t * olen)
 113 {
 114   priv_t * p = (priv_t *)effp->priv;
 115
 116   if (p->last == p->min)
 117     p->min_runs += sqr(p->min_run);
 118   if (p->last == p->max)
 119     p->max_runs += sqr(p->max_run);
 120
 121   (void)obuf, *olen = 0;
 122   return SOX_SUCCESS;
 123 }
 124
 125 static unsigned bit_depth(uint32_t mask, double min, double max, unsigned * x)
 126 {
 127   SOX_SAMPLE_LOCALS;
 128   unsigned result = 32, dummy = 0;
 129
 130   for (; result && !(mask & 1); --result, mask >>= 1);
 131   if (x)
 132     *x = result;
 133   mask = SOX_FLOAT_64BIT_TO_SAMPLE(max, dummy);
 134   if (min < 0)
 135     mask |= ~(SOX_FLOAT_64BIT_TO_SAMPLE(min, dummy) << 1);
 136   for (; result && !(mask & SOX_SAMPLE_MIN); --result, mask <<= 1);
 137   return result;
 138 }
 139
 140 static void output(priv_t const * p, double x)
 141 {
 142   if (p->scale_bits) {
 143     unsigned mult = 1 << (p->scale_bits - 1);
 144     int i;
 145     x = floor(x * mult + .5);
 146     i = min(x, mult - 1.);
 147     if (p->hex_bits)
 148       if (x < 0) {
 149         char buf[30];
 150         sprintf(buf, "%x", -i);
 151         fprintf(stderr, " %*c%s", 9 - (int)strlen(buf), '-', buf);
 152       }
 153       else fprintf(stderr, " %9x", i);
 154     else fprintf(stderr, " %9i", i);
 155   }
 156   else fprintf(stderr, " %9.*f", fabs(p->scale) < 10 ? 6 : 5, p->scale * x);
 157 }
 158
 159 static int stop(sox_effect_t * effp)
 160 {
 161   priv_t * p = (priv_t *)effp->priv;
 162
 163   if (!effp->flow) {
 164     double min_runs = 0, max_count = 0, min = 2, max = -2, max_sigma_x = 0, sigma_x = 0, sigma_x2 = 0, min_sigma_x2 = 2, max_sigma_x2 = 0, avg_peak = 0;
 165     off_t num_samples = 0, min_count = 0, max_runs = 0;
 166     uint32_t mask = 0;
 167     unsigned b1, b2, i, n = effp->flows > 1 ? effp->flows : 0;
 168
 169     for (i = 0; i < effp->flows; ++i) {
 170       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 171       min = min(min, q->min);
 172       max = max(max, q->max);
 173       if (q->num_samples < q->tc_samples)
 174         q->min_sigma_x2 = q->max_sigma_x2 = q->sigma_x2 / q->num_samples;
 175       min_sigma_x2 = min(min_sigma_x2, q->min_sigma_x2);
 176       max_sigma_x2 = max(max_sigma_x2, q->max_sigma_x2);
 177       sigma_x += q->sigma_x;
 178       sigma_x2 += q->sigma_x2;
 179       num_samples += q->num_samples;
 180       mask |= q->mask;
 181       if (fabs(q->sigma_x) > fabs(max_sigma_x))
 182         max_sigma_x = q->sigma_x;
 183       min_count += q->min_count;
 184       min_runs += q->min_runs;
 185       max_count += q->max_count;
 186       max_runs += q->max_runs;
 187       avg_peak += max(-q->min, q->max);
 188     }
 189     avg_peak /= effp->flows;
 190
 191     if (!num_samples) {
 192       lsx_warn("no audio");
 193       return SOX_SUCCESS;
 194     }
 195
 196     if (n == 2)
 197       fprintf(stderr, "             Overall     Left      Right\n");
 198     else if (n) {
 199       fprintf(stderr, "             Overall");
 200       for (i = 0; i < n; ++i)
 201         fprintf(stderr, "     Ch%-3i", i + 1);
 202       fprintf(stderr, "\n");
 203     }
 204
 205     fprintf(stderr, "DC offset ");
 206     output(p, max_sigma_x / p->num_samples);
 207     for (i = 0; i < n; ++i) {
 208       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 209       output(p, q->sigma_x / q->num_samples);
 210     }
 211
 212     fprintf(stderr, "\nMin level ");
 213     output(p, min);
 214     for (i = 0; i < n; ++i) {
 215       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 216       output(p, q->min);
 217     }
 218
 219     fprintf(stderr, "\nMax level ");
 220     output(p, max);
 221     for (i = 0; i < n; ++i) {
 222       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 223       output(p, q->max);
 224     }
 225
 226     fprintf(stderr, "\nPk lev dB %10.2f", linear_to_dB(max(-min, max)));
 227     for (i = 0; i < n; ++i) {
 228       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 229       fprintf(stderr, "%10.2f", linear_to_dB(max(-q->min, q->max)));
 230     }
 231
 232     fprintf(stderr, "\nRMS lev dB%10.2f", linear_to_dB(sqrt(sigma_x2 / num_samples)));
 233     for (i = 0; i < n; ++i) {
 234       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 235       fprintf(stderr, "%10.2f", linear_to_dB(sqrt(q->sigma_x2 / q->num_samples)));
 236     }
 237
 238     fprintf(stderr, "\nRMS Pk dB %10.2f", linear_to_dB(sqrt(max_sigma_x2)));
 239     for (i = 0; i < n; ++i) {
 240       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 241       fprintf(stderr, "%10.2f", linear_to_dB(sqrt(q->max_sigma_x2)));
 242     }
 243
 244     fprintf(stderr, "\nRMS Tr dB ");
 245     if (min_sigma_x2 != 1)
 246       fprintf(stderr, "%10.2f", linear_to_dB(sqrt(min_sigma_x2)));
 247     else fprintf(stderr, "         -");
 248     for (i = 0; i < n; ++i) {
 249       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 250       if (q->min_sigma_x2 != 1)
 251         fprintf(stderr, "%10.2f", linear_to_dB(sqrt(q->min_sigma_x2)));
 252       else fprintf(stderr, "         -");
 253     }
 254
 255     if (effp->flows > 1)
 256       fprintf(stderr, "\nCrest factor       -");
 257     else fprintf(stderr, "\nCrest factor %7.2f", sigma_x2 ? avg_peak / sqrt(sigma_x2 / num_samples) : 1);
 258     for (i = 0; i < n; ++i) {
 259       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 260       fprintf(stderr, "%10.2f", q->sigma_x2? max(-q->min, q->max) / sqrt(q->sigma_x2 / q->num_samples) : 1);
 261     }
 262
 263     fprintf(stderr, "\nFlat factor%9.2f", linear_to_dB((min_runs + max_runs) / (min_count + max_count)));
 264     for (i = 0; i < n; ++i) {
 265       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 266       fprintf(stderr, " %9.2f", linear_to_dB((q->min_runs + q->max_runs) / (q->min_count + q->max_count)));
 267     }
 268
 269     fprintf(stderr, "\nPk count   %9s", lsx_sigfigs3((min_count + max_count) / effp->flows));
 270     for (i = 0; i < n; ++i) {
 271       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 272       fprintf(stderr, " %9s", lsx_sigfigs3((double)(q->min_count + q->max_count)));
 273     }
 274
 275     b1 = bit_depth(mask, min, max, &b2);
 276     fprintf(stderr, "\nBit-depth      %2u/%-2u", b1, b2);
 277     for (i = 0; i < n; ++i) {
 278       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 279       b1 = bit_depth(q->mask, q->min, q->max, &b2);
 280       fprintf(stderr, "     %2u/%-2u", b1, b2);
 281     }
 282
 283     fprintf(stderr, "\nNum samples%9s", lsx_sigfigs3((double)p->num_samples));
 284     fprintf(stderr, "\nLength s   %9.3f", p->num_samples / effp->in_signal.rate);
 285     fprintf(stderr, "\nScale max ");
 286     output(p, 1.);
 287     fprintf(stderr, "\nWindow s   %9.3f", p->time_constant);
 288     fprintf(stderr, "\n");
 289   }
 290   return SOX_SUCCESS;
 291 }
 292
 293 sox_effect_handler_t const * lsx_stats_effect_fn(void)
 294 {
 295   static sox_effect_handler_t handler = {
 296     "stats", "[-b bits|-x bits|-s scale] [-w window-time]", SOX_EFF_MODIFY,
 297     getopts, start, flow, drain, stop, NULL, sizeof(priv_t)};
 298   return &handler;
 299 }