]> git.sesse.net Git - c64tapwav/blobdiff - decode.cpp
Make leveler minimum level configurable, since it is connected to hysteresis limit.
[c64tapwav] / decode.cpp
index 2a27371ac2db191c78c0119e812f8bad01834096..b3b8aa3376f355b9fcad9e8cd3aa4fc9b1bf2b4b 100644 (file)
@@ -9,6 +9,7 @@
 
 #include "audioreader.h"
 #include "interpolate.h"
+#include "level.h"
 #include "tap.h"
 
 #define BUFSIZE 4096
 #define SYNC_PULSE_LENGTH 378.0
 #define SYNC_TEST_TOLERANCE 1.10
 
+// SPSA options
 #define NUM_FILTER_COEFF 32
+#define NUM_ITER 5000
+#define A NUM_ITER/10  // approx
+#define INITIAL_A 0.005 // A bit of trial and error...
+#define INITIAL_C 0.02  // This too.
+#define GAMMA 0.166
+#define ALPHA 1.0
 
 static float hysteresis_limit = 3000.0 / 32768.0;
 static bool do_calibrate = true;
 static bool output_cycles_plot = false;
 static bool use_filter = false;
+static bool do_crop = false;
+static float crop_start = 0.0f, crop_end = HUGE_VAL;
 static float filter_coeff[NUM_FILTER_COEFF] = { 1.0f };  // The rest is filled with 0.
 static bool output_filtered = false;
 static bool quiet = false;
+static bool do_auto_level = false;
+static bool output_leveled = false;
+static std::vector<float> train_snap_points;
+static bool do_train = false;
+
+// The minimum estimated sound level (for do_auto_level) at any given point.
+// If you decrease this, you'll be able to amplify really silent signals
+// by more, but you'll also increase the level of silent (ie. noise-only) segments,
+// possibly caused misdetected pulses in these segments.
+static float min_level = 0.05f;
 
 // between [x,x+1]
 double find_zerocrossing(const std::vector<float> &pcm, int x)
