add more verbose status reporting; pass aopts to analyze_init/finish()
[flac.git] / 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 <assert.h>
20 #include <math.h>
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include "FLAC/all.h"
25 #include "analyze.h"
26
27 typedef struct {
28         int32 residual;
29         unsigned count;
30 } pair_t;
31
32 typedef struct {
33         pair_t buckets[FLAC__MAX_BLOCK_SIZE];
34         int peak_index;
35         unsigned nbuckets;
36         unsigned nsamples;
37         double sum, sos;
38         double variance;
39         double mean;
40         double stddev;
41 } subframe_stats_t;
42
43 static subframe_stats_t all_;
44
45 static void init_stats(subframe_stats_t *stats);
46 static void update_stats(subframe_stats_t *stats, int32 residual, unsigned incr);
47 static void compute_stats(subframe_stats_t *stats);
48 static bool dump_stats(const subframe_stats_t *stats, const char *filename);
49
50 void analyze_init()
51 {
52         init_stats(&all_);
53 }
54
55 void analyze_frame(const FLAC__Frame *frame, unsigned frame_number, analysis_options aopts, FILE *fout)
56 {
57         const unsigned channels = frame->header.channels;
58         char outfilename[1024];
59         subframe_stats_t stats;
60         unsigned i, channel;
61
62         /* do the human-readable part first */
63         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]);
64         for(channel = 0; channel < channels; channel++) {
65                 const FLAC__Subframe *subframe = frame->subframes+channel;
66                 fprintf(fout, "\tsubframe=%u\twasted_bits=%u\ttype=%s", channel, subframe->wasted_bits, FLAC__SubframeTypeString[subframe->type]);
67                 switch(subframe->type) {
68                         case FLAC__SUBFRAME_TYPE_CONSTANT:
69                                 fprintf(fout, "\tvalue=%d\n", subframe->data.constant.value);
70                                 break;
71                         case FLAC__SUBFRAME_TYPE_FIXED:
72                                 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 */
73                                 for(i = 0; i < subframe->data.fixed.order; i++)
74                                         fprintf(fout, "\t\twarmup[%u]=%d\n", i, subframe->data.fixed.warmup[i]);
75                                 if(aopts.do_residual_text) {
76                                         for(i = 0; i < frame->header.blocksize-subframe->data.fixed.order; i++)
77                                                 fprintf(fout, "\t\tresidual[%u]=%d\n", i, subframe->data.fixed.residual[i]);
78                                 }
79                                 break;
80                         case FLAC__SUBFRAME_TYPE_LPC:
81                                 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 */
82                                 for(i = 0; i < subframe->data.lpc.order; i++)
83                                         fprintf(fout, "\t\twarmup[%u]=%d\n", i, subframe->data.lpc.warmup[i]);
84                                 if(aopts.do_residual_text) {
85                                         for(i = 0; i < frame->header.blocksize-subframe->data.lpc.order; i++)
86                                                 fprintf(fout, "\t\tresidual[%u]=%d\n", i, subframe->data.lpc.residual[i]);
87                                 }
88                                 break;
89                         case FLAC__SUBFRAME_TYPE_VERBATIM:
90                                 fprintf(fout, "\n");
91                                 break;
92                 }
93         }
94
95         /* now do the residual distributions if requested */
96         if(aopts.do_residual_gnuplot) {
97                 for(channel = 0; channel < channels; channel++) {
98                         const FLAC__Subframe *subframe = frame->subframes+channel;
99                         unsigned residual_samples;
100
101                         init_stats(&stats);
102
103                         switch(subframe->type) {
104                                 case FLAC__SUBFRAME_TYPE_FIXED:
105                                         residual_samples = frame->header.blocksize - subframe->data.fixed.order;
106                                         for(i = 0; i < residual_samples; i++)
107                                                 update_stats(&stats, subframe->data.fixed.residual[i], 1);
108                                         break;
109                                 case FLAC__SUBFRAME_TYPE_LPC:
110                                         residual_samples = frame->header.blocksize - subframe->data.lpc.order;
111                                         for(i = 0; i < residual_samples; i++)
112                                                 update_stats(&stats, subframe->data.lpc.residual[i], 1);
113                                         break;
114                                 default:
115                                         break;
116                         }
117
118                         /* update all_ */
119                         for(i = 0; i < stats.nbuckets; i++) {
120                                 update_stats(&all_, stats.buckets[i].residual, stats.buckets[i].count);
121                         }
122
123                         /* write the subframe */
124                         sprintf(outfilename, "f%06u.s%u.gp", frame_number, channel);
125                         compute_stats(&stats);
126 if(frame_number<50)//@@@
127                         (void)dump_stats(&stats, outfilename);
128                 }
129         }
130 }
131
132 void analyze_finish()
133 {
134         compute_stats(&all_);
135         (void)dump_stats(&all_, "all");
136 }
137
138 void init_stats(subframe_stats_t *stats)
139 {
140         stats->peak_index = -1;
141         stats->nbuckets = 0;
142         stats->nsamples = 0;
143         stats->sum = 0.0;
144         stats->sos = 0.0;
145 }
146
147 void update_stats(subframe_stats_t *stats, int32 residual, unsigned incr)
148 {
149         unsigned i;
150         const double r = (double)residual, a = r*incr;
151
152         stats->nsamples += incr;
153         stats->sum += a;
154         stats->sos += (a*r);
155
156         for(i = 0; i < stats->nbuckets; i++) {
157                 if(stats->buckets[i].residual == residual) {
158                         stats->buckets[i].count += incr;
159                         goto find_peak;
160                 }
161         }
162         /* not found, make a new bucket */
163         i = stats->nbuckets;
164         stats->buckets[i].residual = residual;
165         stats->buckets[i].count = incr;
166         stats->nbuckets++;
167 find_peak:
168         if(stats->peak_index < 0 || stats->buckets[i].count > stats->buckets[stats->peak_index].count)
169                 stats->peak_index = i;
170 }
171
172 void compute_stats(subframe_stats_t *stats)
173 {
174         stats->mean = stats->sum / (double)stats->nsamples;
175         stats->variance = (stats->sos - (stats->sum * stats->sum / stats->nsamples)) / stats->nsamples;
176         stats->stddev = sqrt(stats->variance);
177 }
178
179 bool dump_stats(const subframe_stats_t *stats, const char *filename)
180 {
181         FILE *outfile;
182         unsigned i;
183         const double m = stats->mean;
184         const double s1 = stats->stddev, s2 = s1*2, s3 = s1*3, s4 = s1*4, s5 = s1*5, s6 = s1*6;
185         const double p = stats->buckets[stats->peak_index].count;
186
187         outfile = fopen(filename, "w");
188
189         if(0 == outfile) {
190                 fprintf(stderr, "ERROR opening %s\n", filename);
191                 return false;
192         }
193
194         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");
195
196         for(i = 0; i < stats->nbuckets; i++) {
197                 fprintf(outfile, "%d %u\n", stats->buckets[i].residual, stats->buckets[i].count);
198         }
199         fprintf(outfile, "e\n");
200
201         fprintf(outfile, "%f %f\ne\n", stats->mean, p);
202         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s1, p*0.8, m+s1, p*0.8);
203         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s2, p*0.7, m+s2, p*0.7);
204         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s3, p*0.6, m+s3, p*0.6);
205         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s4, p*0.5, m+s4, p*0.5);
206         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s5, p*0.4, m+s5, p*0.4);
207         fprintf(outfile, "%f %f\n%f %f\ne\n", m-s6, p*0.3, m+s6, p*0.3);
208
209         fprintf(outfile, "pause -1 'waiting...'\n");
210
211         fclose(outfile);
212         return true;
213 }