restructure files into directories by year and "Misc" to ease testing with current...

[fw/tmflights] / Misc / kaiser-test.5c

diff --git a/Misc/kaiser-test.5c b/Misc/kaiser-test.5c
new file mode 100644 (file)
index 0000000..41fd5e5
--- /dev/null
+++ b/Misc/kaiser-test.5c
@@ -0,0 +1,58 @@
+load "kaiser.5c"
+
+real sinc(real x) = x != 0 ? sin(x)/x : 1;
+
+real sum(real[...] x) { real s = 0; for(int i = 0; i < dim(x); i++) s += x[i]; return s; }
+
+real[...] convolve(real[...] d, real[...] e) {
+       real sample(n) = n < 0 ? d[0] : n >= dim(d) ? d[dim(d)-1] : d[n];
+       real w = floor ((dim(e) - 1) / 2);
+       real c(int center) {
+               real    v = 0;
+               real    s = 0;
+               int     start;
+               int     end;
+
+               start = center - w;
+               if (start < 0)
+                       start = 0;
+               end = center + w;
+               if (end >= dim(d))
+                       end = dim(d) - 1;
+               for (int o = start; o <= end; o++) {
+                       real    e_o = e[o - center + w];
+                       real    part = d[o] * e_o;
+                       v += part;
+                       s += e_o;
+               }
+               v /= s;
+               return v;
+       }
+       return (real[dim(d)]) { [n] = c(n) };
+}
+
+real[...] kaiser_filter(real[...] d, int half_width) {
+#      real[half_width * 2 + 1] fir = { [n] = sinc(2 * pi * n / (2 * half_width)) };
+       real M = half_width * 2 + 1;
+       real[M] fir = { [n] = kaiser(n+1, M+1, 8) };
+       return convolve(d, fir);
+}
+
+
+real[1024] square = { [n] = 2 * ((n >> 8) & 1) - 1 };
+
+real[] cut = low_pass_filter(π*1/128, π*1/64, 1e-10);
+
+real[dim(square)] filtered_square = convolve(square, cut);
+
+for (int i = 0; i < dim(filtered_square); i++)
+    printf ("%d %g %g\n", i, square[i], filtered_square[i]);
+
+exit(0);
+
+int size = 100;
+real[size] raw = { [n] = n };
+real[size] filtered = kaiser_filter(raw, 8);
+
+for (int i = 0; i < size; i++)
+       printf ("%d %g %g\n", i, raw[i], filtered[i]);