Skip to content

Commit 39aef0a

Browse files
8bitbubsy8bitbubsy
committed
sinc interpolation ringing fix
1 parent 3460fcf commit 39aef0a

File tree

2 files changed

+128
-139
lines changed

2 files changed

+128
-139
lines changed

‎src/ft2_header.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212
#endif
1313
#include "ft2_replayer.h"
1414

15-
#define PROG_VER_STR "2.02"
15+
#define PROG_VER_STR "2.03"
1616

1717
// do NOT change these! It will only mess things up...
1818

‎src/mixer/ft2_windowed_sinc.c‎

100644100755
Lines changed: 127 additions & 138 deletions
Original file line numberDiff line numberDiff line change
@@ -1,138 +1,127 @@
1-
// 8-point/16-point windowed-sinc interpolation LUT generator
2-
3-
#include <stdint.h>
4-
#include <stdbool.h>
5-
#include <stdlib.h>
6-
#include <math.h>
7-
#include "ft2_windowed_sinc.h"
8-
#include "../ft2_header.h" // PI
9-
#include "../ft2_video.h" // showErrorMsgBox()
10-
11-
typedef struct
12-
{
13-
double kaiserBeta, sincCutoff;
14-
} sincKernel_t;
15-
16-
// globalized
17-
float *fSinc[SINC_KERNELS], *fSinc8[SINC_KERNELS], *fSinc16[SINC_KERNELS];
18-
uint64_t sincRatio1, sincRatio2;
19-
20-
static sincKernel_t sincKernelConfig[2][SINC_KERNELS] =
21-
{
22-
/* Some notes on the Kaiser-Bessel beta parameter:
23-
** Lower beta = less treble cut off, more aliasing (narrower mainlobe, stronger sidelobe)
24-
** Higher beta = more treble cut off, less aliasing (wider mainlobe, weaker sidelobe)
25-
*/
26-
{ // -- settings for 8-point sinc --
27-
// beta, cutoff
28-
{ 9.20, 1.000 }, // kernel #1 (beta < ~9.2 leads to audible aliasing here)
29-
{ 8.50, 0.750 }, // kernel #2
30-
{ 7.30, 0.425 } // kernel #3
31-
},
32-
33-
{ // -- settings for 16-point sinc --
34-
// beta, cutoff
35-
{ 8.61, 1.000 }, // kernel #1 (beta 8.61 = Blackman-window approximation)
36-
{ 8.50, 0.750 }, // kernel #2
37-
{ 7.30, 0.425 } // kernel #3
38-
}
39-
};
40-
41-
// zeroth-order modified Bessel function of the first kind (series approximation)
42-
static inline double besselI0(double z)
43-
{
44-
double s = 1.0, ds = 1.0, d = 2.0;
45-
const double zz = z * z;
46-
47-
do
48-
{
49-
ds *= zz / (d * d);
50-
s += ds;
51-
d += 2.0;
52-
}
53-
while (ds > s*(1E-12));
54-
55-
return s;
56-
}
57-
58-
static inline double sinc(double x, double cutoff)
59-
{
60-
if (x == 0.0)
61-
{
62-
return cutoff;
63-
}
64-
else
65-
{
66-
x *= PI;
67-
return sin(cutoff * x) / x;
68-
}
69-
}
70-
71-
// note: numPoints/numPhases must be 2^n!
72-
static void makeSincKernel(float *out, int32_t numPoints, int32_t numPhases, double beta, double cutoff)
73-
{
74-
const int32_t kernelLen = numPhases * numPoints;
75-
const int32_t pointBits = (int32_t)log2(numPoints);
76-
const int32_t pointMask = numPoints - 1;
77-
const int32_t centerPoint = (numPoints / 2) - 1;
78-
const double besselI0Beta = 1.0 / besselI0(beta);
79-
const double phaseMul = 1.0 / numPhases;
80-
const double xMul = 1.0 / (numPoints / 2);
81-
82-
for (int32_t i = 0; i < kernelLen; i++)
83-
{
84-
const double x = ((i & pointMask) - centerPoint) - ((i >> pointBits) * phaseMul);
85-
86-
// Kaiser-Bessel window
87-
const double n = x * xMul;
88-
const double window = besselI0(beta * sqrt(1.0 - n * n)) * besselI0Beta;
89-
90-
out[i] = (float)(sinc(x, cutoff) * window);
91-
}
92-
}
93-
94-
bool setupWindowedSincTables(void)
95-
{
96-
sincKernel_t *k;
97-
for (int32_t i = 0; i < SINC_KERNELS; i++, k++)
98-
{
99-
fSinc8[i] = (float *)malloc( 8 * SINC_PHASES * sizeof (float));
100-
fSinc16[i] = (float *)malloc(16 * SINC_PHASES * sizeof (float));
101-
102-
if (fSinc8[i] == NULL || fSinc16[i] == NULL)
103-
{
104-
showErrorMsgBox("Not enough memory!");
105-
return false;
106-
}
107-
108-
k = &sincKernelConfig[0][i];
109-
makeSincKernel(fSinc8[i], 8, SINC_PHASES, k->kaiserBeta, k->sincCutoff);
110-
111-
k = &sincKernelConfig[1][i];
112-
makeSincKernel(fSinc16[i], 16, SINC_PHASES, k->kaiserBeta, k->sincCutoff);
113-
}
114-
115-
// resampling ratios for sinc kernel selection
116-
sincRatio1 = (uint64_t)(1.1875 * MIXER_FRAC_SCALE); // fSinc[0] if <=
117-
sincRatio2 = (uint64_t)(1.5000 * MIXER_FRAC_SCALE); // fSinc[1] if <=, else fSinc[2] if >
118-
119-
return true;
120-
}
121-
122-
void freeWindowedSincTables(void)
123-
{
124-
for (int32_t i = 0; i < SINC_KERNELS; i++)
125-
{
126-
if (fSinc8[i] != NULL)
127-
{
128-
free(fSinc8[i]);
129-
fSinc8[i] = NULL;
130-
}
131-
132-
if (fSinc16[i] != NULL)
133-
{
134-
free(fSinc16[i]);
135-
fSinc16[i] = NULL;
136-
}
137-
}
138-
}
1+
// 8-point/16-point windowed-sinc interpolation LUT generator
2+
3+
#include <stdint.h>
4+
#include <stdbool.h>
5+
#include <stdlib.h>
6+
#include <math.h>
7+
#include "ft2_windowed_sinc.h"
8+
#include "../ft2_header.h" // PI
9+
#include "../ft2_video.h" // showErrorMsgBox()
10+
11+
typedef struct
12+
{
13+
double kaiserBeta, sincCutoff;
14+
} sincKernel_t;
15+
16+
// globalized
17+
float *fSinc[SINC_KERNELS], *fSinc8[SINC_KERNELS], *fSinc16[SINC_KERNELS];
18+
uint64_t sincRatio1, sincRatio2;
19+
20+
static sincKernel_t sincKernelConfig[SINC_KERNELS] =
21+
{
22+
/* Some notes on the Kaiser-Bessel beta parameter:
23+
** Lower beta = less treble cut off, more aliasing (narrower mainlobe, stronger sidelobe)
24+
** Higher beta = more treble cut off, less aliasing (wider mainlobe, weaker sidelobe)
25+
*/
26+
27+
// beta, cutoff
28+
{ 9.6377, 1.000 }, // kernel #1 (lower beta results in audible ringing in some cases)
29+
{ 8.5000, 0.750 }, // kernel #2
30+
{ 7.3000, 0.425 } // kernel #3
31+
};
32+
33+
// zeroth-order modified Bessel function of the first kind (series approximation)
34+
static inline double besselI0(double z)
35+
{
36+
double s = 1.0, ds = 1.0, d = 2.0;
37+
const double zz = z * z;
38+
39+
do
40+
{
41+
ds *= zz / (d * d);
42+
s += ds;
43+
d += 2.0;
44+
}
45+
while (ds > s*(1E-12));
46+
47+
return s;
48+
}
49+
50+
static inline double sinc(double x, double cutoff)
51+
{
52+
if (x == 0.0)
53+
{
54+
return cutoff;
55+
}
56+
else
57+
{
58+
x *= PI;
59+
return sin(cutoff * x) / x;
60+
}
61+
}
62+
63+
// note: numPoints/numPhases must be 2^n!
64+
static void makeSincKernel(float *fOut, int32_t numPoints, int32_t numPhases, double beta, double cutoff)
65+
{
66+
const int32_t kernelLen = numPhases * numPoints;
67+
const int32_t pointBits = (int32_t)log2(numPoints);
68+
const int32_t pointMask = numPoints - 1;
69+
const int32_t centerPoint = (numPoints / 2) - 1;
70+
const double besselI0Beta = 1.0 / besselI0(beta);
71+
const double phaseMul = 1.0 / numPhases;
72+
const double xMul = 1.0 / (numPoints / 2);
73+
74+
for (int32_t i = 0; i < kernelLen; i++)
75+
{
76+
const double x = ((i & pointMask) - centerPoint) - ((i >> pointBits) * phaseMul);
77+
78+
// Kaiser-Bessel window
79+
const double n = x * xMul;
80+
const double window = besselI0(beta * sqrt(1.0 - n * n)) * besselI0Beta;
81+
82+
fOut[i] = (float)(sinc(x, cutoff) * window);
83+
}
84+
}
85+
86+
bool setupWindowedSincTables(void)
87+
{
88+
sincKernel_t *k = sincKernelConfig;
89+
for (int32_t i = 0; i < SINC_KERNELS; i++, k++)
90+
{
91+
fSinc8[i] = (float *)malloc( 8 * SINC_PHASES * sizeof (float));
92+
fSinc16[i] = (float *)malloc(16 * SINC_PHASES * sizeof (float));
93+
94+
if (fSinc8[i] == NULL || fSinc16[i] == NULL)
95+
{
96+
showErrorMsgBox("Not enough memory!");
97+
return false;
98+
}
99+
100+
makeSincKernel( fSinc8[i], 8, SINC_PHASES, k->kaiserBeta, k->sincCutoff);
101+
makeSincKernel(fSinc16[i], 16, SINC_PHASES, k->kaiserBeta, k->sincCutoff);
102+
}
103+
104+
// resampling ratios for sinc kernel selection
105+
sincRatio1 = (uint64_t)(1.1875 * MIXER_FRAC_SCALE); // fSinc[0] if <=
106+
sincRatio2 = (uint64_t)(1.5000 * MIXER_FRAC_SCALE); // fSinc[1] if <=, else fSinc[2] if >
107+
108+
return true;
109+
}
110+
111+
void freeWindowedSincTables(void)
112+
{
113+
for (int32_t i = 0; i < SINC_KERNELS; i++)
114+
{
115+
if (fSinc8[i] != NULL)
116+
{
117+
free(fSinc8[i]);
118+
fSinc8[i] = NULL;
119+
}
120+
121+
if (fSinc16[i] != NULL)
122+
{
123+
free(fSinc16[i]);
124+
fSinc16[i] = NULL;
125+
}
126+
}
127+
}

0 commit comments

Comments
 (0)