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