1 /* flac - Command-line FLAC encoder/decoder
2 * Copyright (C) 2000-2009 Josh Coalson
3 * Copyright (C) 2011-2023 Xiph.Org Foundation
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * as published by the Free Software Foundation; either version 2
8 * of the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License along
16 * with this program; if not, write to the Free Software Foundation, Inc.,
17 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19
20 #ifdef HAVE_CONFIG_H
21 # include <config.h>
22 #endif
23
24 #include <errno.h>
25 #include <math.h>
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29
30 #include "FLAC/all.h"
31 #include "analyze.h"
32
33 #include "share/compat.h"
34
35 typedef struct {
36 FLAC__int32 residual;
37 uint32_t count;
38 } pair_t;
39
40 typedef struct {
41 pair_t buckets[FLAC__MAX_BLOCK_SIZE];
42 int peak_index;
43 uint32_t nbuckets;
44 uint32_t nsamples;
45 double sum, sos;
46 double variance;
47 double mean;
48 double stddev;
49 } subframe_stats_t;
50
51 static subframe_stats_t all_;
52
53 static void init_stats(subframe_stats_t *stats);
54 static void update_stats(subframe_stats_t *stats, FLAC__int32 residual, uint32_t incr);
55 static void compute_stats(subframe_stats_t *stats);
56 static FLAC__bool dump_stats(const subframe_stats_t *stats, const char *filename);
57
flac__analyze_init(analysis_options aopts)58 void flac__analyze_init(analysis_options aopts)
59 {
60 if(aopts.do_residual_gnuplot) {
61 init_stats(&all_);
62 }
63 }
64
flac__analyze_frame(const FLAC__Frame * frame,uint32_t frame_number,FLAC__uint64 frame_offset,FLAC__uint64 frame_bytes,analysis_options aopts,FILE * fout)65 void flac__analyze_frame(const FLAC__Frame *frame, uint32_t frame_number, FLAC__uint64 frame_offset, FLAC__uint64 frame_bytes, analysis_options aopts, FILE *fout)
66 {
67 const uint32_t channels = frame->header.channels;
68 char outfilename[1024];
69 subframe_stats_t stats;
70 uint32_t i, channel, partitions;
71
72 /* do the human-readable part first */
73 fprintf(fout, "frame=%u\toffset=%" PRIu64 "\tbits=%" PRIu64 "\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]);
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 uint32_t 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=%" PRId64 "\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]=%" PRId64 "\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 uint32_t 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]=%" PRId64 "\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 uint32_t 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 uint32_t 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 if(stats.nsamples > 0) {
156 /* write the subframe */
157 flac_snprintf(outfilename, sizeof (outfilename), "f%06u.s%u.gp", frame_number, channel);
158 compute_stats(&stats);
159
160 (void)dump_stats(&stats, outfilename);
161 }
162 }
163 }
164 }
165
flac__analyze_finish(analysis_options aopts)166 void flac__analyze_finish(analysis_options aopts)
167 {
168 if(aopts.do_residual_gnuplot && all_.nsamples > 0) {
169 compute_stats(&all_);
170 (void)dump_stats(&all_, "all");
171 }
172 }
173
init_stats(subframe_stats_t * stats)174 void init_stats(subframe_stats_t *stats)
175 {
176 stats->peak_index = -1;
177 stats->nbuckets = 0;
178 stats->nsamples = 0;
179 stats->sum = 0.0;
180 stats->sos = 0.0;
181 }
182
update_stats(subframe_stats_t * stats,FLAC__int32 residual,uint32_t incr)183 void update_stats(subframe_stats_t *stats, FLAC__int32 residual, uint32_t incr)
184 {
185 uint32_t i;
186 const double r = (double)residual, a = r*incr;
187
188 stats->nsamples += incr;
189 stats->sum += a;
190 stats->sos += (a*r);
191
192 for(i = 0; i < stats->nbuckets; i++) {
193 if(stats->buckets[i].residual == residual) {
194 stats->buckets[i].count += incr;
195 goto find_peak;
196 }
197 }
198 /* not found, make a new bucket */
199 i = stats->nbuckets;
200 stats->buckets[i].residual = residual;
201 stats->buckets[i].count = incr;
202 stats->nbuckets++;
203 find_peak:
204 if(stats->peak_index < 0 || stats->buckets[i].count > stats->buckets[stats->peak_index].count)
205 stats->peak_index = i;
206 }
207
compute_stats(subframe_stats_t * stats)208 void compute_stats(subframe_stats_t *stats)
209 {
210 stats->mean = stats->sum / (double)stats->nsamples;
211 stats->variance = (stats->sos - (stats->sum * stats->sum / stats->nsamples)) / stats->nsamples;
212 stats->stddev = sqrt(stats->variance);
213 }
214
dump_stats(const subframe_stats_t * stats,const char * filename)215 FLAC__bool dump_stats(const subframe_stats_t *stats, const char *filename)
216 {
217 FILE *outfile;
218 uint32_t i;
219 const double m = stats->mean;
220 const double s1 = stats->stddev, s2 = s1*2, s3 = s1*3, s4 = s1*4, s5 = s1*5, s6 = s1*6;
221 const double p = stats->buckets[stats->peak_index].count;
222
223 outfile = flac_fopen(filename, "w");
224
225 if(0 == outfile) {
226 fprintf(stderr, "ERROR opening %s: %s\n", filename, strerror(errno));
227 return false;
228 }
229
230 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");
231
232 for(i = 0; i < stats->nbuckets; i++) {
233 fprintf(outfile, "%d %u\n", stats->buckets[i].residual, stats->buckets[i].count);
234 }
235 fprintf(outfile, "e\n");
236
237 fprintf(outfile, "%f %f\ne\n", stats->mean, p);
238 fprintf(outfile, "%f %f\n%f %f\ne\n", m-s1, p*0.8, m+s1, p*0.8);
239 fprintf(outfile, "%f %f\n%f %f\ne\n", m-s2, p*0.7, m+s2, p*0.7);
240 fprintf(outfile, "%f %f\n%f %f\ne\n", m-s3, p*0.6, m+s3, p*0.6);
241 fprintf(outfile, "%f %f\n%f %f\ne\n", m-s4, p*0.5, m+s4, p*0.5);
242 fprintf(outfile, "%f %f\n%f %f\ne\n", m-s5, p*0.4, m+s5, p*0.4);
243 fprintf(outfile, "%f %f\n%f %f\ne\n", m-s6, p*0.3, m+s6, p*0.3);
244
245 fprintf(outfile, "pause -1 'waiting...'\n");
246
247 fclose(outfile);
248 #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
249 unlink(filename);
250 #endif
251 return true;
252 }
253