1 /* Code taken from the OpenMPT project, which has a BSD
2 ** license which is compatible with this project.
3 **
4 ** The code has been slightly modified.
5 */
6
7 #include <stdint.h>
8 #include <stdbool.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include "ft2_windowed_sinc.h"
12
13 // globalized
14 float *fKaiserSinc = NULL;
15 float *fDownSample1 = NULL;
16 float *fDownSample2 = NULL;
17
Izero(double y)18 static double Izero(double y) // Compute Bessel function Izero(y) using a series approximation
19 {
20 double s = 1.0, ds = 1.0, d = 0.0;
21
22 const double epsilon = 1E-9; // 8bb: 1E-7 -> 1E-9 for added precision (still fast to calculate)
23
24 do
25 {
26 d = d + 2.0;
27 ds = ds * (y * y) / (d * d);
28 s = s + ds;
29 }
30 while (ds > epsilon * s);
31
32 return s;
33 }
34
getSinc(float * fLUTPtr,const double beta,const double cutoff)35 static void getSinc(float *fLUTPtr, const double beta, const double cutoff)
36 {
37 const double izeroBeta = Izero(beta);
38 const double kPi = 4.0 * atan(1.0) * cutoff; // M_PI can't be trusted
39
40 for (int32_t i = 0; i < SINC_LUT_LEN; i++)
41 {
42 double dSinc;
43 int32_t ix = (SINC_TAPS-1) - (i & (SINC_TAPS-1));
44
45 ix = (ix * SINC_PHASES) + (i >> SINC_TAPS_BITS);
46 if (ix == SINC_MID_TAP)
47 {
48 dSinc = 1.0;
49 }
50 else
51 {
52 const double x = (ix - SINC_MID_TAP) * (1.0 / SINC_PHASES);
53 const double xPi = x * kPi;
54
55 const double xMul = 1.0 / ((SINC_TAPS/2) * (SINC_TAPS/2));
56 dSinc = sin(xPi) * Izero(beta * sqrt(1.0 - x * x * xMul)) / (izeroBeta * xPi); // Kaiser window
57 }
58
59 fLUTPtr[i] = (float)(dSinc * cutoff);
60 }
61 }
62
calcWindowedSincTables(void)63 bool calcWindowedSincTables(void)
64 {
65 fKaiserSinc = (float *)malloc(SINC_LUT_LEN * sizeof (float));
66 fDownSample1 = (float *)malloc(SINC_LUT_LEN * sizeof (float));
67 fDownSample2 = (float *)malloc(SINC_LUT_LEN * sizeof (float));
68
69 if (fKaiserSinc == NULL || fDownSample1 == NULL || fDownSample2 == NULL)
70 return false;
71
72 getSinc(fKaiserSinc, 9.6377, 0.97);
73 getSinc(fDownSample1, 8.5, 0.5);
74 getSinc(fDownSample2, 2.7625, 0.425);
75
76 return true;
77 }
78
freeWindowedSincTables(void)79 void freeWindowedSincTables(void)
80 {
81 if (fKaiserSinc != NULL)
82 {
83 free(fKaiserSinc);
84 fKaiserSinc = NULL;
85 }
86
87 if (fDownSample1 != NULL)
88 {
89 free(fDownSample1);
90 fDownSample1 = NULL;
91 }
92
93 if (fDownSample2 != NULL)
94 {
95 free(fDownSample2);
96 fDownSample2 = NULL;
97 }
98 }
99