From 6a20d03691c7c34dcf4a126909586199ef26cf94 Mon Sep 17 00:00:00 2001 From: Josh Coalson Date: Tue, 20 Mar 2001 22:55:06 +0000 Subject: [PATCH] generate gnuplot files for subframe residuals --- src/flac/analyze.c | 119 ++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 90 insertions(+), 29 deletions(-) diff --git a/src/flac/analyze.c b/src/flac/analyze.c index 5ae7856..f47d38c 100644 --- a/src/flac/analyze.c +++ b/src/flac/analyze.c @@ -17,6 +17,7 @@ */ #include +#include #include #include #include @@ -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; -- 2.7.4