+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' },
+ {"filter", required_argument, 0, 'f' },
+ {"rc-filter", required_argument, 0, 'r' },
+ {"output-filtered", 0, 0, 'F' },
+ {"crop", required_argument, 0, 'c' },
+ {"quiet", 0, 0, 'q' },
+ {"help", 0, 0, 'h' },
+ {0, 0, 0, 0 }
+};
+
+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..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 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");
+ 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);
+}
+
+void parse_options(int argc, char **argv)
+{
+ for ( ;; ) {
+ int option_index = 0;
+ int c = getopt_long(argc, argv, "ab:Am:spl:f:r:Fc:t:qh", long_options, &option_index);
+ if (c == -1)
+ break;
+
+ switch (c) {
+ case 'a':
+ 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);
+ break;
+
+ case 's':
+ do_calibrate = false;
+ break;
+
+ case 'p':
+ output_cycles_plot = true;
+ break;
+
+ 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, ": ");
+ int coeff_index = 0;
+ while (coeff_index < NUM_FILTER_COEFF && coeffstr != NULL) {
+ filter_coeff[coeff_index++] = atof(coeffstr);
+ coeffstr = strtok(NULL, ": ");
+ }
+ use_fir_filter = true;
+ break;
+ }
+
+ case 'r':
+ use_rc_filter = true;
+ rc_filter_freq = atof(optarg);
+ break;
+
+ case 'F':
+ 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;
+
+ case 'h':
+ default:
+ help();
+ exit(1);
+ }
+ }
+}
+
+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);
+}
+
+std::vector<float> do_fir_filter(const std::vector<float>& pcm, const float* filter)
+{
+ std::vector<float> filtered_pcm;
+ 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[i] = s;
+ }
+
+ if (output_filtered) {
+ FILE *fp = fopen("filtered.raw", "wb");
+ fwrite(filtered_pcm.data(), filtered_pcm.size() * sizeof(filtered_pcm[0]), 1, fp);
+ fclose(fp);
+ }
+
+ return filtered_pcm;
+}
+
+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());
+ 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;