src/flac/analyze.c

   1 /* flac - Command-line FLAC encoder/decoder
   2  * Copyright (C) 2000,2001  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
  15  * along with this program; if not, write to the Free Software
  16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  17  */
  18
  19 #include <math.h>
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include "FLAC/all.h"
  24 #include "analyze.h"
  25
  26 typedef struct {
  27         FLAC__int32 residual;
  28         unsigned count;
  29 } pair_t;
  30
  31 typedef struct {
  32         pair_t buckets[FLAC__MAX_BLOCK_SIZE];
  33         int peak_index;
  34         unsigned nbuckets;
  35         unsigned nsamples;
  36         double sum, sos;
  37         double variance;
  38         double mean;
  39         double stddev;
  40 } subframe_stats_t;
  41
  42 static subframe_stats_t all_;
  43
  44 static void init_stats(subframe_stats_t *stats);
  45 static void update_stats(subframe_stats_t *stats, FLAC__int32 residual, unsigned incr);
  46 static void compute_stats(subframe_stats_t *stats);
  47 static FLAC__bool dump_stats(const subframe_stats_t *stats, const char *filename);
  48
  49 void flac__analyze_init(analysis_options aopts)
  50 {
  51         if(aopts.do_residual_gnuplot) {
  52                 init_stats(&all_);
  53         }
  54 }
  55
  56 void flac__analyze_frame(const FLAC__Frame *frame, unsigned frame_number, analysis_options aopts, FILE *fout)
  57 {
  58         const unsigned channels = frame->header.channels;
  59         char outfilename[1024];
  60         subframe_stats_t stats;
  61         unsigned i, channel;
  62
  63         /* do the human-readable part first */
  64         fprintf(fout, "frame=%u\tblocksize=%u\tsample_rate=%u\tchannels=%u\tchannel_assignment=%s\n", frame_number, frame->header.blocksize, frame->header.sample_rate, channels, FLAC__ChannelAssignmentString[frame->header.channel_assignment]);
  65         for(channel = 0; channel < channels; channel++) {
  66                 const FLAC__Subframe *subframe = frame->subframes+channel;
  67                 fprintf(fout, "\tsubframe=%u\twasted_bits=%u\ttype=%s", channel, subframe->wasted_bits, FLAC__SubframeTypeString[subframe->type]);
  68                 switch(subframe->type) {
  69                         case FLAC__SUBFRAME_TYPE_CONSTANT:
  70                                 fprintf(fout, "\tvalue=%d\n", subframe->data.constant.value);
  71                                 break;
  72                         case FLAC__SUBFRAME_TYPE_FIXED:
  73                                 fprintf(fout, "\torder=%u\tpartition_order=%u\n", subframe->data.fixed.order, subframe->data.fixed.entropy_coding_method.data.partitioned_rice.order); /*@@@ assumes method is partitioned-rice */
  74                                 for(i = 0; i < subframe->data.fixed.order; i++)
  75                                         fprintf(fout, "\t\twarmup[%u]=%d\n", i, subframe->data.fixed.warmup[i]);
  76                                 if(aopts.do_residual_text) {
  77                                         const unsigned partitions = (1u << subframe->data.fixed.entropy_coding_method.data.partitioned_rice.order);
  78                                         for(i = 0; i < partitions; i++) {
  79                                                 unsigned parameter = subframe->data.fixed.entropy_coding_method.data.partitioned_rice.parameters[i];
  80                                                 if(parameter == FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE_ESCAPE_PARAMETER)
  81                                                         fprintf(fout, "\t\tparameter[%u]=ESCAPE, raw_bits=%u\n", i, subframe->data.fixed.entropy_coding_method.data.partitioned_rice.raw_bits[i]);
  82                                                 else
  83                                                         fprintf(fout, "\t\tparameter[%u]=%u\n", i, parameter);
  84                                         }
  85                                         for(i = 0; i < frame->header.blocksize-subframe->data.fixed.order; i++)
  86                                                 fprintf(fout, "\t\tresidual[%u]=%d\n", i, subframe->data.fixed.residual[i]);
  87                                 }
  88                                 break;
  89                         case FLAC__SUBFRAME_TYPE_LPC:
  90                                 fprintf(fout, "\torder=%u\tpartition_order=%u\tqlp_coeff_precision=%u\tquantization_level=%d\n", subframe->data.lpc.order, subframe->data.lpc.entropy_coding_method.data.partitioned_rice.order, subframe->data.lpc.qlp_coeff_precision, subframe->data.lpc.quantization_level); /*@@@ assumes method is partitioned-rice */
  91                                 for(i = 0; i < subframe->data.lpc.order; i++)
  92                                         fprintf(fout, "\t\twarmup[%u]=%d\n", i, subframe->data.lpc.warmup[i]);
  93                                 if(aopts.do_residual_text) {
  94                                         const unsigned partitions = (1u << subframe->data.lpc.entropy_coding_method.data.partitioned_rice.order);
  95                                         for(i = 0; i < partitions; i++) {
  96                                                 unsigned parameter = subframe->data.lpc.entropy_coding_method.data.partitioned_rice.parameters[i];
  97                                                 if(parameter == FLAC__ENTROPY_CODING_METHOD_PARTITIONED_RICE_ESCAPE_PARAMETER)
  98                                                         fprintf(fout, "\t\tparameter[%u]=ESCAPE, raw_bits=%u\n", i, subframe->data.lpc.entropy_coding_method.data.partitioned_rice.raw_bits[i]);
  99                                                 else
 100                                                         fprintf(fout, "\t\tparameter[%u]=%u\n", i, parameter);
 101                                         }
 102                                         for(i = 0; i < frame->header.blocksize-subframe->data.lpc.order; i++)
 103                                                 fprintf(fout, "\t\tresidual[%u]=%d\n", i, subframe->data.lpc.residual[i]);
 104                                 }
 105                                 break;
 106                         case FLAC__SUBFRAME_TYPE_VERBATIM:
 107                                 fprintf(fout, "\n");
 108                                 break;
 109                 }
 110         }
 111
 112         /* now do the residual distributions if requested */
 113         if(aopts.do_residual_gnuplot) {
 114                 for(channel = 0; channel < channels; channel++) {
 115                         const FLAC__Subframe *subframe = frame->subframes+channel;
 116                         unsigned residual_samples;
 117
 118                         init_stats(&stats);
 119
 120                         switch(subframe->type) {
 121                                 case FLAC__SUBFRAME_TYPE_FIXED:
 122                                         residual_samples = frame->header.blocksize - subframe->data.fixed.order;
 123                                         for(i = 0; i < residual_samples; i++)
 124                                                 update_stats(&stats, subframe->data.fixed.residual[i], 1);
 125                                         break;
 126                                 case FLAC__SUBFRAME_TYPE_LPC:
 127                                         residual_samples = frame->header.blocksize - subframe->data.lpc.order;
 128                                         for(i = 0; i < residual_samples; i++)
 129                                                 update_stats(&stats, subframe->data.lpc.residual[i], 1);
 130                                         break;
 131                                 default:
 132                                         break;
 133                         }
 134
 135                         /* update all_ */
 136                         for(i = 0; i < stats.nbuckets; i++) {
 137                                 update_stats(&all_, stats.buckets[i].residual, stats.buckets[i].count);
 138                         }
 139
 140                         /* write the subframe */
 141                         sprintf(outfilename, "f%06u.s%u.gp", frame_number, channel);
 142                         compute_stats(&stats);
 143
 144                         (void)dump_stats(&stats, outfilename);
 145                 }
 146         }
 147 }
 148
 149 void flac__analyze_finish(analysis_options aopts)
 150 {
 151         if(aopts.do_residual_gnuplot) {
 152                 compute_stats(&all_);
 153                 (void)dump_stats(&all_, "all");
 154         }
 155 }
 156
 157 void init_stats(subframe_stats_t *stats)
 158 {
 159         stats->peak_index = -1;
 160         stats->nbuckets = 0;
 161         stats->nsamples = 0;
 162         stats->sum = 0.0;
 163         stats->sos = 0.0;
 164 }
 165
 166 void update_stats(subframe_stats_t *stats, FLAC__int32 residual, unsigned incr)
 167 {
 168         unsigned i;
 169         const double r = (double)residual, a = r*incr;
 170
 171         stats->nsamples += incr;
 172         stats->sum += a;
 173         stats->sos += (a*r);
 174
 175         for(i = 0; i < stats->nbuckets; i++) {
 176                 if(stats->buckets[i].residual == residual) {
 177                         stats->buckets[i].count += incr;
 178                         goto find_peak;
 179                 }
 180         }
 181         /* not found, make a new bucket */
 182         i = stats->nbuckets;
 183         stats->buckets[i].residual = residual;
 184         stats->buckets[i].count = incr;
 185         stats->nbuckets++;
 186 find_peak:
 187         if(stats->peak_index < 0 || stats->buckets[i].count > stats->buckets[stats->peak_index].count)
 188                 stats->peak_index = i;
 189 }
 190
 191 void compute_stats(subframe_stats_t *stats)
 192 {
 193         stats->mean = stats->sum / (double)stats->nsamples;
 194         stats->variance = (stats->sos - (stats->sum * stats->sum / stats->nsamples)) / stats->nsamples;
 195         stats->stddev = sqrt(stats->variance);
 196 }
 197
 198 FLAC__bool dump_stats(const subframe_stats_t *stats, const char *filename)
 199 {
 200         FILE *outfile;
 201         unsigned i;
 202         const double m = stats->mean;
 203         const double s1 = stats->stddev, s2 = s1*2, s3 = s1*3, s4 = s1*4, s5 = s1*5, s6 = s1*6;
 204         const double p = stats->buckets[stats->peak_index].count;
 205
 206         outfile = fopen(filename, "w");
 207
 208         if(0 == outfile) {
 209                 fprintf(stderr, "ERROR opening %s\n", filename);
 210                 return false;
 211         }
 212
 213         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");
 214
 215         for(i = 0; i < stats->nbuckets; i++) {
 216                 fprintf(outfile, "%d %u\n", stats->buckets[i].residual, stats->buckets[i].count);
 217         }
 218         fprintf(outfile, "e\n");
 219
 220         fprintf(outfile, "%f %f\ne\n", stats->mean, p);
 221         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s1, p*0.8, m+s1, p*0.8);
 222         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s2, p*0.7, m+s2, p*0.7);
 223         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s3, p*0.6, m+s3, p*0.6);
 224         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s4, p*0.5, m+s4, p*0.5);
 225         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s5, p*0.4, m+s5, p*0.4);
 226         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s6, p*0.3, m+s6, p*0.3);
 227
 228         fprintf(outfile, "pause -1 'waiting...'\n");
 229
 230         fclose(outfile);
 231         return true;
 232 }