@@ -149,11 +169,14 @@ void output_tap(const std::vector<pulse>& pulses, double calibration_factor)
 }
 
 static struct option long_options[] = {
+       {"auto-level",       0,                 0, 'a' },
+       {"output-leveled",   0,                 0, 'A' },
        {"no-calibrate",     0,                 0, 's' },
        {"plot-cycles",      0,                 0, 'p' },
        {"hysteresis-limit", required_argument, 0, 'l' },
        {"filter",           required_argument, 0, 'f' },
        {"output-filtered",  0,                 0, 'F' },
+       {"crop",             required_argument, 0, 'c' },
        {"quiet",            0,                 0, 'q' },
        {"help",             0,                 0, 'h' },
        {0,                  0,                 0, 0   }
@@ -163,11 +186,17 @@ void help()
 {
        fprintf(stderr, "decode [OPTIONS] AUDIO-FILE > TAP-FILE\n");
        fprintf(stderr, "\n");
+       fprintf(stderr, "  -a, --auto-level             automatically adjust amplitude levels throughout the file\n");
+       fprintf(stderr, "  -A, --output-leveled         output leveled waveform to leveled.raw\n");
+       fprintf(stderr, "  -m, --min-level              minimum estimated sound level (0..32768) for --auto-level\n");
        fprintf(stderr, "  -s, --no-calibrate           do not try to calibrate on sync pulse length\n");
        fprintf(stderr, "  -p, --plot-cycles            output debugging info to cycles.plot\n");
        fprintf(stderr, "  -l, --hysteresis-limit VAL   change amplitude threshold for ignoring pulses (0..32768)\n");
        fprintf(stderr, "  -f, --filter C1:C2:C3:...    specify FIR filter (up to %d coefficients)\n", NUM_FILTER_COEFF);
        fprintf(stderr, "  -F, --output-filtered        output filtered waveform to filtered.raw\n");
+       fprintf(stderr, "  -c, --crop START[:END]       use only the given part of the file\n");
+       fprintf(stderr, "  -t, --train LEN1:LEN2:...    train a filter for detecting any of the given number of cycles\n");
+       fprintf(stderr, "                               (implies --no-calibrate and --quiet unless overridden)\n");
        fprintf(stderr, "  -q, --quiet                  suppress some informational messages\n");
        fprintf(stderr, "  -h, --help                   display this help, then exit\n");
        exit(1);
@@ -177,11 +206,23 @@ void parse_options(int argc, char **argv)
 {
        for ( ;; ) {
                int option_index = 0;
-               int c = getopt_long(argc, argv, "spl:f:Fqh", long_options, &option_index);
+               int c = getopt_long(argc, argv, "aAm:spl:f:Fc:t:qh", long_options, &option_index);
                if (c == -1)
                        break;
 
                switch (c) {
+               case 'a':
+                       do_auto_level = true;
+                       break;
+
+               case 'A':
+                       output_leveled = true;
+                       break;
+
+               case 'm':
+                       min_level = atof(optarg) / 32768.0;
+                       break;
+
                case 's':
                        do_calibrate = false;
                        break;
@@ -195,11 +236,11 @@ void parse_options(int argc, char **argv)
                        break;
 
                case 'f': {
-                       const char *coeffstr = strtok(optarg, ":");
+                       const char *coeffstr = strtok(optarg, ": ");
                        int coeff_index = 0;
                        while (coeff_index < NUM_FILTER_COEFF && coeffstr != NULL) {
                                filter_coeff[coeff_index++] = atof(coeffstr);
-                               coeffstr = strtok(NULL, ":");
+                               coeffstr = strtok(NULL, ": ");
                        }
                        use_filter = true;
                        break;
@@ -209,6 +250,33 @@ void parse_options(int argc, char **argv)
                        output_filtered = true;
                        break;
 
+               case 'c': {
+                       const char *cropstr = strtok(optarg, ":");
+                       crop_start = atof(cropstr);
+                       cropstr = strtok(NULL, ":");
+                       if (cropstr == NULL) {
+                               crop_end = HUGE_VAL;
+                       } else {
+                               crop_end = atof(cropstr);
+                       }
+                       do_crop = true;
+                       break;
+               }
+
+               case 't': {
+                       const char *cyclestr = strtok(optarg, ":");
+                       while (cyclestr != NULL) {
+                               train_snap_points.push_back(atof(cyclestr));
+                               cyclestr = strtok(NULL, ":");
+                       }
+                       do_train = true;
+
+                       // Set reasonable defaults (can be overridden later on the command line).
+                       do_calibrate = false;
+                       quiet = true;
+                       break;
+               }
+
                case 'q':
                        quiet = true;
                        break;
@@ -221,6 +289,18 @@ void parse_options(int argc, char **argv)
        }
 }
 
+std::vector<float> crop(const std::vector<float>& pcm, float crop_start, float crop_end, int sample_rate)
+{
+       size_t start_sample, end_sample;
+       if (crop_start >= 0.0f) {
+               start_sample = std::min<size_t>(lrintf(crop_start * sample_rate), pcm.size());
+       }
+       if (crop_end >= 0.0f) {
+               end_sample = std::min<size_t>(lrintf(crop_end * sample_rate), pcm.size());
+       }
+       return std::vector<float>(pcm.begin() + start_sample, pcm.begin() + end_sample);
+}
+
 // TODO: Support AVX here.
 std::vector<float> do_filter(const std::vector<float>& pcm, const float* filter)
 {
@@ -276,7 +356,7 @@ std::vector<pulse> detect_pulses(const std::vector<float> &pcm, int sample_rate)
                        } 
 
                        // down-flank!
-                       double t = find_zerocrossing(pcm, i - 1) * (1.0 / sample_rate);
+                       double t = find_zerocrossing(pcm, i - 1) * (1.0 / sample_rate) + crop_start;
                        if (last_downflank > 0) {
                                pulse p;
                                p.time = t;
@@ -290,6 +370,86 @@ std::vector<pulse> detect_pulses(const std::vector<float> &pcm, int sample_rate)
        return pulses;
 }
 
+void output_cycle_plot(const std::vector<pulse> &pulses, double calibration_factor)
+{
+       FILE *fp = fopen("cycles.plot", "w");
+       for (unsigned i = 0; i < pulses.size(); ++i) {
+               double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
+               fprintf(fp, "%f %f\n", pulses[i].time, cycles);
+       }
+       fclose(fp);
+}
+
+float eval_badness(const std::vector<pulse>& pulses, double calibration_factor)
+{
+       double sum_badness = 0.0;
+       for (unsigned i = 0; i < pulses.size(); ++i) {
+               double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
+               if (cycles > 2000.0) cycles = 2000.0;  // Don't make pauses arbitrarily bad.
+               double badness = (cycles - train_snap_points[0]) * (cycles - train_snap_points[0]);
+               for (unsigned j = 1; j < train_snap_points.size(); ++j) {
+                       badness = std::min(badness, (cycles - train_snap_points[j]) * (cycles - train_snap_points[j]));
+               }
+               sum_badness += badness;
+       }
+       return sqrt(sum_badness / (pulses.size() - 1));
+}
+
+void spsa_train(std::vector<float> &pcm, int sample_rate)
+{
+       // Train!
+       float filter[NUM_FILTER_COEFF] = { 1.0f };  // The rest is filled with 0.
+
+       float start_c = INITIAL_C;
+       double best_badness = HUGE_VAL;
+
+       for (int n = 1; n < NUM_ITER; ++n) {
+               float a = INITIAL_A * pow(n + A, -ALPHA);
+               float c = start_c * pow(n, -GAMMA);
+
+               // find a random perturbation
+               float p[NUM_FILTER_COEFF];
+               float filter1[NUM_FILTER_COEFF], filter2[NUM_FILTER_COEFF];
+               for (int i = 0; i < NUM_FILTER_COEFF; ++i) {
+                       p[i] = (rand() % 2) ? 1.0 : -1.0;
+                       filter1[i] = std::max(std::min(filter[i] - c * p[i], 1.0f), -1.0f);
+                       filter2[i] = std::max(std::min(filter[i] + c * p[i], 1.0f), -1.0f);
+               }
+
+               std::vector<pulse> pulses1 = detect_pulses(do_filter(pcm, filter1), sample_rate);
+               std::vector<pulse> pulses2 = detect_pulses(do_filter(pcm, filter2), sample_rate);
+               float badness1 = eval_badness(pulses1, 1.0);
+               float badness2 = eval_badness(pulses2, 1.0);
+
+               // Find the gradient estimator
+               float g[NUM_FILTER_COEFF];
+               for (int i = 0; i < NUM_FILTER_COEFF; ++i) {
+                       g[i] = (badness2 - badness1) / (2.0 * c * p[i]);
+                       filter[i] -= a * g[i];
+                       filter[i] = std::max(std::min(filter[i], 1.0f), -1.0f);
+               }
+               if (badness2 < badness1) {
+                       std::swap(badness1, badness2);
+                       std::swap(filter1, filter2);
+                       std::swap(pulses1, pulses2);
+               }
+               if (badness1 < best_badness) {
+                       printf("\nNew best filter (badness=%f):", badness1);
+                       for (int i = 0; i < NUM_FILTER_COEFF; ++i) {
+                               printf(" %.5f", filter1[i]);
+                       }
+                       best_badness = badness1;
+                       printf("\n");
+
+                       if (output_cycles_plot) {
+                               output_cycle_plot(pulses1, 1.0);
+                       }
+               }
+               printf("%d ", n);
+               fflush(stdout);
+       }
+}
+
 int main(int argc, char **argv)
 {
        parse_options(argc, argv);
@@ -301,10 +461,23 @@ int main(int argc, char **argv)
                exit(1);
        }
 
+       if (do_crop) {
+               pcm = crop(pcm, crop_start, crop_end, sample_rate);
+       }
+
        if (use_filter) {
                pcm = do_filter(pcm, filter_coeff);
        }
 
+       if (do_auto_level) {
+               pcm = level_samples(pcm, min_level, sample_rate);
+               if (output_leveled) {
+                       FILE *fp = fopen("leveled.raw", "wb");
+                       fwrite(pcm.data(), pcm.size() * sizeof(pcm[0]), 1, fp);
+                       fclose(fp);
+               }
+       }
+
 #if 0
        for (int i = 0; i < LEN; ++i) {
                in[i] += rand() % 10000;
@@ -317,6 +490,11 @@ int main(int argc, char **argv)
        }
 #endif
 
+       if (do_train) {
+               spsa_train(pcm, sample_rate);
+               exit(0);
+       }
+
        std::vector<pulse> pulses = detect_pulses(pcm, sample_rate);
 
        double calibration_factor = 1.0;
@@ -325,12 +503,7 @@ int main(int argc, char **argv)
        }
 
        if (output_cycles_plot) {
-               FILE *fp = fopen("cycles.plot", "w");
-               for (unsigned i = 0; i < pulses.size(); ++i) {
-                       double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
-                       fprintf(fp, "%f %f\n", pulses[i].time, cycles);
-               }
-               fclose(fp);
+               output_cycle_plot(pulses, calibration_factor);
        }
 
        output_tap(pulses, calibration_factor);