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   min = -fmax(fabs(min), fabs(max));
 134   mask = SOX_FLOAT_64BIT_TO_SAMPLE(min, dummy) << 1;
 135   for (; result && (mask & SOX_SAMPLE_MIN); --result, mask <<= 1);
 136   return result;
 137 }
 138
 139 static void output(priv_t const * p, double x)
 140 {
 141   if (p->scale_bits) {
 142     unsigned mult = 1 << (p->scale_bits - 1);
 143     int i;
 144     x = floor(x * mult + .5);
 145     i = min(x, mult - 1.);
 146     if (p->hex_bits)
 147       if (x < 0) {
 148         char buf[30];
 149         sprintf(buf, "%x", -i);
 150         fprintf(stderr, " %*c%s", 9 - (int)strlen(buf), '-', buf);
 151       }
 152       else fprintf(stderr, " %9x", i);
 153     else fprintf(stderr, " %9i", i);
 154   }
 155   else fprintf(stderr, " %9.*f", fabs(p->scale) < 10 ? 6 : 5, p->scale * x);
 156 }
 157
 158 static int stop(sox_effect_t * effp)
 159 {
 160   priv_t * p = (priv_t *)effp->priv;
 161
 162   if (!effp->flow) {
 163     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;
 164     off_t num_samples = 0, min_count = 0, max_runs = 0;
 165     uint32_t mask = 0;
 166     unsigned b1, b2, i, n = effp->flows > 1 ? effp->flows : 0;
 167
 168     for (i = 0; i < effp->flows; ++i) {
 169       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 170       min = min(min, q->min);
 171       max = max(max, q->max);
 172       if (q->num_samples < q->tc_samples)
 173         q->min_sigma_x2 = q->max_sigma_x2 = q->sigma_x2 / q->num_samples;
 174       min_sigma_x2 = min(min_sigma_x2, q->min_sigma_x2);
 175       max_sigma_x2 = max(max_sigma_x2, q->max_sigma_x2);
 176       sigma_x += q->sigma_x;
 177       sigma_x2 += q->sigma_x2;
 178       num_samples += q->num_samples;
 179       mask |= q->mask;
 180       if (fabs(q->sigma_x) > fabs(max_sigma_x))
 181         max_sigma_x = q->sigma_x;
 182       min_count += q->min_count;
 183       min_runs += q->min_runs;
 184       max_count += q->max_count;
 185       max_runs += q->max_runs;
 186       avg_peak += max(-q->min, q->max);
 187     }
 188     avg_peak /= effp->flows;
 189
 190     if (!num_samples) {
 191       lsx_warn("no audio");
 192       return SOX_SUCCESS;
 193     }
 194
 195     if (n == 2)
 196       fprintf(stderr, "             Overall     Left      Right\n");
 197     else if (n) {
 198       fprintf(stderr, "             Overall");
 199       for (i = 0; i < n; ++i)
 200         fprintf(stderr, "     Ch%-3i", i + 1);
 201       fprintf(stderr, "\n");
 202     }
 203
 204     fprintf(stderr, "DC offset ");
 205     output(p, max_sigma_x / p->num_samples);
 206     for (i = 0; i < n; ++i) {
 207       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 208       output(p, q->sigma_x / q->num_samples);
 209     }
 210
 211     fprintf(stderr, "\nMin level ");
 212     output(p, min);
 213     for (i = 0; i < n; ++i) {
 214       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 215       output(p, q->min);
 216     }
 217
 218     fprintf(stderr, "\nMax level ");
 219     output(p, max);
 220     for (i = 0; i < n; ++i) {
 221       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 222       output(p, q->max);
 223     }
 224
 225     fprintf(stderr, "\nPk lev dB %10.2f", linear_to_dB(max(-min, max)));
 226     for (i = 0; i < n; ++i) {
 227       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 228       fprintf(stderr, "%10.2f", linear_to_dB(max(-q->min, q->max)));
 229     }
 230
 231     fprintf(stderr, "\nRMS lev dB%10.2f", linear_to_dB(sqrt(sigma_x2 / num_samples)));
 232     for (i = 0; i < n; ++i) {
 233       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 234       fprintf(stderr, "%10.2f", linear_to_dB(sqrt(q->sigma_x2 / q->num_samples)));
 235     }
 236
 237     fprintf(stderr, "\nRMS Pk dB %10.2f", linear_to_dB(sqrt(max_sigma_x2)));
 238     for (i = 0; i < n; ++i) {
 239       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 240       fprintf(stderr, "%10.2f", linear_to_dB(sqrt(q->max_sigma_x2)));
 241     }
 242
 243     fprintf(stderr, "\nRMS Tr dB ");
 244     if (min_sigma_x2 != 1)
 245       fprintf(stderr, "%10.2f", linear_to_dB(sqrt(min_sigma_x2)));
 246     else fprintf(stderr, "         -");
 247     for (i = 0; i < n; ++i) {
 248       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 249       if (q->min_sigma_x2 != 1)
 250         fprintf(stderr, "%10.2f", linear_to_dB(sqrt(q->min_sigma_x2)));
 251       else fprintf(stderr, "         -");
 252     }
 253
 254     if (effp->flows > 1)
 255       fprintf(stderr, "\nCrest factor       -");
 256     else fprintf(stderr, "\nCrest factor %7.2f", sigma_x2 ? avg_peak / sqrt(sigma_x2 / num_samples) : 1);
 257     for (i = 0; i < n; ++i) {
 258       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 259       fprintf(stderr, "%10.2f", q->sigma_x2? max(-q->min, q->max) / sqrt(q->sigma_x2 / q->num_samples) : 1);
 260     }
 261
 262     fprintf(stderr, "\nFlat factor%9.2f", linear_to_dB((min_runs + max_runs) / (min_count + max_count)));
 263     for (i = 0; i < n; ++i) {
 264       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 265       fprintf(stderr, " %9.2f", linear_to_dB((q->min_runs + q->max_runs) / (q->min_count + q->max_count)));
 266     }
 267
 268     fprintf(stderr, "\nPk count   %9s", lsx_sigfigs3((min_count + max_count) / effp->flows));
 269     for (i = 0; i < n; ++i) {
 270       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 271       fprintf(stderr, " %9s", lsx_sigfigs3((double)(q->min_count + q->max_count)));
 272     }
 273
 274     b1 = bit_depth(mask, min, max, &b2);
 275     fprintf(stderr, "\nBit-depth      %2u/%-2u", b1, b2);
 276     for (i = 0; i < n; ++i) {
 277       priv_t * q = (priv_t *)(effp - effp->flow + i)->priv;
 278       b1 = bit_depth(q->mask, q->min, q->max, &b2);
 279       fprintf(stderr, "     %2u/%-2u", b1, b2);
 280     }
 281
 282     fprintf(stderr, "\nNum samples%9s", lsx_sigfigs3((double)p->num_samples));
 283     fprintf(stderr, "\nLength s   %9.3f", p->num_samples / effp->in_signal.rate);
 284     fprintf(stderr, "\nScale max ");
 285     output(p, 1.);
 286     fprintf(stderr, "\nWindow s   %9.3f", p->time_constant);
 287     fprintf(stderr, "\n");
 288   }
 289   return SOX_SUCCESS;
 290 }
 291
 292 sox_effect_handler_t const * lsx_stats_effect_fn(void)
 293 {
 294   static sox_effect_handler_t handler = {
 295     "stats", "[-b bits|-x bits|-s scale] [-w window-time]", SOX_EFF_MODIFY,
 296     getopts, start, flow, drain, stop, NULL, sizeof(priv_t)};
 297   return &handler;
 298 }