Fix NaN in polyphase FIR generator

james34602 · Nov 22, 2023 · 5855bb6 · 5855bb6
1 parent be98515
commit 5855bb6
Show file tree

Hide file tree

Showing 7 changed files with 32 additions and 170 deletions.
diff --git a/Main/DSPManager/jni/libsamplerate/common.h b/Main/DSPManager/jni/libsamplerate/common.h
@@ -14,36 +14,6 @@
 #include <stdbool.h>
 #endif
 
-#if defined(__x86_64__) || defined(_M_X64)
-#   define HAVE_SSE2_INTRINSICS
-#elif defined(ENABLE_SSE2_LRINT) && (defined(_M_IX86) || defined(__i386__))
-#   if defined(_MSC_VER)
-#       define HAVE_SSE2_INTRINSICS
-#   elif defined(__clang__)
-#       ifdef __SSE2__
-#           define HAVE_SSE2_INTRINSICS
-#       elif (__has_attribute(target))
-#           define HAVE_SSE2_INTRINSICS
-#           define USE_TARGET_ATTRIBUTE
-#       endif
-#   elif defined(__GNUC__)
-#       ifdef __SSE2__
-#           define HAVE_SSE2_INTRINSICS
-#       elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 9))
-#           define HAVE_SSE2_INTRINSICS
-#           define USE_TARGET_ATTRIBUTE
-#       endif
-#   endif
-#endif
-
-#ifdef HAVE_SSE2_INTRINSICS
-#ifdef HAVE_IMMINTRIN_H
-#include <immintrin.h>
-#else
-#include <emmintrin.h>
-#endif
-#endif /* HAVE_SSE2_INTRINSICS */
-
 #include <math.h>
 
 #ifdef HAVE_VISIBILITY
@@ -185,40 +155,15 @@ const char* linear_get_description (int src_enum) ;
 
 SRC_STATE *linear_state_new (int channels, SRC_ERROR *error) ;
 
-/*----------------------------------------------------------
-** SIMD optimized math functions.
-*/
-
-#ifdef HAVE_SSE2_INTRINSICS
-static inline int
-#ifdef USE_TARGET_ATTRIBUTE
-__attribute__((target("sse2")))
-#endif
-psf_lrintf (float x)
+static inline int psf_lrintf(float x)
 {
-	return _mm_cvtss_si32 (_mm_load_ss (&x)) ;
-}
-static inline int
-#ifdef USE_TARGET_ATTRIBUTE
-__attribute__((target("sse2")))
-#endif
-psf_lrint (double x)
-{
-	return _mm_cvtsd_si32 (_mm_load_sd (&x)) ;
-}
-
-#else
-
-static inline int psf_lrintf (float x)
-{
-	return lrintf (x) ;
+	return lrintf(x);
 } /* psf_lrintf */
 
-static inline int psf_lrint (double x)
+static inline int psf_lrint(double x)
 {
-	return lrint (x) ;
+	return lrint(x);
 } /* psf_lrint */
-#endif
 
 /*----------------------------------------------------------
 **	Common static inline functions.
@@ -241,5 +186,4 @@ is_bad_src_ratio (double ratio)
 } /* is_bad_src_ratio */
 
 
-#endif	/* COMMON_H_INCLUDED */
-
+#endif	/* COMMON_H_INCLUDED */
diff --git a/Main/DSPManager/jni/libsamplerate/src_sinc.c b/Main/DSPManager/jni/libsamplerate/src_sinc.c
@@ -29,24 +29,6 @@ double sinc(double x)
 	else
 		return 1.0;
 }
