resampler: better bessel function
Applied fix from RBJ to Bessel polynomial.
Benchmarked various implementations.
Chose HyperbolicCosineWindow because it is faster to compute and cleaner.
diff --git a/src/flowgraph/resampler/BesselFunctionRBJ.h b/src/flowgraph/resampler/BesselFunctionRBJ.h
index 8dc6aaa..de6200b 100644
--- a/src/flowgraph/resampler/BesselFunctionRBJ.h
+++ b/src/flowgraph/resampler/BesselFunctionRBJ.h
@@ -31,12 +31,12 @@
* Higher K gives higher precision.
* @param K number of terms in the polynomial expansion
*/
- BesselFunctionRBJ(int K = 8)
+ BesselFunctionRBJ(int K = 20)
: mK(K) {
mCoefficients = std::make_unique<double[]>(K + 1);
mCoefficients[0] = 1.0;
for (int64_t k = 1; k <= K; k++) {
- mCoefficients[k] = mCoefficients[k - 1] * -0.25 / (k * k);
+ mCoefficients[k] = mCoefficients[k - 1] * 0.25 / (k * k);
}
}
@@ -53,8 +53,8 @@
}
// Direct summation of the first K terms.
- double direct(double x, int K) {
- double z = x * x * -0.25;
+ double direct(double x, int K = 20) {
+ double z = x * x * 0.25;
double sum = 0.0;
for (int k = K; k >= 1; k--) {
int n = k;
diff --git a/src/flowgraph/resampler/KaiserWindow.h b/src/flowgraph/resampler/KaiserWindow.h
index f45aced..9232917 100644
--- a/src/flowgraph/resampler/KaiserWindow.h
+++ b/src/flowgraph/resampler/KaiserWindow.h
@@ -17,11 +17,10 @@
#ifndef RESAMPLER_KAISER_WINDOW_H
#define RESAMPLER_KAISER_WINDOW_H
-#include "Math.h"
+#include <math.h>
namespace resampler {
-
/**
* Calculate a Kaiser window centered at 0.
*/
@@ -62,6 +61,10 @@
return bessel(w) * mInverseBesselBeta;
}
+ // Approximation of a
+ // modified zero order Bessel function of the first kind.
+ // Based on:
+ // https://dsp.stackexchange.com/questions/37714/kaiser-window-approximation
static double bessel(double x) {
double y = cosh(0.970941817426052 * x);
y += cosh(0.8854560256532099 * x);
diff --git a/src/flowgraph/resampler/MultiChannelResampler.cpp b/src/flowgraph/resampler/MultiChannelResampler.cpp
index 2eca588..c35cf04 100644
--- a/src/flowgraph/resampler/MultiChannelResampler.cpp
+++ b/src/flowgraph/resampler/MultiChannelResampler.cpp
@@ -112,12 +112,6 @@
}
}
-float MultiChannelResampler::hammingWindow(float radians, float spread) {
- const float alpha = 0.54f;
- const float windowPhase = radians / spread;
- return (float) (alpha + ((1.0 - alpha) * cosf(windowPhase)));
-}
-
float MultiChannelResampler::sinc(float radians) {
if (abs(radians) < 1.0e-9) return 1.0f; // avoid divide by zero
return sinf(radians) / radians; // Sinc function
@@ -139,14 +133,19 @@
? ((float)outputRate / inputRate)
: ((float)inputRate / outputRate));
const int numTapsHalf = getNumTaps() / 2; // numTaps must be even.
+ const float numTapsHalfInverse = 1.0f / numTapsHalf;
for (int i = 0; i < numRows; i++) {
float tapPhase = phase - numTapsHalf;
float gain = 0.0; // sum of raw coefficients
int gainCursor = coefficientIndex;
for (int tap = 0; tap < getNumTaps(); tap++) {
float radians = tapPhase * M_PI;
- float window = mCoshWindow(tapPhase / numTapsHalf);
- // float window = hammingWindow(radians, numTapsHalf);
+
+#if MCR_USE_KAISER
+ float window = mKaiserWindow(tapPhase * numTapsHalfInverse);
+#else
+ float window = mCoshWindow(tapPhase * numTapsHalfInverse);
+#endif
float coefficient = sinc(radians * cutoffScaler) * window;
mCoefficients.at(coefficientIndex++) = coefficient;
gain += coefficient;
@@ -161,7 +160,6 @@
float gainCorrection = 1.0 / gain; // normalize the gain
for (int tap = 0; tap < getNumTaps(); tap++) {
mCoefficients.at(gainCursor + tap) *= gainCorrection;
-// printf("%d, %8.5f\n", tap, mCoefficients.at(gainCursor + tap));
}
}
}
diff --git a/src/flowgraph/resampler/MultiChannelResampler.h b/src/flowgraph/resampler/MultiChannelResampler.h
index 932af41..31ff2c0 100644
--- a/src/flowgraph/resampler/MultiChannelResampler.h
+++ b/src/flowgraph/resampler/MultiChannelResampler.h
@@ -22,7 +22,15 @@
#include <sys/types.h>
#include <unistd.h>
+#ifndef MCR_USE_KAISER
+#define MCR_USE_KAISER 0
+#endif
+
+#if MCR_USE_KAISER
+#include "KaiserWindow.h"
+#else
#include "HyperbolicCosineWindow.h"
+#endif
namespace resampler {
@@ -191,12 +199,6 @@
explicit MultiChannelResampler(const MultiChannelResampler::Builder &builder);
/**
- * @param phase between 0.0 and 2*spread // TODO use centered phase, maybe
- * @return windowedSinc
- */
- // static float calculateWindowedSinc(float phase, int spread); // TODO remove
-
- /**
* Write a frame containing N samples.
* Call advanceWrite() after calling this.
* @param frame pointer to the first sample in a frame
@@ -251,7 +253,11 @@
private:
+#if MCR_USE_KAISER
+ KaiserWindow mKaiserWindow;
+#else
HyperbolicCosineWindow mCoshWindow;
+#endif
// max coefficients for polyphase filter
static constexpr float kDefaultNormalizedCutoff = 0.70f;
