Clean up interpolation a bit.

diff --git a/Makefile b/Makefile
index 2f1b6728df609e8472b8940abd55e0ef7d7471eb..03d97ff26db0e74b02cfafdebc8d9998673efa9e 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-CXXFLAGS=-O2 -g -Wall
+CXXFLAGS=-O2 -ffast-math -g -Wall

 
 all: synth decode sync
 
 
 all: synth decode sync
 

diff --git a/decode.cpp b/decode.cpp
index 20c748ec6ab48bc68d3a81ee6fe32e21f9a48f07..ca0c675e411fbddfca3037cddbbd7fc1e69083e7 100644 (file)
--- a/decode.cpp
+++ b/decode.cpp
@@ -6,7 +6,8 @@
 #include <vector>
 #include <algorithm>
 
 #include <vector>
 #include <algorithm>
 

-#define LANCZOS_RADIUS 30
+#include "interpolate.h"
+

 #define BUFSIZE 4096
 #define HYSTERESIS_LIMIT 3000
 #define SAMPLE_RATE 44100
 #define BUFSIZE 4096
 #define HYSTERESIS_LIMIT 3000
 #define SAMPLE_RATE 44100


@@ -25,45 +26,6 @@ struct tap_header {
        unsigned int data_len;
 };
 
        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<short> &pcm, double i)
-{
-       int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
-       int upper = std::min<int>(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<short> &pcm, int x)
 {
 // between [x,x+1]
 double find_zerocrossing(const std::vector<short> &pcm, int x)
 {


@@ -81,7 +43,7 @@ double find_zerocrossing(const std::vector<short> &pcm, int x)
        double upper = x + 1;
        while (upper - lower > 1e-6) {
                double mid = 0.5f * (upper + lower);
        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;
                        upper = mid;
                } else {
                        lower = mid;

diff --git a/interpolate.h b/interpolate.h
new file mode 100644 (file)
index 0000000..2b4842a
--- /dev/null
+++ b/interpolate.h
@@ -0,0 +1,78 @@
+#ifndef _INTERPOLATE_H
+#define _INTERPOLATE_H 1
+
+#include <math.h>
+
+#include <algorithm>
+#include <vector>
+
+#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<class T>
+inline double lanczos_interpolate(const std::vector<T> &pcm, double i)
+{
+       int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
+       int upper = std::min<int>(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<class T>
+inline double lanczos_interpolate_right(const std::vector<T> &pcm, double i)
+{
+       int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
+       int upper = std::min<int>(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<class T>
+inline double linear_interpolate(const std::vector<T> &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<class T>
+inline double linear_interpolate_right(const std::vector<T> &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 4709b84271a2dde47ceed5326c2b037b6d7f82ad..b0fb6a6a0d68a04a8a34efca5ed2d646a5af6c5c 100644 (file)
--- a/sync.cpp
+++ b/sync.cpp
@@ -11,7 +11,8 @@
 #include <vector>
 #include <algorithm>
 
 #include <vector>
 #include <algorithm>
 

-#define LANCZOS_RADIUS 30
+#include "interpolate.h"
+

 #define BUFSIZE 4096
 
 struct stereo_sample {
 #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<double_stereo_sample> &pcm, double i)
-{
-       int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
-       int upper = std::min<int>(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<double> &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<class T>
-inline double interpolate(const std::vector<T> &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;
 struct hypothesis {
        int id;
        double cost;


@@ -171,11 +110,12 @@ int main(int argc, char **argv)
                bases[h].prev = NULL;
        }
 
                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) {
        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<size_t>(i + BUFSIZE, total_end);
        
                hypothesis *hyp = new hypothesis[NUM_THEORIES];
                size_t end = std::min<size_t>(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 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;
                        }
                                double diff = (right - left) * (right - left);
                                sum += diff;
                        }


@@ -209,6 +149,7 @@ int main(int argc, char **argv)
 
                prev_hyp = hyp;
        }
 
                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;
 
        // best winner
        double best_cost = HUGE_VAL;


@@ -236,14 +177,14 @@ int main(int argc, char **argv)
        }
        fclose(fp);
 
        }
        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<stereo_sample> aligned_pcm;
        std::vector<short> mono_pcm;
        for (unsigned i = 0; i < total_end; ++i) {
        double inv_sd = sqrt(2.0) / sqrt(var_left + var_right);
        std::vector<stereo_sample> aligned_pcm;
        std::vector<short> 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 left = pcm[i].left;

-               int right = interpolate(pcm, i + d);
+               int right = linear_interpolate_right(pcm, i + d);

 
                stereo_sample ss;
                ss.left = left;
 
                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))));
                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");
 
        fprintf(stderr, "Writing realigned stereo channels to aligned.raw...\n");
        fp = fopen("aligned.raw", "wb");