diff --git a/audio/sinc.c b/audio/sinc.c index fb912a3371..1d68cea627 100644 --- a/audio/sinc.c +++ b/audio/sinc.c @@ -14,7 +14,7 @@ */ // Bog-standard windowed SINC implementation. -// Only suitable as an upsampler, as there is no low-pass filter stage. +// Only suitable as an upsampler, as cutoff frequency isn't dynamically configurable (yet). #include "resampler.h" #include "../performance.h" @@ -33,16 +33,45 @@ #include #endif -#ifdef SINC_LOWER_QUALITY +// Rough SNR values for upsampling: +// LOWEST: 40 dB +// LOWER: 55 dB +// NORMAL: 70 dB +// HIGHER: 110 dB +// HIGHEST: 140 dB + +// TODO, make all this more configurable. +#if defined(SINC_LOWEST_QUALITY) +#define SINC_WINDOW_LANCZOS +#define CUTOFF 0.98 +#define PHASE_BITS 11 +#define SIDELOBES 2 +#define ENABLE_AVX 0 +#elif defined(SINC_LOWER_QUALITY) +#define SINC_WINDOW_LANCZOS +#define CUTOFF 0.98 #define PHASE_BITS 12 #define SIDELOBES 4 #define ENABLE_AVX 0 #elif defined(SINC_HIGHER_QUALITY) +#define SINC_WINDOW_KAISER +#define SINC_WINDOW_KAISER_BETA 10.5 +#define CUTOFF 0.90 #define PHASE_BITS 16 #define SIDELOBES 32 #define ENABLE_AVX 1 -#else +#elif defined(SINC_HIGHEST_QUALITY) +#define SINC_WINDOW_KAISER +#define SINC_WINDOW_KAISER_BETA 14.5 +#define CUTOFF 0.95 #define PHASE_BITS 16 +#define SIDELOBES 128 +#define ENABLE_AVX 1 +#else +#define SINC_WINDOW_KAISER +#define SINC_WINDOW_KAISER_BETA 5.5 +#define CUTOFF 0.825 +#define PHASE_BITS 14 #define SIDELOBES 8 #define ENABLE_AVX 0 #endif @@ -60,7 +89,6 @@ #define PHASES (1 << (PHASE_BITS + SUBPHASE_BITS)) #define TAPS (SIDELOBES * 2) -#define CUTOFF 0.98 typedef struct rarch_sinc_resampler { @@ -80,23 +108,62 @@ static inline double sinc(double val) return sin(val) / val; } -static inline double lanzcos(double index) +#if defined(SINC_WINDOW_LANCZOS) +static inline double window_function(double index) { - return sinc(index); + return sinc(M_PI * index); +} +#elif defined(SINC_WINDOW_KAISER) +// Modified Bessel function of first order. +// Check Wiki for mathematical definition ... +static inline double besseli0(double x) +{ + double sum = 0.0; + + double factorial = 1.0; + double factorial_mult = 0.0; + double x_pow = 1.0; + double two_div_pow = 1.0; + + // Approximate. This is an infinite sum. + // Luckily, it converges rather fast. + for (unsigned i = 0; i < 12; i++) + { + sum += x_pow * two_div_pow / (factorial * factorial); + + factorial_mult += 1.0; + x_pow *= x * x; + two_div_pow *= 0.25; + factorial *= factorial_mult; + } + + return sum; } -static void init_sinc_table(rarch_sinc_resampler_t *resamp) +static inline double window_function(double index) { - // Sinc phases: [..., p + 3, p + 2, p + 1, p + 0, p - 1, p - 2, p - 3, p - 4, ...] - for (int i = 0; i < (1 << PHASE_BITS); i++) - { - for (int j = 0; j < TAPS; j++) - { - double p = (double)i / (1 << PHASE_BITS); - double sinc_phase = M_PI * (p + (SIDELOBES - 1 - j)); + return besseli0(SINC_WINDOW_KAISER_BETA * sqrt(1 - index * index)); +} +#else +#error "No SINC window function defined." +#endif - float val = CUTOFF * sinc(CUTOFF * sinc_phase) * lanzcos(sinc_phase / SIDELOBES); - resamp->phase_table[i][j] = val; +static void init_sinc_table(rarch_sinc_resampler_t *resamp, double cutoff, + float *phase_table, int phases, int taps) +{ + double window_mod = window_function(0.0); // Need to normalize w(0) to 1.0. + + for (int i = 0; i < phases; i++) + { + for (int j = 0; j < taps; j++) + { + int n = j * phases + i; + double window_phase = (double)n / (((1 << PHASE_BITS) * TAPS) - 1.0); // [0, 1]. + window_phase = 2.0 * window_phase - 1.0; // [-1, 1] + double sinc_phase = SIDELOBES * window_phase; + + float val = cutoff * sinc(M_PI * sinc_phase * cutoff) * window_function(window_phase) / window_mod; + phase_table[i * taps + j] = val; } } } @@ -263,9 +330,9 @@ static void resampler_sinc_process(void *re_, struct resampler_data *data) { while (frames && re->time >= PHASES) { + re->ptr = (re->ptr - 1) & (TAPS - 1); re->buffer_l[re->ptr + TAPS] = re->buffer_l[re->ptr] = *input++; re->buffer_r[re->ptr + TAPS] = re->buffer_r[re->ptr] = *input++; - re->ptr = (re->ptr + 1) & (TAPS - 1); re->time -= PHASES; frames--; @@ -296,7 +363,7 @@ static void *resampler_sinc_new(void) memset(re, 0, sizeof(*re)); - init_sinc_table(re); + init_sinc_table(re, CUTOFF, &re->phase_table[0][0], 1 << PHASE_BITS, TAPS); #if defined(__AVX__) && ENABLE_AVX RARCH_LOG("Sinc resampler [AVX]\n"); diff --git a/audio/test/Makefile b/audio/test/Makefile index 55ea3498ea..1e149ffacc 100644 --- a/audio/test/Makefile +++ b/audio/test/Makefile @@ -1,4 +1,15 @@ -TESTS := test-hermite test-sinc test-snr-sinc test-snr-hermite +TESTS := test-hermite \ + test-snr-hermite \ + test-sinc-lowest \ + test-snr-sinc-lowest \ + test-sinc-lower \ + test-snr-sinc-lower \ + test-sinc \ + test-snr-sinc \ + test-sinc-higher \ + test-snr-sinc-higher \ + test-sinc-highest \ + test-snr-sinc-highest CFLAGS += -O3 -g -Wall -pedantic -std=gnu99 -DRESAMPLER_TEST LDFLAGS += -lm @@ -8,12 +19,6 @@ all: $(TESTS) test-hermite: ../hermite.o ../utils.o main.o resampler-hermite.o $(CC) -o $@ $^ $(LDFLAGS) -test-sinc: ../sinc.o ../utils.o main.o ../hermite.o resampler-sinc.o - $(CC) -o $@ $^ $(LDFLAGS) - -test-snr-sinc: ../sinc.o ../utils.o snr.o ../hermite.o resampler-sinc.o - $(CC) -o $@ $^ $(LDFLAGS) - test-snr-hermite: ../hermite.o ../utils.o snr.o resampler-hermite.o $(CC) -o $@ $^ $(LDFLAGS) @@ -23,6 +28,51 @@ resampler-sinc.o: ../resampler.c resampler-hermite.o: ../resampler.c $(CC) -c -o $@ $< $(CFLAGS) +sinc-lowest.o: ../sinc.c + $(CC) -c -o $@ $< $(CFLAGS) -DHAVE_SINC -DSINC_LOWEST_QUALITY + +sinc-lower.o: ../sinc.c + $(CC) -c -o $@ $< $(CFLAGS) -DHAVE_SINC -DSINC_LOWER_QUALITY + +sinc.o: ../sinc.c + $(CC) -c -o $@ $< $(CFLAGS) -DHAVE_SINC + +sinc-higher.o: ../sinc.c + $(CC) -c -o $@ $< $(CFLAGS) -DHAVE_SINC -DSINC_HIGHER_QUALITY + +sinc-highest.o: ../sinc.c + $(CC) -c -o $@ $< $(CFLAGS) -DHAVE_SINC -DSINC_HIGHEST_QUALITY + +test-sinc-lowest: sinc-lowest.o ../utils.o main.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-snr-sinc-lowest: sinc-lowest.o ../utils.o snr.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-sinc-lower: sinc-lower.o ../utils.o main.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-snr-sinc-lower: sinc-lower.o ../utils.o snr.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-sinc: sinc.o ../utils.o main.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-snr-sinc: sinc.o ../utils.o snr.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-sinc-higher: sinc-higher.o ../utils.o main.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-snr-sinc-higher: sinc-higher.o ../utils.o snr.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-sinc-highest: sinc-highest.o ../utils.o main.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + +test-snr-sinc-highest: sinc-highest.o ../utils.o snr.o ../hermite.o resampler-sinc.o + $(CC) -o $@ $^ $(LDFLAGS) + %.o: %.c $(CC) -c -o $@ $< $(CFLAGS) diff --git a/audio/test/sinc_test.m b/audio/test/sinc_test.m new file mode 100644 index 0000000000..24584cc5aa --- /dev/null +++ b/audio/test/sinc_test.m @@ -0,0 +1,38 @@ +% MATLAB test case for RetroArch SINC upsampler. +close all; + +% 4-tap and 8-tap are Lanczos windowed, but include here for completeness. +sidelobes = [2 4 8 32 128]; +taps = sidelobes * 2; + +phases = 256; +cutoffs = [0.65 0.75 0.825 0.90 0.95]; +betas = [2.0 3.0 5.5 10.5 14.5]; + +freqs = 0.05 : 0.05 : 0.99; + +filters = length(taps); +for i = 1 : filters + filter_length = taps(i) * phases; + + % Generate SINC. + sinc_indices = 2 * ((0 : (filter_length - 1)) / (filter_length - 1)) - 1; + s = cutoffs(i) * sinc(cutoffs(i) * sinc_indices * sidelobes(i)); + win = kaiser(filter_length, betas(i))'; + filter = s .* win; + + impulse_response_half = 0.5 * upfirdn(1, filter, phases, phases / 2); + figure('name', sprintf('Response SINC: %d taps', taps(i))); + freqz(impulse_response_half); + ylim([-200 0]); + + signal = zeros(1, 4001); + for freq = freqs + signal = signal + sin(pi * freq * (0 : 4000)); + end + + resampled = upfirdn(signal, filter, phases, round(phases * 44100 / 48000)); + figure('name', sprintf('Kaiser SINC: %d taps, w = %.f', taps(i), freq)); + freqz(resampled .* kaiser(length(resampled), 40.0)'); + ylim([-80 70]); +end \ No newline at end of file