generate gnuplot files for subframe residuals

diff --git a/src/flac/analyze.c b/src/flac/analyze.c
index 5ae7856..f47d38c 100644 (file)
--- a/src/flac/analyze.c
+++ b/src/flac/analyze.c
@@ -17,6 +17,7 @@
  */
 
 #include <assert.h>
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -28,23 +29,35 @@ typedef struct {
        unsigned count;
 } pair_t;
 
-static pair_t super_buckets[FLAC__MAX_BLOCK_SIZE];
-static unsigned nsuper_buckets;
-
-static void update_buckets(int32 residual, pair_t *buckets, unsigned *nbuckets, unsigned incr);
-static bool print_buckets(const pair_t buckets[], const unsigned nbuckets, const char *filename);
+typedef struct {
+       pair_t buckets[FLAC__MAX_BLOCK_SIZE];
+       int peak_index;
+       unsigned nbuckets;
+       unsigned nsamples;
+       double sum, sos;
+       double variance;
+       double mean;
+       double stddev;
+} subframe_stats_t;
+
+static subframe_stats_t all_;
+
+static void init_stats(subframe_stats_t *stats);
+static void update_stats(subframe_stats_t *stats, int32 residual, unsigned incr);
+static void compute_stats(subframe_stats_t *stats);
+static bool dump_stats(const subframe_stats_t *stats, const char *filename);
 
 void analyze_init()
 {
-       nsuper_buckets = 0;
+       init_stats(&all_);
 }
 
 void analyze_frame(const FLAC__Frame *frame, unsigned frame_number, analysis_options aopts, FILE *fout)
 {
        const unsigned channels = frame->header.channels;
        char outfilename[1024];
-       pair_t buckets[FLAC__MAX_BLOCK_SIZE];
-       unsigned i, channel, nbuckets;
+       subframe_stats_t stats;
+       unsigned i, channel;
 
        /* do the human-readable part first */
        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]);
