From 3e3d9ec3af12e7c4cf492035bb6878f3b002ee6c Mon Sep 17 00:00:00 2001
From: Themaister <maister@archlinux.us>
Date: Thu, 14 Feb 2013 12:58:50 +0100
Subject: [PATCH] Use higher precision Kaiser.

Add test case for Kaiser window in MATLAB.
---
 audio/sinc.c                 |  2 +-
 audio/test/kaiser_window.m   |  6 ++++++
 audio/test/modified_bessel.m | 21 +++++++++++++++++++++
 audio/test/sinc_test.m       |  9 +++++++++
 4 files changed, 37 insertions(+), 1 deletion(-)
 create mode 100644 audio/test/kaiser_window.m
 create mode 100644 audio/test/modified_bessel.m

diff --git a/audio/sinc.c b/audio/sinc.c
index 1200c4b0a1..483e148a73 100644
--- a/audio/sinc.c
+++ b/audio/sinc.c
@@ -146,7 +146,7 @@ static inline double besseli0(double x)
 
    // Approximate. This is an infinite sum.
    // Luckily, it converges rather fast.
-   for (unsigned i = 0; i < 12; i++)
+   for (unsigned i = 0; i < 18; i++)
    {
       sum += x_pow * two_div_pow / (factorial * factorial);
 
diff --git a/audio/test/kaiser_window.m b/audio/test/kaiser_window.m
new file mode 100644
index 0000000000..a55517f08e
--- /dev/null
+++ b/audio/test/kaiser_window.m
@@ -0,0 +1,6 @@
+function win = kaiser_window(N, beta)
+% Create an N-point kaiser window with given beta.
+
+indices = 2 * (0 : N - 1) / (N - 1) - 1;
+mod = modified_bessel(beta);
+win = modified_bessel(beta * sqrt(1 - indices.^2)) / mod;
\ No newline at end of file
diff --git a/audio/test/modified_bessel.m b/audio/test/modified_bessel.m
new file mode 100644
index 0000000000..3ffa6d7c03
--- /dev/null
+++ b/audio/test/modified_bessel.m
@@ -0,0 +1,21 @@
+function val = modified_bessel(x)
+
+% Mirrors operation as done in RetroArch. Verify accuracy against Matlab's
+% implementation.
+
+sum = zeros(size(x));
+factorial = ones(size(x));
+factorial_mult = zeros(size(x));
+x_pow = ones(size(x));
+two_div_pow = ones(size(x));
+x_sqr = x .* x;
+
+for i = 0 : 17
+   sum = sum + x_pow .* two_div_pow ./ (factorial .* factorial);
+   factorial_mult = factorial_mult + 1.0;
+   x_pow = x_pow .* x_sqr;
+   two_div_pow = two_div_pow * 0.25;
+   factorial = factorial .* factorial_mult;
+end
+
+val = sum;
\ No newline at end of file
diff --git a/audio/test/sinc_test.m b/audio/test/sinc_test.m
index bbc0200696..57173ce833 100644
--- a/audio/test/sinc_test.m
+++ b/audio/test/sinc_test.m
@@ -1,6 +1,15 @@
 % MATLAB test case for RetroArch SINC upsampler.
 close all;
 
+%%
+% Test RetroArch's kaiser function.
+real_kaiser = kaiser(1024, 10.0)';
+rarch_kaiser = kaiser_window(1024, 10.0);
+figure('name', 'Bessel function test');
+subplot(2, 1, 1), plot(rarch_kaiser), title('RetroArch kaiser');
+subplot(2, 1, 2), plot(rarch_kaiser - real_kaiser), title('Error');
+
+%%
 % 4-tap and 8-tap are Lanczos windowed, but include here for completeness.
 phases = 256;
 bw = 0.375;