Merge branch 'master' of /srv/git.sesse.net/www/c64tapwav
authorSteinar H. Gunderson <sgunderson@bigfoot.com>
Thu, 17 Mar 2016 23:36:28 +0000 (00:36 +0100)
committerSteinar H. Gunderson <sgunderson@bigfoot.com>
Thu, 17 Mar 2016 23:36:28 +0000 (00:36 +0100)
Makefile
decode.cpp
level.cpp
level.h

index 28ab7bc..2756d4f 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-CXXFLAGS=--std=gnu++0x -O2 -fno-math-errno -g -Wall
+CXXFLAGS=--std=gnu++0x -O2 -fno-math-errno -march=native -g -Wall
 LDLIBS=-lavcodec -lavformat -lavutil -lswresample
 
 all: synth decode sync cleaner
index d1b291f..d2d6fc8 100644 (file)
@@ -7,6 +7,9 @@
 #include <assert.h>
 #include <limits.h>
 #include <getopt.h>
+#ifdef __AVX__
+#include <immintrin.h>
+#endif
 #include <vector>
 #include <algorithm>
 
@@ -25,6 +28,7 @@
 
 // SPSA options
 #define NUM_FILTER_COEFF 32
+#define NUM_SPSA_VALS (NUM_FILTER_COEFF + 2)
 #define NUM_ITER 5000
 #define A NUM_ITER/10  // approx
 #define INITIAL_A 0.005 // A bit of trial and error...
@@ -32,7 +36,8 @@
 #define GAMMA 0.166
 #define ALPHA 1.0
 
-static float hysteresis_limit = 3000.0 / 32768.0;
+static float hysteresis_upper_limit = 0.1;
+static float hysteresis_lower_limit = -0.1;
 static bool do_calibrate = true;
 static bool output_cycles_plot = false;
 static bool do_crop = false;
@@ -50,37 +55,39 @@ static bool output_leveled = false;
 static std::vector<float> train_snap_points;
 static bool do_train = false;
 
+// The frequency to filter on (for do_auto_level), in Hertz.
+// Larger values makes the compressor react faster, but if it is too large,
+// you'll ruin the waveforms themselves.
+static float auto_level_freq = 200.0;
+
 // 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)
