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