-double bessi0(double x)
-{
-	double xh, sum, pow, ds;
-	int k;
-	xh = 0.5 * x;
-	sum = 1.0;
-	pow = 1.0;
-	k = 0;
-	ds = 1.0;
-	while (ds > sum * DBL_EPSILON)
-	{
-		++k;
-		pow = pow * (xh / k);
-		ds = pow * pow;
-		sum = sum + ds;
-	}
-	return sum;
-}
 double chbevl(double x, double array[], int n)
 {
 	double b0, b1, b2, *p;
@@ -156,7 +138,10 @@ unsigned int getcoeff(double ratio, double atten, unsigned long long N, unsigned
 	y[0] = 1.0;
 	for (m = 1; m < ((N - 1) >> 1) + 1; m++)
 	{
-		double w = i0(bta * sqrt(1.0 - (double)(4 * m * m) * xind)) * bes;
+		double v = 1.0 - (double)(4 * m * m) * xind;
+		if (v < 0.0)
+			v = 0.0;
+		double w = i0(bta * sqrt(v)) * bes;
 		double f = sinc(m / ratio) * w;
 		y[m] = f;
 		sum += f;

diff --git a/Main/libjamesdsp/jni/jamesdsp/Android.mk b/Main/libjamesdsp/jni/jamesdsp/Android.mk
@@ -91,17 +91,17 @@ LOCAL_SRC_FILES := \
 # terminator
 LOCAL_LDLIBS := -llog
 ifeq ($(TARGET_ARCH_ABI), armeabi-v7a)
-LOCAL_CPPFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv7-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
-LOCAL_CFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv7-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
+LOCAL_CPPFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv7-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
+LOCAL_CFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv7-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
 else ifeq ($(TARGET_ARCH_ABI), arm64-v8a)
-LOCAL_CPPFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv8-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
-LOCAL_CFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv8-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
+LOCAL_CPPFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv8-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
+LOCAL_CFLAGS += -Wall -Wextra -ffunction-sections -fdata-sections -Ofast -march=armv8-a -mfpu=neon -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
 else ifeq ($(TARGET_ARCH_ABI), x86)
-LOCAL_CPPFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
-LOCAL_CFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
+LOCAL_CPPFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
+LOCAL_CFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
 else ifeq ($(TARGET_ARCH_ABI), armeabi)
-LOCAL_CPPFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
-LOCAL_CFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL -DDEBUG
+LOCAL_CPPFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
+LOCAL_CFLAGS += -ffunction-sections -fdata-sections -Ofast -ftree-vectorize -fvisibility=hidden -DJAMESDSP_REFERENCE_IMPL #-DDEBUG -g
 endif
 LOCAL_LDFLAGS += -Wl,--gc-sections,--exclude-libs,ALL
 include $(BUILD_SHARED_LIBRARY)
diff --git a/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/crossfeed.c b/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/crossfeed.c
@@ -18,16 +18,19 @@ void CrossfeedDestructor(JamesDSPLib *jdsp)
 		free(jdsp->advXF.conv[0]);
 		free(jdsp->advXF.conv[1]);
 		free(jdsp->advXF.conv[2]);
+		jdsp->advXF.conv[0] = 0;
 	}
 	if (jdsp->advXF.convLong_T_S)
 	{
 		TwoStageFFTConvolver2x4x2Free(jdsp->advXF.convLong_T_S);
 		free(jdsp->advXF.convLong_T_S);
+		jdsp->advXF.convLong_T_S = 0;
 	}
 	if (jdsp->advXF.convLong_S_S)
 	{
 		FFTConvolver2x4x2Free(jdsp->advXF.convLong_S_S);
 		free(jdsp->advXF.convLong_S_S);
+		jdsp->advXF.convLong_S_S = 0;
 	}
 }
 void CrossfeedProcessFFTConvolver2x4x2(JamesDSPLib *jdsp, size_t n)
@@ -42,25 +45,7 @@ void CrossfeedEnable(JamesDSPLib *jdsp, char enable)
 {
 	if (jdsp->crossfeedForceRefresh || !jdsp->advXF.conv[0] || !jdsp->advXF.conv[1] || !jdsp->advXF.conv[2] || !jdsp->advXF.convLong_S_S || !jdsp->advXF.convLong_T_S)
 	{
-		if (jdsp->advXF.conv[0])
-		{
-			FFTConvolver2x4x2Free(jdsp->advXF.conv[0]);
-			FFTConvolver2x4x2Free(jdsp->advXF.conv[1]);
-			FFTConvolver2x4x2Free(jdsp->advXF.conv[2]);
-			free(jdsp->advXF.conv[0]);
-			free(jdsp->advXF.conv[1]);
-			free(jdsp->advXF.conv[2]);
-		}
-		if (jdsp->advXF.convLong_T_S)
-		{
-			TwoStageFFTConvolver2x4x2Free(jdsp->advXF.convLong_T_S);
-			free(jdsp->advXF.convLong_T_S);
-		}
-		if (jdsp->advXF.convLong_S_S)
-		{
-			FFTConvolver2x4x2Free(jdsp->advXF.convLong_S_S);
-			free(jdsp->advXF.convLong_S_S);
-		}
+		CrossfeedDestructor(jdsp);
 		jdsp->advXF.conv[0] = (FFTConvolver2x4x2 *)malloc(sizeof(FFTConvolver2x4x2));
 		FFTConvolver2x4x2Init(jdsp->advXF.conv[0]);
 		FFTConvolver2x4x2LoadImpulseResponse(jdsp->advXF.conv[0], (unsigned int)jdsp->blockSize, jdsp->blobsCh1[0], jdsp->blobsCh2[0], jdsp->blobsCh3[0], jdsp->blobsCh4[0], jdsp->blobsResampledLen);

diff --git a/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/dynamic.c b/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/dynamic.c
@@ -1656,7 +1656,7 @@ void CompressorSetGain(JamesDSPLib *jdsp, double *freq, double *gains, char cpy)
 	cm->gains2[0] = cm->gains2[1];
 	cm->freq2[NUMPTS_DRS + 1] = 24000.0;
 	cm->gains2[NUMPTS_DRS + 1] = cm->gains2[NUMPTS_DRS];
-	makima(&cm->pch, cm->freq2, cm->gains2, NUMPTS_DRS + 2, 1, 1);
+	pchip(&cm->pch, cm->freq2, cm->gains2, NUMPTS_DRS + 2, 1, 1);
 	unsigned int specLen = *((unsigned int *)(cm->octaveSmooth));
 	float reciprocal = *((float *)(cm->octaveSmooth + sizeof(unsigned int)));
 	unsigned int lpLen = *((unsigned int *)(cm->octaveSmooth + sizeof(unsigned int) + sizeof(float)));

diff --git a/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/eel2/numericSys/libsamplerate/common.h b/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/eel2/numericSys/libsamplerate/common.h
@@ -14,36 +14,6 @@
 #include <stdbool.h>
 #endif
 
-#if defined(__x86_64__) || defined(_M_X64)
-#   define HAVE_SSE2_INTRINSICS
-#elif defined(ENABLE_SSE2_LRINT) && (defined(_M_IX86) || defined(__i386__))
-#   if defined(_MSC_VER)
-#       define HAVE_SSE2_INTRINSICS
-#   elif defined(__clang__)
-#       ifdef __SSE2__
-#           define HAVE_SSE2_INTRINSICS
-#       elif (__has_attribute(target))
-#           define HAVE_SSE2_INTRINSICS
-#           define USE_TARGET_ATTRIBUTE
-#       endif
-#   elif defined(__GNUC__)
-#       ifdef __SSE2__
-#           define HAVE_SSE2_INTRINSICS
-#       elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 9))
-#           define HAVE_SSE2_INTRINSICS
-#           define USE_TARGET_ATTRIBUTE
-#       endif
-#   endif
-#endif
-
-#ifdef HAVE_SSE2_INTRINSICS
-#ifdef HAVE_IMMINTRIN_H
-#include <immintrin.h>
-#else
-#include <emmintrin.h>
-#endif
-#endif /* HAVE_SSE2_INTRINSICS */
-
 #include <math.h>
 
 #ifdef HAVE_VISIBILITY
@@ -185,40 +155,15 @@ const char* linear_get_description (int src_enum) ;
 
 SRC_STATE *linear_state_new (int channels, SRC_ERROR *error) ;
 
-/*----------------------------------------------------------
-** SIMD optimized math functions.
-*/
-
-#ifdef HAVE_SSE2_INTRINSICS
-static inline int
-#ifdef USE_TARGET_ATTRIBUTE
-__attribute__((target("sse2")))
-#endif
-psf_lrintf (float x)
+static inline int psf_lrintf(float x)
 {
-	return _mm_cvtss_si32 (_mm_load_ss (&x)) ;
-}
-static inline int
-#ifdef USE_TARGET_ATTRIBUTE
-__attribute__((target("sse2")))
-#endif
-psf_lrint (double x)
-{
-	return _mm_cvtsd_si32 (_mm_load_sd (&x)) ;
-}
-
-#else
-
-static inline int psf_lrintf (float x)
-{
-	return lrintf (x) ;
+	return lrintf(x);
 } /* psf_lrintf */
 
-static inline int psf_lrint (double x)
+static inline int psf_lrint(double x)
 {
-	return lrint (x) ;
+	return lrint(x);
 } /* psf_lrint */
-#endif
 
 /*----------------------------------------------------------
 **	Common static inline functions.
@@ -241,5 +186,4 @@ is_bad_src_ratio (double ratio)
 } /* is_bad_src_ratio */
 
 
-#endif	/* COMMON_H_INCLUDED */
-
+#endif	/* COMMON_H_INCLUDED */
diff --git a/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/eel2/numericSys/libsamplerate/src_sinc.c b/Main/libjamesdsp/jni/jamesdsp/jdsp/Effects/eel2/numericSys/libsamplerate/src_sinc.c
@@ -1,5 +1,6 @@
 /*
 ** Copyright (c) 2002-2021, Erik de Castro Lopo <[email protected]>
+** Copyright (c) 2023, James Fung <[email protected]>
 ** All rights reserved.
 **
 ** This code is released under 2-clause BSD license. Please see the
@@ -132,7 +133,10 @@ unsigned int getcoeff(double ratio, double atten, unsigned long long N, unsigned
 	y[0] = 1.0;
 	for (m = 1; m < ((N - 1) >> 1) + 1; m++)
 	{
-		double w = i0(bta * sqrt(1.0 - (double)(4 * m * m) * xind)) * bes;
+		double v = 1.0 - (double)(4 * m * m) * xind;
+		if (v < 0.0)
+			v = 0.0;
+		double w = i0(bta * sqrt(v)) * bes;
 		double f = sinc(m / ratio) * w;
 		y[m] = f;
 		sum += f;