resampler: use more accurate bessel function for Kaiser window
diff --git a/src/flowgraph/resampler/BesselFunctionRBJ.h b/src/flowgraph/resampler/BesselFunctionRBJ.h
new file mode 100644
index 0000000..5605db5
--- /dev/null
+++ b/src/flowgraph/resampler/BesselFunctionRBJ.h
@@ -0,0 +1,59 @@
+/*
+ * Copyright 2019 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef RESAMPLER_BESSEL_FUNCTION_H
+#define RESAMPLER_BESSEL_FUNCTION_H
+
+#include <memory>
+
+/**
+ * Approximate a zero order Bessel function of the first kind.
+ * Use Horner's method to retain numeric precision.
+ *
+ * Based on a posting by Robert Bristow-Johnson at:
+ * https://dsp.stackexchange.com/questions/37714/kaiser-window-approximation
+ */
+template <class T>
+class BesselFunctionRBJ {
+public:
+ /**
+ * Higher K gives higher precision.
+ * @param K number of terms in the polynomial expansion
+ */
+ BesselFunctionRBJ(int K = 8)
+ : mK(K) {
+ mCoefficients = std::make_unique<T[]>(K + 1);
+ mCoefficients[0] = 1.0;
+ for (int64_t k = 1; k <= K; k++) {
+ mCoefficients[k] = mCoefficients[k - 1] * -1 / ((k * k) * 4);
+ }
+ }
+
+ T operator()(T x) {
+ T x2 = x * x;
+ T y = x2 * mCoefficients[mK];
+ for (int k = mK - 1; k >= 0; k--) {
+ y = mCoefficients[k] + x2 * y;
+ }
+ return y;
+ }
+
+private:
+ std::unique_ptr<T[]> mCoefficients;
+ int mK;
+};
+
+#endif //RESAMPLER_BESSEL_FUNCTION_H
diff --git a/src/flowgraph/resampler/KaiserWindow.h b/src/flowgraph/resampler/KaiserWindow.h
new file mode 100644
index 0000000..6a5d086
--- /dev/null
+++ b/src/flowgraph/resampler/KaiserWindow.h
@@ -0,0 +1,64 @@
+/*
+ * Copyright 2019 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef RESAMPLER_KAISER_WINDOW_H
+#define RESAMPLER_KAISER_WINDOW_H
+
+#include "Math.h"
+
+/**
+ * Calculate a Kaiser window centered at 0.
+ */
+class KaiserWindow {
+public:
+ KaiserWindow() {
+ setStopBandAttenuation(90);
+ }
+
+ /**
+ * @param attenuation typical values range from 30 to 90 dB
+ */
+ double setStopBandAttenuation(double attenuation) {
+ double beta = 0.0;
+ if (attenuation > 50) {
+ beta = 0.1102 * (attenuation - 8.7);
+ } else if (attenuation >= 21) {
+ double a21 = attenuation - 21;
+ beta = 0.5842 * pow(a21, 0.4) + (0.07886 * a21);
+ }
+ setBeta(beta);
+ return beta;
+ }
+
+ void setBeta(double beta) {
+ mBeta = beta;
+ mInverseBesselBeta = 1.0 / Math::bessel(beta);
+ }
+
+ /**
+ * @param x ranges from -1.0 to +1.0
+ */
+ double operator()(double x) {
+ double w = mBeta * sqrt(1.0 - (x * x));
+ return Math::bessel(w) * mInverseBesselBeta;
+ }
+
+private:
+ double mBeta = 0.0;
+ double mInverseBesselBeta = 1.0;
+};
+
+#endif //RESAMPLER_KAISER_WINDOW_H
diff --git a/src/flowgraph/resampler/Math.h b/src/flowgraph/resampler/Math.h
new file mode 100644
index 0000000..06bd7bd
--- /dev/null
+++ b/src/flowgraph/resampler/Math.h
@@ -0,0 +1,36 @@
+/*
+ * Copyright 2019 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef RESAMPLER_MATH_H
+#define RESAMPLER_MATH_H
+
+class Math {
+public:
+static double bessel(double x) {
+ double y = cosh(0.970941817426052 * x);
+ y += cosh(0.8854560256532099 * x);
+ y += cosh(0.7485107481711011 * x);
+ y += cosh(0.5680647467311558 * x);
+ y += cosh(0.3546048870425356 * x);
+ y += cosh(0.120536680255323 * x);
+ y *= 2;
+ y += cosh(x);
+ y /= 13;
+ return y;
+}
+
+};
+#endif //RESAMPLER_MATH_H
diff --git a/src/flowgraph/resampler/MultiChannelResampler.cpp b/src/flowgraph/resampler/MultiChannelResampler.cpp
index e461975..cbd0693 100644
--- a/src/flowgraph/resampler/MultiChannelResampler.cpp
+++ b/src/flowgraph/resampler/MultiChannelResampler.cpp
@@ -112,23 +112,17 @@
}
}
-float MultiChannelResampler::hammingWindow(float radians, int spread) {
+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) < 0.00000001) return 1.0f; // avoid divide by zero
+ if (abs(radians) < 1.0e-9) return 1.0f; // avoid divide by zero
return sinf(radians) / radians; // Sinc function
}
-// Unoptimized calculation used to construct lookup tables.
-//float MultiChannelResampler::calculateWindowedSinc(float radians, int spread) {
-// return sinc(radians) * hammingWindow(radians, spread); // TODO try Kaiser window
-//}
-
-
// Generate coefficients in the order they will be used by readFrame().
// This is more complicated but readFrame() is called repeatedly and should be optimized.
void MultiChannelResampler::generateCoefficients(int32_t inputRate,
diff --git a/src/flowgraph/resampler/MultiChannelResampler.h b/src/flowgraph/resampler/MultiChannelResampler.h
index d4e3aa4..f903d31 100644
--- a/src/flowgraph/resampler/MultiChannelResampler.h
+++ b/src/flowgraph/resampler/MultiChannelResampler.h
@@ -180,7 +180,7 @@
return mChannelCount;
}
- static float hammingWindow(float radians, int spread);
+ static float hammingWindow(float radians, float spread);
static float sinc(float radians);
