generate gnuplot files for subframe residuals

This commit is contained in:
Josh Coalson 2001-03-20 22:55:06 +00:00
parent 73a9935105
commit 6a20d03691

View File

@ -17,6 +17,7 @@
*/ */
#include <assert.h> #include <assert.h>
#include <math.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
@ -28,23 +29,35 @@ typedef struct {
unsigned count; unsigned count;
} pair_t; } pair_t;
static pair_t super_buckets[FLAC__MAX_BLOCK_SIZE]; typedef struct {
static unsigned nsuper_buckets; 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 void update_buckets(int32 residual, pair_t *buckets, unsigned *nbuckets, unsigned incr); static subframe_stats_t all_;
static bool print_buckets(const pair_t buckets[], const unsigned nbuckets, const char *filename);
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() void analyze_init()
{ {
nsuper_buckets = 0; init_stats(&all_);
} }
void analyze_frame(const FLAC__Frame *frame, unsigned frame_number, analysis_options aopts, FILE *fout) void analyze_frame(const FLAC__Frame *frame, unsigned frame_number, analysis_options aopts, FILE *fout)
{ {
const unsigned channels = frame->header.channels; const unsigned channels = frame->header.channels;
char outfilename[1024]; char outfilename[1024];
pair_t buckets[FLAC__MAX_BLOCK_SIZE]; subframe_stats_t stats;
unsigned i, channel, nbuckets; unsigned i, channel;
/* do the human-readable part first */ /* 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]); 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) { if(aopts.do_residual_gnuplot) {
for(channel = 0; channel < channels; channel++) { for(channel = 0; channel < channels; channel++) {
const FLAC__Subframe *subframe = frame->subframes+channel; const FLAC__Subframe *subframe = frame->subframes+channel;
unsigned residual_samples;
nbuckets = 0; init_stats(&stats);
switch(subframe->type) { switch(subframe->type) {
case FLAC__SUBFRAME_TYPE_FIXED: case FLAC__SUBFRAME_TYPE_FIXED:
for(i = 0; i < frame->header.blocksize-subframe->data.fixed.order; i++) residual_samples = frame->header.blocksize - subframe->data.fixed.order;
update_buckets(subframe->data.fixed.residual[i], buckets, &nbuckets, 1); for(i = 0; i < residual_samples; i++)
update_stats(&stats, subframe->data.fixed.residual[i], 1);
break; break;
case FLAC__SUBFRAME_TYPE_LPC: case FLAC__SUBFRAME_TYPE_LPC:
for(i = 0; i < frame->header.blocksize-subframe->data.lpc.order; i++) residual_samples = frame->header.blocksize - subframe->data.lpc.order;
update_buckets(subframe->data.lpc.residual[i], buckets, &nbuckets, 1); for(i = 0; i < residual_samples; i++)
update_stats(&stats, subframe->data.lpc.residual[i], 1);
break; break;
default: default:
break; break;
} }
/* update super_buckets */ /* update all_ */
for(i = 0; i < nbuckets; i++) { for(i = 0; i < stats.nbuckets; i++) {
update_buckets(buckets[i].residual, super_buckets, &nsuper_buckets, buckets[i].count); update_stats(&all_, stats.buckets[i].residual, stats.buckets[i].count);
} }
/* write the subframe */ /* write the subframe */
sprintf(outfilename, "f%06u.s%u.gp", frame_number, channel); 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 analyze_finish()
{ {
(void)print_buckets(super_buckets, nsuper_buckets, "all"); compute_stats(&all_);
(void)dump_stats(&all_, "all");
} }
void update_buckets(int32 residual, pair_t *buckets, unsigned *nbuckets, unsigned incr) 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_stats(subframe_stats_t *stats, int32 residual, unsigned incr)
{ {
unsigned i; unsigned i;
const double r = (double)residual, a = r*incr;
for(i = 0; i < *nbuckets; i++) { stats->nsamples += incr;
if(buckets[i].residual == residual) { stats->sum += a;
buckets[i].count += incr; stats->sos += (a*r);
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; /* not found, make a new bucket */
buckets[*nbuckets].count = incr; i = stats->nbuckets;
(*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;
} }
bool print_buckets(const pair_t buckets[], const unsigned nbuckets, const char *filename) 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 dump_stats(const subframe_stats_t *stats, const char *filename)
{ {
FILE *outfile; FILE *outfile;
unsigned i; 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"); outfile = fopen(filename, "w");
@ -143,9 +191,22 @@ bool print_buckets(const pair_t buckets[], const unsigned nbuckets, const char *
return false; return false;
} }
for(i = 0; i < nbuckets; i++) { 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");
fprintf(outfile, "%d %u\n", buckets[i].residual, buckets[i].count);
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); fclose(outfile);
return true; return true;