+// search for the value <limit> between [x,x+1]
+template<bool fast>
+double find_crossing(const std::vector<float> &pcm, int x, float limit)
 {
-       if (pcm[x] == 0) {
-               return x;
-       }
-       if (pcm[x + 1] == 0) {
-               return x + 1;
-       }
-
-       assert(pcm[x + 1] < 0);
-       assert(pcm[x] > 0);
-
-       double upper = x;
-       double lower = x + 1;
-       while (lower - upper > 1e-3) {
-               double mid = 0.5f * (upper + lower);
-               if (lanczos_interpolate(pcm, mid) > 0) {
-                       upper = mid;
-               } else {
-                       lower = mid;
+       if (fast) {
+               // Do simple linear interpolation.
+               return x + (limit - pcm[x]) / (pcm[x + 1] - pcm[x]);
+       } else {
+               // Binary search for the zero crossing as given by Lanczos interpolation.
+               double upper = x;
+               double lower = x + 1;
+               while (lower - upper > 1e-3) {
+                       double mid = 0.5f * (upper + lower);
+                       if (lanczos_interpolate(pcm, mid) > limit) {
+                               upper = mid;
+                       } else {
+                               lower = mid;
+                       }
                }
-       }
 
-       return 0.5f * (upper + lower);
+               return 0.5f * (upper + lower);
+       }
 }
 
 struct pulse {
@@ -178,7 +185,9 @@ void output_tap(const std::vector<pulse>& pulses, double calibration_factor)
 
 static struct option long_options[] = {
        {"auto-level",       0,                 0, 'a' },
+       {"auto-level-freq",  required_argument, 0, 'b' },
        {"output-leveled",   0,                 0, 'A' },
+       {"min-level",        required_argument, 0, 'm' },
        {"no-calibrate",     0,                 0, 's' },
        {"plot-cycles",      0,                 0, 'p' },
        {"hysteresis-limit", required_argument, 0, 'l' },
@@ -186,6 +195,7 @@ static struct option long_options[] = {
        {"rc-filter",        required_argument, 0, 'r' },
        {"output-filtered",  0,                 0, 'F' },
        {"crop",             required_argument, 0, 'c' },
+       {"train",            required_argument, 0, 't' },
        {"quiet",            0,                 0, 'q' },
        {"help",             0,                 0, 'h' },
        {0,                  0,                 0, 0   }
@@ -196,11 +206,12 @@ 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, "  -b, --auto-level-freq        minimum frequency in Hertz of corrected level changes (default 200 Hz)\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, "  -m, --min-level              minimum estimated sound level (0..1) 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, "  -l, --hysteresis-limit U[:L] change amplitude threshold for ignoring pulses (-1..1)\n");
        fprintf(stderr, "  -f, --filter C1:C2:C3:...    specify FIR filter (up to %d coefficients)\n", NUM_FILTER_COEFF);
        fprintf(stderr, "  -r, --rc-filter FREQ         send signal through a highpass RC filter with given frequency (in Hertz)\n");
        fprintf(stderr, "  -F, --output-filtered        output filtered waveform to filtered.raw\n");
@@ -216,7 +227,7 @@ void parse_options(int argc, char **argv)
 {
        for ( ;; ) {
                int option_index = 0;
-               int c = getopt_long(argc, argv, "aAm:spl:f:r:Fc:t:qh", long_options, &option_index);
+               int c = getopt_long(argc, argv, "ab:Am:spl:f:r:Fc:t:qh", long_options, &option_index);
                if (c == -1)
                        break;
 
@@ -225,12 +236,16 @@ void parse_options(int argc, char **argv)
                        do_auto_level = true;
                        break;
 
+               case 'b':
+                       auto_level_freq = atof(optarg);
+                       break;
+
                case 'A':
                        output_leveled = true;
                        break;
 
                case 'm':
-                       min_level = atof(optarg) / 32768.0;
+                       min_level = atof(optarg);
                        break;
 
                case 's':
@@ -241,9 +256,17 @@ void parse_options(int argc, char **argv)
                        output_cycles_plot = true;
                        break;
 
-               case 'l':
-                       hysteresis_limit = atof(optarg) / 32768.0;
+               case 'l': {
+                       const char *hyststr = strtok(optarg, ": ");
+                       hysteresis_upper_limit = atof(hyststr);
+                       hyststr = strtok(NULL, ": ");
+                       if (hyststr == NULL) {
+                               hysteresis_lower_limit = -hysteresis_upper_limit;
+                       } else {
+                               hysteresis_lower_limit = atof(hyststr);
+                       }
                        break;
+               }
 
                case 'f': {
                        const char *coeffstr = strtok(optarg, ": ");
@@ -316,17 +339,30 @@ std::vector<float> crop(const std::vector<float>& pcm, float crop_start, float c
        return std::vector<float>(pcm.begin() + start_sample, pcm.begin() + end_sample);
 }
 
-// TODO: Support AVX here.
 std::vector<float> do_fir_filter(const std::vector<float>& pcm, const float* filter)
 {
        std::vector<float> filtered_pcm;
-       filtered_pcm.reserve(pcm.size());
-       for (unsigned i = NUM_FILTER_COEFF; i < pcm.size(); ++i) {
+       filtered_pcm.resize(pcm.size());
+       unsigned i = NUM_FILTER_COEFF;
+#ifdef __AVX__
+       unsigned avx_end = i + ((pcm.size() - i) & ~7);
+       for ( ; i < avx_end; i += 8) {
+               __m256 s = _mm256_setzero_ps();
+               for (int j = 0; j < NUM_FILTER_COEFF; ++j) {
+                       __m256 f = _mm256_set1_ps(filter[j]);
+                       s = _mm256_fmadd_ps(f, _mm256_load_ps(&pcm[i - j]), s);
+               }
+               _mm256_storeu_ps(&filtered_pcm[i], s);
+       }
+#endif
+       // Do what we couldn't do with AVX (which is everything for non-AVX machines)
+       // as scalar code.
+       for (; i < pcm.size(); ++i) {
                float s = 0.0f;
                for (int j = 0; j < NUM_FILTER_COEFF; ++j) {
                        s += filter[j] * pcm[i - j];
                }
-               filtered_pcm.push_back(s);
+               filtered_pcm[i] = s;
        }
 
        if (output_filtered) {
@@ -340,11 +376,23 @@ std::vector<float> do_fir_filter(const std::vector<float>& pcm, const float* fil
 
 std::vector<float> do_rc_filter(const std::vector<float>& pcm, float freq, int sample_rate)
 {
+       // This is only a 6 dB/oct filter, which seemingly works better
+       // than the Filter class, which is a standard biquad (12 dB/oct).
+       // The b/c calculations come from libnyquist (atone.c);
+       // I haven't checked, but I suppose they fall out of the bilinear
+       // transform of the transfer function H(s) = s/(s + w).
        std::vector<float> filtered_pcm;
        filtered_pcm.resize(pcm.size());
-       Filter filter = Filter::hpf(M_PI * freq / sample_rate);
+       const float b = 2.0f - cos(2.0 * M_PI * freq / sample_rate);
+       const float c = b - sqrt(b * b - 1.0f);
+       float prev_in = 0.0f;
+       float prev_out = 0.0f;
        for (unsigned i = 0; i < pcm.size(); ++i) {
-               filtered_pcm[i] = filter.update(pcm[i]);
+               float in = pcm[i];
+               float out = c * (prev_out + in - prev_in);
+               filtered_pcm[i] = out;
+               prev_in = in;
+               prev_out = out;
        }
 
        if (output_filtered) {
@@ -356,49 +404,31 @@ std::vector<float> do_rc_filter(const std::vector<float>& pcm, float freq, int s
        return filtered_pcm;
 }
 
-std::vector<pulse> detect_pulses(const std::vector<float> &pcm, int sample_rate)
+template<bool fast>
+std::vector<pulse> detect_pulses(const std::vector<float> &pcm, float hysteresis_upper_limit, float hysteresis_lower_limit, int sample_rate)
 {
        std::vector<pulse> pulses;
 
        // Find the flanks.
-       int last_bit = -1;
+       enum State { START, ABOVE, BELOW } state = START;
        double last_downflank = -1;
        for (unsigned i = 0; i < pcm.size(); ++i) {
-               int bit = (pcm[i] > 0) ? 1 : 0;
-               if (bit == 0 && last_bit == 1) {
-                       // Check if we ever go up above <hysteresis_limit> before we dip down again.
-                       bool true_pulse = false;
-                       unsigned j;
-                       int min_level_after = 32767;
-                       for (j = i; j < pcm.size(); ++j) {
-                               min_level_after = std::min<int>(min_level_after, pcm[j]);
-                               if (pcm[j] > 0) break;
-                               if (pcm[j] < -hysteresis_limit) {
-                                       true_pulse = true;
-                                       break;
+               if (pcm[i] > hysteresis_upper_limit) {
+                       state = ABOVE;
+               } else if (pcm[i] < hysteresis_lower_limit) {
+                       if (state == ABOVE) {
+                               // down-flank!
+                               double t = find_crossing<fast>(pcm, i - 1, hysteresis_lower_limit) * (1.0 / sample_rate) + crop_start;
+                               if (last_downflank > 0) {
+                                       pulse p;
+                                       p.time = t;
+                                       p.len = t - last_downflank;
+                                       pulses.push_back(p);
                                }
+                               last_downflank = t;
                        }
-
-                       if (!true_pulse) {
-#if 0
-                               fprintf(stderr, "Ignored down-flank at %.6f seconds due to hysteresis (%d < %d).\n",
-                                       double(i) / sample_rate, -min_level_after, hysteresis_limit);
-#endif
-                               i = j;
-                               continue;
-                       } 
-
-                       // down-flank!
-                       double t = find_zerocrossing(pcm, i - 1) * (1.0 / sample_rate) + crop_start;
-                       if (last_downflank > 0) {
-                               pulse p;
-                               p.time = t;
-                               p.len = t - last_downflank;
-                               pulses.push_back(p);
-                       }
-                       last_downflank = t;
+                       state = BELOW;
                }
-               last_bit = bit;
        }
        return pulses;
 }
@@ -465,7 +495,7 @@ void find_kmeans(const std::vector<pulse> &pulses, double calibration_factor, co
                bool any_moved = false;
                for (unsigned i = 0; i < initial_centers.size(); ++i) {
                        if (num[i] == 0) {
-                               printf("K-means broke down, can't output new reference training points\n");
+                               fprintf(stderr, "K-means broke down, can't output new reference training points\n");
                                return;
                        }
                        float new_center = sums[i] / num[i];
@@ -478,16 +508,16 @@ void find_kmeans(const std::vector<pulse> &pulses, double calibration_factor, co
                        break;
                }
        }
-       printf("New reference training points:");
+       fprintf(stderr, "New reference training points:");
        for (unsigned i = 0; i < last_centers.size(); ++i) {
-               printf(" %.3f", last_centers[i]);
+               fprintf(stderr, " %.3f", last_centers[i]);
        }
-       printf("\n");
+       fprintf(stderr, "\n");
 }
 
 void spsa_train(const std::vector<float> &pcm, int sample_rate)
 {
-       float filter[NUM_FILTER_COEFF] = { 1.0f };  // The rest is filled with 0.
+       float vals[NUM_SPSA_VALS] = { hysteresis_upper_limit, hysteresis_lower_limit, 1.0f };  // The rest is filled with 0.
 
        float start_c = INITIAL_C;
        double best_badness = HUGE_VAL;
@@ -497,38 +527,38 @@ void spsa_train(const std::vector<float> &pcm, int sample_rate)
                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) {
+               float p[NUM_SPSA_VALS];
+               float vals1[NUM_SPSA_VALS], vals2[NUM_SPSA_VALS];
+               for (int i = 0; i < NUM_SPSA_VALS; ++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);
+                       vals1[i] = std::max(std::min(vals[i] - c * p[i], 1.0f), -1.0f);
+                       vals2[i] = std::max(std::min(vals[i] + c * p[i], 1.0f), -1.0f);
                }
 
-               std::vector<pulse> pulses1 = detect_pulses(do_fir_filter(pcm, filter1), sample_rate);
-               std::vector<pulse> pulses2 = detect_pulses(do_fir_filter(pcm, filter2), sample_rate);
+               std::vector<pulse> pulses1 = detect_pulses<true>(do_fir_filter(pcm, vals1 + 2), vals1[0], vals1[1], sample_rate);
+               std::vector<pulse> pulses2 = detect_pulses<true>(do_fir_filter(pcm, vals2 + 2), vals2[0], vals2[1], 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) {
+               float g[NUM_SPSA_VALS];
+               for (int i = 0; i < NUM_SPSA_VALS; ++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);
+                       vals[i] -= a * g[i];
+                       vals[i] = std::max(std::min(vals[i], 1.0f), -1.0f);
                }
                if (badness2 < badness1) {
                        std::swap(badness1, badness2);
-                       std::swap(filter1, filter2);
+                       std::swap(vals1, vals2);
                        std::swap(pulses1, pulses2);
                }
                if (badness1 < best_badness) {
-                       printf("\nNew best filter (badness=%f):", badness1);
+                       fprintf(stderr, "\nNew best filter (badness=%f):", badness1);
                        for (int i = 0; i < NUM_FILTER_COEFF; ++i) {
-                               printf(" %.5f", filter1[i]);
+                               fprintf(stderr, " %.5f", vals1[i + 2]);
                        }
+                       fprintf(stderr, ", hysteresis limits = %f %f\n", vals1[0], vals1[1]);
                        best_badness = badness1;
-                       printf("\n");
 
                        find_kmeans(pulses1, 1.0, train_snap_points);
 
@@ -536,8 +566,8 @@ void spsa_train(const std::vector<float> &pcm, int sample_rate)
                                output_cycle_plot(pulses1, 1.0);
                        }
                }
-               printf("%d ", n);
-               fflush(stdout);
+               fprintf(stderr, "%d ", n);
+               fflush(stderr);
        }
 }
 
@@ -565,7 +595,7 @@ int main(int argc, char **argv)
        }
 
        if (do_auto_level) {
-               pcm = level_samples(pcm, min_level, sample_rate);
+               pcm = level_samples(pcm, min_level, auto_level_freq, sample_rate);
                if (output_leveled) {
                        FILE *fp = fopen("leveled.raw", "wb");
                        fwrite(pcm.data(), pcm.size() * sizeof(pcm[0]), 1, fp);
@@ -590,7 +620,7 @@ int main(int argc, char **argv)
                exit(0);
        }
 
-       std::vector<pulse> pulses = detect_pulses(pcm, sample_rate);
+       std::vector<pulse> pulses = detect_pulses<false>(pcm, hysteresis_upper_limit, hysteresis_lower_limit, sample_rate);
 
        double calibration_factor = 1.0;
        if (do_calibrate) {
index 6e63d26..504cf0b 100644 (file)
--- a/level.cpp
+++ b/level.cpp
@@ -8,11 +8,6 @@
 
 #include "filter.h"
 
-// The frequency to filter on, in Hertz. Larger values makes the
-// compressor react faster, but if it is too large, you'll
-// ruin the waveforms themselves.
-#define LPFILTER_FREQ 50.0
-
 // A final scalar to get the audio within approximately the right range.
 // Increase to _lower_ overall volume.
 #define DAMPENING_FACTOR 5.0
@@ -20,7 +15,7 @@
 // 6dB/oct per round.
 #define FILTER_DEPTH 4
 
-std::vector<float> level_samples(const std::vector<float> &pcm, float min_level, int sample_rate)
+std::vector<float> level_samples(const std::vector<float> &pcm, float min_level, float lpfilter_freq, int sample_rate)
 {
        // filter forwards, then backwards (perfect phase filtering)
        std::vector<float> filtered_samples, refiltered_samples, leveled_samples;
@@ -28,7 +23,7 @@ std::vector<float> level_samples(const std::vector<float> &pcm, float min_level,
        refiltered_samples.resize(pcm.size());
        leveled_samples.resize(pcm.size());
 
-       Filter filter = Filter::lpf(M_PI * LPFILTER_FREQ / sample_rate);
+       Filter filter = Filter::lpf(2.0 * M_PI * lpfilter_freq / sample_rate);
        for (unsigned i = 0; i < pcm.size(); ++i) {
                filtered_samples[i] = filter.update(fabs(pcm[i]));
        }
diff --git a/level.h b/level.h
index 763e617..3a7f39d 100644 (file)
--- a/level.h
+++ b/level.h
@@ -3,6 +3,6 @@
 
 #include <vector>
 
-std::vector<float> level_samples(const std::vector<float> &pcm, float min_level, int sample_rate);
+std::vector<float> level_samples(const std::vector<float> &pcm, float min_level, float lpfilter_freq, int sample_rate);
 
 #endif  // !defined(_LEVEL_H)