From 52dc441eb7db99f2231478474f228a62d642eacf Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sat, 4 May 2013 20:59:21 +0200 Subject: [PATCH] Clean up interpolation a bit. --- Makefile | 2 +- decode.cpp | 44 ++------------------------- interpolate.h | 78 ++++++++++++++++++++++++++++++++++++++++++++++++ sync.cpp | 82 +++++++++------------------------------------------ 4 files changed, 96 insertions(+), 110 deletions(-) create mode 100644 interpolate.h diff --git a/Makefile b/Makefile index 2f1b672..03d97ff 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -CXXFLAGS=-O2 -g -Wall +CXXFLAGS=-O2 -ffast-math -g -Wall all: synth decode sync diff --git a/decode.cpp b/decode.cpp index 20c748e..ca0c675 100644 --- a/decode.cpp +++ b/decode.cpp @@ -6,7 +6,8 @@ #include #include -#define LANCZOS_RADIUS 30 +#include "interpolate.h" + #define BUFSIZE 4096 #define HYSTERESIS_LIMIT 3000 #define SAMPLE_RATE 44100 @@ -25,45 +26,6 @@ struct tap_header { unsigned int data_len; }; -double sinc(double x) -{ - if (fabs(x) < 1e-6) { - return 1.0f - fabs(x); - } else { - return sin(x) / x; - } -} - -#if 0 -double weight(double x) -{ - if (fabs(x) > LANCZOS_RADIUS) { - return 0.0f; - } - return sinc(M_PI * x) * sinc(M_PI * x / LANCZOS_RADIUS); -} -#else -double weight(double x) -{ - if (fabs(x) > 1.0f) { - return 0.0f; - } - return 1.0f - fabs(x); -} -#endif - -double interpolate(const std::vector &pcm, double i) -{ - int lower = std::max(ceil(i - LANCZOS_RADIUS), 0); - int upper = std::min(floor(i + LANCZOS_RADIUS), pcm.size() - 1); - double sum = 0.0f; - - for (int x = lower; x <= upper; ++x) { - sum += pcm[x] * weight(i - x); - } - return sum; -} - // between [x,x+1] double find_zerocrossing(const std::vector &pcm, int x) { @@ -81,7 +43,7 @@ double find_zerocrossing(const std::vector &pcm, int x) double upper = x + 1; while (upper - lower > 1e-6) { double mid = 0.5f * (upper + lower); - if (interpolate(pcm, mid) > 0) { + if (lanczos_interpolate(pcm, mid) > 0) { upper = mid; } else { lower = mid; diff --git a/interpolate.h b/interpolate.h new file mode 100644 index 0000000..2b4842a --- /dev/null +++ b/interpolate.h @@ -0,0 +1,78 @@ +#ifndef _INTERPOLATE_H +#define _INTERPOLATE_H 1 + +#include + +#include +#include + +#define LANCZOS_RADIUS 30 + +inline double sinc(double x) +{ + if (fabs(x) < 1e-6) { + return 1.0f - fabs(x); + } else { + return sin(x) / x; + } +} + +inline double lanczos_weight(double x) +{ + if (fabs(x) > LANCZOS_RADIUS) { + return 0.0f; + } + return sinc(M_PI * x) * sinc(M_PI * x / LANCZOS_RADIUS); +} + +template +inline double lanczos_interpolate(const std::vector &pcm, double i) +{ + int lower = std::max(ceil(i - LANCZOS_RADIUS), 0); + int upper = std::min(floor(i + LANCZOS_RADIUS), pcm.size() - 1); + double sum = 0.0f; + + for (int x = lower; x <= upper; ++x) { + sum += pcm[x] * lanczos_weight(i - x); + } + return sum; +} + +template +inline double lanczos_interpolate_right(const std::vector &pcm, double i) +{ + int lower = std::max(ceil(i - LANCZOS_RADIUS), 0); + int upper = std::min(floor(i + LANCZOS_RADIUS), pcm.size() - 1); + double sum = 0.0f; + + for (int x = lower; x <= upper; ++x) { + sum += pcm[x].right * lanczos_weight(i - x); + } + return sum; +} + +template +inline double linear_interpolate(const std::vector &pcm, double i) +{ + int ii = int(i); + if (ii < 0 || ii >= int(pcm.size() - 1)) { + return 0.0; + } + double frac = i - ii; + + return pcm[ii] + frac * (pcm[ii + 1] - pcm[ii]); +} + +template +inline double linear_interpolate_right(const std::vector &pcm, double i) +{ + int ii = int(i); + if (ii < 0 || ii >= int(pcm.size() - 1)) { + return 0.0; + } + double frac = i - ii; + + return pcm[ii].right + frac * (pcm[ii + 1].right - pcm[ii].right); +} + +#endif // !defined(_INTERPOLATE_H) diff --git a/sync.cpp b/sync.cpp index 4709b84..b0fb6a6 100644 --- a/sync.cpp +++ b/sync.cpp @@ -11,7 +11,8 @@ #include #include -#define LANCZOS_RADIUS 30 +#include "interpolate.h" + #define BUFSIZE 4096 struct stereo_sample { @@ -32,68 +33,6 @@ inline short clip(int x) } } -double sinc(double x) -{ - if (fabs(x) < 1e-6) { - return 1.0f - fabs(x); - } else { - return sin(x) / x; - } -} - -#if 0 -double weight(double x) -{ - if (fabs(x) > LANCZOS_RADIUS) { - return 0.0f; - } - return sinc(M_PI * x) * sinc(M_PI * x / LANCZOS_RADIUS); -} - -double interpolate(const std::vector &pcm, double i) -{ - int lower = std::max(ceil(i - LANCZOS_RADIUS), 0); - int upper = std::min(floor(i + LANCZOS_RADIUS), pcm.size() - 1); - double sum = 0.0f; - - for (int x = lower; x <= upper; ++x) { - sum += pcm[x].right * weight(i - x); - } - return sum; -} -#else -double weight(double x) -{ - if (fabs(x) > 1.0f) { - return 0.0f; - } - return 1.0f - fabs(x); -} - -inline double interpolate(const std::vector &pcm, double i) -{ - int ii = int(i); - if (ii < 0 || ii >= int(pcm.size() - 1)) { - return 0.0; - } - double frac = i - ii; - - return pcm[ii] + frac * (pcm[ii + 1] - pcm[ii]); -} - -template -inline double interpolate(const std::vector &pcm, double i) -{ - int ii = int(i); - if (ii < 0 || ii >= int(pcm.size() - 1)) { - return 0.0; - } - double frac = i - ii; - - return pcm[ii].right + frac * (pcm[ii + 1].right - pcm[ii].right); -} -#endif - struct hypothesis { int id; double cost; @@ -171,11 +110,12 @@ int main(int argc, char **argv) bases[h].prev = NULL; } + fprintf(stderr, "Matching blocks... %7.2f", 0.0); hypothesis *prev_hyp = bases; size_t total_end = pcm.size(); //size_t total_end = 4410000; for (unsigned i = 0; i < total_end; i += BUFSIZE) { - fprintf(stderr, "%.3f\n", i / 44100.0); + fprintf(stderr, "\b\b\b\b\b\b\b%7.2f", i / 44100.0); size_t end = std::min(i + BUFSIZE, total_end); hypothesis *hyp = new hypothesis[NUM_THEORIES]; @@ -186,7 +126,7 @@ int main(int argc, char **argv) double d = delays[h]; for (unsigned s = i; s < end; ++s) { double left = norm[s].left; - double right = interpolate(norm, s + d); + double right = linear_interpolate_right(norm, s + d); double diff = (right - left) * (right - left); sum += diff; } @@ -209,6 +149,7 @@ int main(int argc, char **argv) prev_hyp = hyp; } + fprintf(stderr, "\b\b\b\b\b\b\b%7.2f", total_end / 44100.0); // best winner double best_cost = HUGE_VAL; @@ -236,14 +177,14 @@ int main(int argc, char **argv) } fclose(fp); - // readjust and write out + fprintf(stderr, "Stretching right channel to match left... %7.2f%%", 0.0); double inv_sd = sqrt(2.0) / sqrt(var_left + var_right); std::vector aligned_pcm; std::vector mono_pcm; for (unsigned i = 0; i < total_end; ++i) { - double d = interpolate(best_path, i / double(BUFSIZE)); + double d = linear_interpolate(best_path, i / double(BUFSIZE)); int left = pcm[i].left; - int right = interpolate(pcm, i + d); + int right = linear_interpolate_right(pcm, i + d); stereo_sample ss; ss.left = left; @@ -251,7 +192,12 @@ int main(int argc, char **argv) aligned_pcm.push_back(ss); mono_pcm.push_back(clip(lrintf(inv_sd * 4096.0 * (left + right)))); + + if (i % 4096 == 0) { + fprintf(stderr, "\b\b\b\b\b\b\b\b%7.2f%%", 100.0 * i / total_end); + } } + fprintf(stderr, "\b\b\b\b\b\b\b%7.2f%%\n", 100.0); fprintf(stderr, "Writing realigned stereo channels to aligned.raw...\n"); fp = fopen("aligned.raw", "wb"); -- 2.39.2