src/flac/analyze.c

   1 /* flac - Command-line FLAC encoder/decoder
   2  * Copyright (C) 2000,2001,2002,2003,2004,2005,2006,2007,2008,2009  Josh Coalson
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License
   6  * as published by the Free Software Foundation; either version 2
   7  * of the License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License along
  15  * with this program; if not, write to the Free Software Foundation, Inc.,
  16  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  17  */
  18
  19 #if HAVE_CONFIG_H
  20 #  include <config.h>
  21 #endif
  22
  23 #include <errno.h>
  24 #include <math.h>
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28
  29 #include "FLAC/all.h"
  30 #include "analyze.h"
  31
  32 #include "share/compat.h"
  33
  34 typedef struct {
  35         FLAC__int32 residual;
  36         unsigned count;
  37 } pair_t;
  38
  39 typedef struct {
  40         pair_t buckets[FLAC__MAX_BLOCK_SIZE];
  41         int peak_index;
  42         unsigned nbuckets;
  43         unsigned nsamples;
  44         double sum, sos;
  45         double variance;
  46         double mean;
  47         double stddev;
  48 } subframe_stats_t;
  49
  50 static subframe_stats_t all_;
  51
  52 static void init_stats(subframe_stats_t *stats);
  53 static void update_stats(subframe_stats_t *stats, FLAC__int32 residual, unsigned incr);
  54 static void compute_stats(subframe_stats_t *stats);
  55 static FLAC__bool dump_stats(const subframe_stats_t *stats, const char *filename);
  56
  57 void flac__analyze_init(analysis_options aopts)
  58 {
  59         if(aopts.do_residual_gnuplot) {
  60                 init_stats(&all_);
  61         }
  62 }
  63
  64 void flac__analyze_frame(const FLAC__Frame *frame, unsigned frame_number, FLAC__uint64 frame_offset, unsigned frame_bytes, analysis_options aopts, FILE *fout)
  65 {
  66         const unsigned channels = frame->header.channels;
  67         char outfilename[1024];
  68         subframe_stats_t stats;
  69         unsigned i, channel, partitions;
  70
  71         /* do the human-readable part first */
  72         fprintf(fout, "frame=%u\toffset=%" PRIu64 "\tbits=%u\tblocksize=%u\tsample_rate=%u\tchannels=%u\tchannel_assignment=%s\n", frame_number, frame_offset, frame_bytes*8, frame->header.blocksize, frame->header.sample_rate, channels, FLAC__ChannelAssignmentString[frame->header.channel_assignment]);
  73         for(channel = 0; channel < channels; channel++) {
  74                 const FLAC__Subframe *subframe = frame->subframes+channel;
  75                 const FLAC__bool is_rice2 = subframe->data.fixed.entropy_coding_method.type == FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE2;
  76                 const unsigned pesc = is_rice2? FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE2_ESCAPE_PARAMETER : FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE_ESCAPE_PARAMETER;
  77                 fprintf(fout, "\tsubframe=%u\twasted_bits=%u\ttype=%s", channel, subframe->wasted_bits, FLAC__SubframeTypeString[subframe->type]);
  78                 switch(subframe->type) {
  79                         case FLAC__SUBFRAME_TYPE_CONSTANT:
  80                                 fprintf(fout, "\tvalue=%d\n", subframe->data.constant.value);
  81                                 break;
  82                         case FLAC__SUBFRAME_TYPE_FIXED:
  83                                 FLAC__ASSERT(subframe->data.fixed.entropy_coding_method.type <= FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE2);
  84                                 fprintf(fout, "\torder=%u\tresidual_type=%s\tpartition_order=%u\n", subframe->data.fixed.order, is_rice2? "RICE2":"RICE", subframe->data.fixed.entropy_coding_method.data.partitioned_rice.order);
  85                                 for(i = 0; i < subframe->data.fixed.order; i++)
  86                                         fprintf(fout, "\t\twarmup[%u]=%d\n", i, subframe->data.fixed.warmup[i]);
  87                                 partitions = (1u << subframe->data.fixed.entropy_coding_method.data.partitioned_rice.order);
  88                                 for(i = 0; i < partitions; i++) {
  89                                         unsigned parameter = subframe->data.fixed.entropy_coding_method.data.partitioned_rice.contents->parameters[i];
  90                                         if(parameter == pesc)
  91                                                 fprintf(fout, "\t\tparameter[%u]=ESCAPE, raw_bits=%u\n", i, subframe->data.fixed.entropy_coding_method.data.partitioned_rice.contents->raw_bits[i]);
  92                                         else
  93                                                 fprintf(fout, "\t\tparameter[%u]=%u\n", i, parameter);
  94                                 }
  95                                 if(aopts.do_residual_text) {
  96                                         for(i = 0; i < frame->header.blocksize-subframe->data.fixed.order; i++)
  97                                                 fprintf(fout, "\t\tresidual[%u]=%d\n", i, subframe->data.fixed.residual[i]);
  98                                 }
  99                                 break;
 100                         case FLAC__SUBFRAME_TYPE_LPC:
 101                                 FLAC__ASSERT(subframe->data.lpc.entropy_coding_method.type <= FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE2);
 102                                 fprintf(fout, "\torder=%u\tqlp_coeff_precision=%u\tquantization_level=%d\tresidual_type=%s\tpartition_order=%u\n", subframe->data.lpc.order, subframe->data.lpc.qlp_coeff_precision, subframe->data.lpc.quantization_level, is_rice2? "RICE2":"RICE", subframe->data.lpc.entropy_coding_method.data.partitioned_rice.order);
 103                                 for(i = 0; i < subframe->data.lpc.order; i++)
 104                                         fprintf(fout, "\t\tqlp_coeff[%u]=%d\n", i, subframe->data.lpc.qlp_coeff[i]);
 105                                 for(i = 0; i < subframe->data.lpc.order; i++)
 106                                         fprintf(fout, "\t\twarmup[%u]=%d\n", i, subframe->data.lpc.warmup[i]);
 107                                 partitions = (1u << subframe->data.lpc.entropy_coding_method.data.partitioned_rice.order);
 108                                 for(i = 0; i < partitions; i++) {
 109                                         unsigned parameter = subframe->data.lpc.entropy_coding_method.data.partitioned_rice.contents->parameters[i];
 110                                         if(parameter == pesc)
 111                                                 fprintf(fout, "\t\tparameter[%u]=ESCAPE, raw_bits=%u\n", i, subframe->data.lpc.entropy_coding_method.data.partitioned_rice.contents->raw_bits[i]);
 112                                         else
 113                                                 fprintf(fout, "\t\tparameter[%u]=%u\n", i, parameter);
 114                                 }
 115                                 if(aopts.do_residual_text) {
 116                                         for(i = 0; i < frame->header.blocksize-subframe->data.lpc.order; i++)
 117                                                 fprintf(fout, "\t\tresidual[%u]=%d\n", i, subframe->data.lpc.residual[i]);
 118                                 }
 119                                 break;
 120                         case FLAC__SUBFRAME_TYPE_VERBATIM:
 121                                 fprintf(fout, "\n");
 122                                 break;
 123                 }
 124         }
 125
 126         /* now do the residual distributions if requested */
 127         if(aopts.do_residual_gnuplot) {
 128                 for(channel = 0; channel < channels; channel++) {
 129                         const FLAC__Subframe *subframe = frame->subframes+channel;
 130                         unsigned residual_samples;
 131
 132                         init_stats(&stats);
 133
 134                         switch(subframe->type) {
 135                                 case FLAC__SUBFRAME_TYPE_FIXED:
 136                                         residual_samples = frame->header.blocksize - subframe->data.fixed.order;
 137                                         for(i = 0; i < residual_samples; i++)
 138                                                 update_stats(&stats, subframe->data.fixed.residual[i], 1);
 139                                         break;
 140                                 case FLAC__SUBFRAME_TYPE_LPC:
 141                                         residual_samples = frame->header.blocksize - subframe->data.lpc.order;
 142                                         for(i = 0; i < residual_samples; i++)
 143                                                 update_stats(&stats, subframe->data.lpc.residual[i], 1);
 144                                         break;
 145                                 default:
 146                                         break;
 147                         }
 148
 149                         /* update all_ */
 150                         for(i = 0; i < stats.nbuckets; i++) {
 151                                 update_stats(&all_, stats.buckets[i].residual, stats.buckets[i].count);
 152                         }
 153
 154                         /* write the subframe */
 155                         flac_snprintf(outfilename, sizeof (outfilename), "f%06u.s%u.gp", frame_number, channel);
 156                         compute_stats(&stats);
 157
 158                         (void)dump_stats(&stats, outfilename);
 159                 }
 160         }
 161 }
 162
 163 void flac__analyze_finish(analysis_options aopts)
 164 {
 165         if(aopts.do_residual_gnuplot) {
 166                 compute_stats(&all_);
 167                 (void)dump_stats(&all_, "all");
 168         }
 169 }
 170
 171 void init_stats(subframe_stats_t *stats)
 172 {
 173         stats->peak_index = -1;
 174         stats->nbuckets = 0;
 175         stats->nsamples = 0;
 176         stats->sum = 0.0;
 177         stats->sos = 0.0;
 178 }
 179
 180 void update_stats(subframe_stats_t *stats, FLAC__int32 residual, unsigned incr)
 181 {
 182         unsigned i;
 183         const double r = (double)residual, a = r*incr;
 184
 185         stats->nsamples += incr;
 186         stats->sum += a;
 187         stats->sos += (a*r);
 188
 189         for(i = 0; i < stats->nbuckets; i++) {
 190                 if(stats->buckets[i].residual == residual) {
 191                         stats->buckets[i].count += incr;
 192                         goto find_peak;
 193                 }
 194         }
 195         /* not found, make a new bucket */
 196         i = stats->nbuckets;
 197         stats->buckets[i].residual = residual;
 198         stats->buckets[i].count = incr;
 199         stats->nbuckets++;
 200 find_peak:
 201         if(stats->peak_index < 0 || stats->buckets[i].count > stats->buckets[stats->peak_index].count)
 202                 stats->peak_index = i;
 203 }
 204
 205 void compute_stats(subframe_stats_t *stats)
 206 {
 207         stats->mean = stats->sum / (double)stats->nsamples;
 208         stats->variance = (stats->sos - (stats->sum * stats->sum / stats->nsamples)) / stats->nsamples;
 209         stats->stddev = sqrt(stats->variance);
 210 }
 211
 212 FLAC__bool dump_stats(const subframe_stats_t *stats, const char *filename)
 213 {
 214         FILE *outfile;
 215         unsigned i;
 216         const double m = stats->mean;
 217         const double s1 = stats->stddev, s2 = s1*2, s3 = s1*3, s4 = s1*4, s5 = s1*5, s6 = s1*6;
 218         const double p = stats->buckets[stats->peak_index].count;
 219
 220         outfile = fopen(filename, "w");
 221
 222         if(0 == outfile) {
 223                 fprintf(stderr, "ERROR opening %s: %s\n", filename, strerror(errno));
 224                 return false;
 225         }
 226
 227         fprintf(outfile, "plot '-' title 'PDF', '-' title 'mean' with impulses, '-' title '1-stddev' with histeps, '-' title '2-stddev' with histeps, '-' title '3-stddev' with histeps, '-' title '4-stddev' with histeps, '-' title '5-stddev' with histeps, '-' title '6-stddev' with histeps\n");
 228
 229         for(i = 0; i < stats->nbuckets; i++) {
 230                 fprintf(outfile, "%d %u\n", stats->buckets[i].residual, stats->buckets[i].count);
 231         }
 232         fprintf(outfile, "e\n");
 233
 234         fprintf(outfile, "%f %f\ne\n", stats->mean, p);
 235         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s1, p*0.8, m+s1, p*0.8);
 236         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s2, p*0.7, m+s2, p*0.7);
 237         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s3, p*0.6, m+s3, p*0.6);
 238         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s4, p*0.5, m+s4, p*0.5);
 239         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s5, p*0.4, m+s5, p*0.4);
 240         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s6, p*0.3, m+s6, p*0.3);
 241
 242         fprintf(outfile, "pause -1 'waiting...'\n");
 243
 244         fclose(outfile);
 245         return true;
 246 }