@@ -83,58 +96,93 @@ void analyze_frame(const FLAC__Frame *frame, unsigned frame_number, analysis_opt
        if(aopts.do_residual_gnuplot) {
                for(channel = 0; channel < channels; channel++) {
                        const FLAC__Subframe *subframe = frame->subframes+channel;
+                       unsigned residual_samples;
 
-                       nbuckets = 0;
+                       init_stats(&stats);
 
                        switch(subframe->type) {
                                case FLAC__SUBFRAME_TYPE_FIXED:
-                                       for(i = 0; i < frame->header.blocksize-subframe->data.fixed.order; i++)
-                                               update_buckets(subframe->data.fixed.residual[i], buckets, &nbuckets, 1);
+                                       residual_samples = frame->header.blocksize - subframe->data.fixed.order;
+                                       for(i = 0; i < residual_samples; i++)
+                                               update_stats(&stats, subframe->data.fixed.residual[i], 1);
                                        break;
                                case FLAC__SUBFRAME_TYPE_LPC:
-                                       for(i = 0; i < frame->header.blocksize-subframe->data.lpc.order; i++)
-                                               update_buckets(subframe->data.lpc.residual[i], buckets, &nbuckets, 1);
+                                       residual_samples = frame->header.blocksize - subframe->data.lpc.order;
+                                       for(i = 0; i < residual_samples; i++)
+                                               update_stats(&stats, subframe->data.lpc.residual[i], 1);
                                        break;
                                default:
                                        break;
                        }
 
-                       /* update super_buckets */
-                       for(i = 0; i < nbuckets; i++) {
-                               update_buckets(buckets[i].residual, super_buckets, &nsuper_buckets, buckets[i].count);
+                       /* update all_ */
+                       for(i = 0; i < stats.nbuckets; i++) {
+                               update_stats(&all_, stats.buckets[i].residual, stats.buckets[i].count);
                        }
 
                        /* write the subframe */
                        sprintf(outfilename, "f%06u.s%u.gp", frame_number, channel);
-                       (void)print_buckets(buckets, nbuckets, outfilename);
+                       compute_stats(&stats);
+if(frame_number<50)//@@@
+                       (void)dump_stats(&stats, outfilename);
                }
        }
 }
 
 void analyze_finish()
 {
-       (void)print_buckets(super_buckets, nsuper_buckets, "all");
+       compute_stats(&all_);
+       (void)dump_stats(&all_, "all");
+}
+
+void init_stats(subframe_stats_t *stats)
+{
+       stats->peak_index = -1;
+       stats->nbuckets = 0;
+       stats->nsamples = 0;
+       stats->sum = 0.0;
+       stats->sos = 0.0;
 }
 
-void update_buckets(int32 residual, pair_t *buckets, unsigned *nbuckets, unsigned incr)
+void update_stats(subframe_stats_t *stats, int32 residual, unsigned incr)
 {
        unsigned i;
+       const double r = (double)residual, a = r*incr;
+
+       stats->nsamples += incr;
+       stats->sum += a;
+       stats->sos += (a*r);
 
-       for(i = 0; i < *nbuckets; i++) {
-               if(buckets[i].residual == residual) {
-                       buckets[i].count += incr;
-                       return;
+       for(i = 0; i < stats->nbuckets; i++) {
+               if(stats->buckets[i].residual == residual) {
+                       stats->buckets[i].count += incr;
+                       goto find_peak;
                }
        }
-       buckets[*nbuckets].residual = residual;
-       buckets[*nbuckets].count = incr;
-       (*nbuckets)++;
+       /* not found, make a new bucket */
+       i = stats->nbuckets;
+       stats->buckets[i].residual = residual;
+       stats->buckets[i].count = incr;
+       stats->nbuckets++;
+find_peak:
+       if(stats->peak_index < 0 || stats->buckets[i].count > stats->buckets[stats->peak_index].count)
+               stats->peak_index = i;
+}
+
+void compute_stats(subframe_stats_t *stats)
+{
+       stats->mean = stats->sum / (double)stats->nsamples;
+       stats->variance = (stats->sos - (stats->sum * stats->sum / stats->nsamples)) / stats->nsamples;
+       stats->stddev = sqrt(stats->variance);
 }
 
-bool print_buckets(const pair_t buckets[], const unsigned nbuckets, const char *filename)
+bool dump_stats(const subframe_stats_t *stats, const char *filename)
 {
        FILE *outfile;
        unsigned i;
+       const double m = stats->mean;
+       const double s1 = stats->stddev, s2 = s1*2, s3 = s1*3, s4 = s1*4, s5 = s1*5, s6 = s1*6;
+       const double p = stats->buckets[stats->peak_index].count;
 
        outfile = fopen(filename, "w");
 
@@ -143,9 +191,22 @@ bool print_buckets(const pair_t buckets[], const unsigned nbuckets, const char *
                return false;
        }
 
-       for(i = 0; i < nbuckets; i++) {
-               fprintf(outfile, "%d %u\n", buckets[i].residual, buckets[i].count);
+       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");
+
+       for(i = 0; i < stats->nbuckets; i++) {
+               fprintf(outfile, "%d %u\n", stats->buckets[i].residual, stats->buckets[i].count);
        }
+       fprintf(outfile, "e\n");
+
+       fprintf(outfile, "%f %f\ne\n", stats->mean, p);
+       fprintf(outfile, "%f %f\n%f %f\ne\n", m-s1, p*0.8, m+s1, p*0.8);
+       fprintf(outfile, "%f %f\n%f %f\ne\n", m-s2, p*0.7, m+s2, p*0.7);
+       fprintf(outfile, "%f %f\n%f %f\ne\n", m-s3, p*0.6, m+s3, p*0.6);
+       fprintf(outfile, "%f %f\n%f %f\ne\n", m-s4, p*0.5, m+s4, p*0.5);
+       fprintf(outfile, "%f %f\n%f %f\ne\n", m-s5, p*0.4, m+s5, p*0.4);
+       fprintf(outfile, "%f %f\n%f %f\ne\n", m-s6, p*0.3, m+s6, p*0.3);
+
+       fprintf(outfile, "pause -1 'waiting...'\n");
 
        fclose(outfile);
        return true;