diff --git a/cmake/SfizzDeps.cmake b/cmake/SfizzDeps.cmake index 1348b558..3dcb1916 100644 --- a/cmake/SfizzDeps.cmake +++ b/cmake/SfizzDeps.cmake @@ -96,6 +96,12 @@ target_include_directories(sfizz_kissfft PUBLIC "src/external/kiss_fft" PUBLIC "src/external/kiss_fft/tools") +# The cephes library +add_library(sfizz_cephes STATIC + "external/cephes/src/chbevl.c" + "external/cephes/src/i0.c") +add_library(sfizz::cephes ALIAS sfizz_cephes) + # The cpuid library add_library(sfizz_cpuid STATIC "src/external/cpuid/src/cpuid/cpuinfo.cpp" diff --git a/external/cephes/LICENSE.txt b/external/cephes/LICENSE.txt new file mode 100644 index 00000000..666ddb13 --- /dev/null +++ b/external/cephes/LICENSE.txt @@ -0,0 +1,119 @@ +==== NOTE ==== +The actual cephes library, shipped with and wrapped by this package, is available on The Netlib at http://www.netlib.org/cephes/ . It does not have any license specified. However, its original authors, Stephen Moshier, has kindly granted permission for inclusion in a BSD-licensed package. See email snippet below for reference. + +Return-Path: +X-Original-To: julien@cornebise.com +Delivered-To: julien@cornebise.com +Received: from atl4mhob11.myregisteredsite.com (atl4mhob11.myregisteredsite.com [209.17.115.49]) + by cornebise.com (Postfix) with ESMTP id D47B139FC0 + for ; Fri, 25 Oct 2013 16:32:40 +0200 (CEST) +Received: from mailpod1.hostingplatform.com ([10.30.71.116]) + by atl4mhob11.myregisteredsite.com (8.14.4/8.14.4) with ESMTP id r9PEWcwQ003543 + for ; Fri, 25 Oct 2013 10:32:38 -0400 +Received: (qmail 11948 invoked by uid 0); 25 Oct 2013 12:36:20 -0000 +X-TCPREMOTEIP: 76.24.25.74 +X-Authenticated-UID: steve@moshier.net +Received: from unknown (HELO d510.local) (steve@moshier.net@76.24.25.74) + by 0 with ESMTPA; 25 Oct 2013 12:36:20 -0000 +Date: Fri, 25 Oct 2013 08:36:19 -0400 (EDT) +From: Stephen Moshier +X-X-Sender: steve@d510 +To: Julien Cornebise +Subject: Re: Cephes: permission to wrap+distribute for Lua +In-Reply-To: <52653AD3.1010004@cornebise.com> +Message-ID: +References: <52653AD3.1010004@cornebise.com> +User-Agent: Alpine 2.02 (DEB 1266 2009-07-14) +MIME-Version: 1.0 +Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed + + +Julien, thank you for writing. +BSD license is fine, modification is OK. +There are more build scripts available in the web site distributions than +there are on the Netlib. I think there is an update to Planck's radiation +function that I haven't sent to Netlib yet. But Netlib is a more stable +site, so it is better to cite that as a reference. + + +On Mon, 21 Oct 2013, Julien Cornebise wrote: + +> -----BEGIN PGP SIGNED MESSAGE----- +> Hash: SHA1 +> +> Dear Mr Moshier +> +> I am a researcher in mathematics and machine learning in London, and +> am writing about your awesome Cephes library, whom I found at the +> heart of Scipy. +> +> It is so useful that, with your permission, I would like to wrap it +> for Lua and Torch (a machine learning overlay to Lua, specialized in +> neural nets, see http://www.torch.ch). I would like to distribute it +> as a package for Torch, including your source code along the wrapping +> code. +> This wouldbe a public package, distributed under BSD License. I have +> put a first draft on github: +> https://github.com/jucor/torch-cephes +> +> Hence my three questions, please: +> +> 1/ How would you like to be acknowledged, beyond the comments that are +> already in your code? Do you have any standard header/disclaimer that +> I could add to the documentation? +> +> 2/ At the moment, your code is left untouched. However, if I ever need +> to modify bits of the code, what are the conditions/restrictions? +> Nothing huge -- I definitely do not want to mess with it: I was +> planning to use the natural completion of some functions on the +> completed real line (e.g. CDF returing 1 when called with "infinity", +> or quantiles returning -Infinity when called with 0), either natively +> if supported, or by setting a specific flag via mtherr(). +> +> 3/ I am currently using the source from Netlib. Do you recommend using +> the source from your website instead ? +> +> Thank you very much for your attention, +> and, more importantly, for the time and effort your poured into Cephes. +> +> Best regards, +> +> Julien Cornebise, Ph.D. +> London, UK +> http://www.cornebise.com/julien +> -----BEGIN PGP SIGNATURE----- +> Version: GnuPG v1.4.14 (Darwin) +> Comment: GPGTools - http://gpgtools.org +> Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/ +> +> iEYEARECAAYFAlJlOtEACgkQKYR3gC0rw/gIpQCfZKu6+iDh9ghhm6QfsLXnldKN +> BuIAn2zZHu1c/IrRAevhjM7N7xGg0LHO +> =WeP5 +> -----END PGP SIGNATURE----- + + +==== LICENSE ==== +Copyright (c) 2013, Julien Cornebise +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the organization nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/external/cephes/src/chbevl.c b/external/cephes/src/chbevl.c new file mode 100644 index 00000000..53938816 --- /dev/null +++ b/external/cephes/src/chbevl.c @@ -0,0 +1,82 @@ +/* chbevl.c + * + * Evaluate Chebyshev series + * + * + * + * SYNOPSIS: + * + * int N; + * double x, y, coef[N], chebevl(); + * + * y = chbevl( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates the series + * + * N-1 + * - ' + * y = > coef[i] T (x/2) + * - i + * i=0 + * + * of Chebyshev polynomials Ti at argument x/2. + * + * Coefficients are stored in reverse order, i.e. the zero + * order term is last in the array. Note N is the number of + * coefficients, not the order. + * + * If coefficients are for the interval a to b, x must + * have been transformed to x -> 2(2x - b - a)/(b-a) before + * entering the routine. This maps x from (a, b) to (-1, 1), + * over which the Chebyshev polynomials are defined. + * + * If the coefficients are for the inverted interval, in + * which (a, b) is mapped to (1/b, 1/a), the transformation + * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity, + * this becomes x -> 4a/x - 1. + * + * + * + * SPEED: + * + * Taking advantage of the recurrence properties of the + * Chebyshev polynomials, the routine requires one more + * addition per loop than evaluating a nested polynomial of + * the same degree. + * + */ + /* chbevl.c */ + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1985, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +double chbevl( x, array, n ) +double x; +double array[]; +int n; +{ +double b0, b1, b2, *p; +int i; + +p = array; +b0 = *p++; +b1 = 0.0; +i = n - 1; + +do + { + b2 = b1; + b1 = b0; + b0 = x * b1 - b2 + *p++; + } +while( --i ); + +return( 0.5*(b0-b2) ); +} diff --git a/external/cephes/src/i0.c b/external/cephes/src/i0.c new file mode 100644 index 00000000..307398e9 --- /dev/null +++ b/external/cephes/src/i0.c @@ -0,0 +1,193 @@ +/* i0.c + * + * Modified Bessel function of order zero + * + * + * + * SYNOPSIS: + * + * double x, y, i0(); + * + * y = i0( x ); + * + * + * + * DESCRIPTION: + * + * Returns modified Bessel function of order zero of the + * argument. + * + * The function is defined as i0(x) = j0( ix ). + * + * The range is partitioned into the two intervals [0,8] and + * (8, infinity). Chebyshev polynomial expansions are employed + * in each interval. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0,30 6000 8.2e-17 1.9e-17 + * IEEE 0,30 30000 5.8e-16 1.4e-16 + * + */ + /* i0e.c + * + * Modified Bessel function of order zero, + * exponentially scaled + * + * + * + * SYNOPSIS: + * + * double x, y, i0e(); + * + * y = i0e( x ); + * + * + * + * DESCRIPTION: + * + * Returns exponentially scaled modified Bessel function + * of order zero of the argument. + * + * The function is defined as i0e(x) = exp(-|x|) j0( ix ). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 30000 5.4e-16 1.2e-16 + * See i0(). + * + */ + +/* i0.c */ + + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1987, 2000 by Stephen L. Moshier +*/ + +#include + +/* Chebyshev coefficients for exp(-x) I0(x) + * in the interval [0,8]. + * + * lim(x->0){ exp(-x) I0(x) } = 1. + */ + +static double A[] = +{ +-4.41534164647933937950E-18, + 3.33079451882223809783E-17, +-2.43127984654795469359E-16, + 1.71539128555513303061E-15, +-1.16853328779934516808E-14, + 7.67618549860493561688E-14, +-4.85644678311192946090E-13, + 2.95505266312963983461E-12, +-1.72682629144155570723E-11, + 9.67580903537323691224E-11, +-5.18979560163526290666E-10, + 2.65982372468238665035E-9, +-1.30002500998624804212E-8, + 6.04699502254191894932E-8, +-2.67079385394061173391E-7, + 1.11738753912010371815E-6, +-4.41673835845875056359E-6, + 1.64484480707288970893E-5, +-5.75419501008210370398E-5, + 1.88502885095841655729E-4, +-5.76375574538582365885E-4, + 1.63947561694133579842E-3, +-4.32430999505057594430E-3, + 1.05464603945949983183E-2, +-2.37374148058994688156E-2, + 4.93052842396707084878E-2, +-9.49010970480476444210E-2, + 1.71620901522208775349E-1, +-3.04682672343198398683E-1, + 6.76795274409476084995E-1 +}; + + +/* Chebyshev coefficients for exp(-x) sqrt(x) I0(x) + * in the inverted interval [8,infinity]. + * + * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). + */ + +static double B[] = +{ +-7.23318048787475395456E-18, +-4.83050448594418207126E-18, + 4.46562142029675999901E-17, + 3.46122286769746109310E-17, +-2.82762398051658348494E-16, +-3.42548561967721913462E-16, + 1.77256013305652638360E-15, + 3.81168066935262242075E-15, +-9.55484669882830764870E-15, +-4.15056934728722208663E-14, + 1.54008621752140982691E-14, + 3.85277838274214270114E-13, + 7.18012445138366623367E-13, +-1.79417853150680611778E-12, +-1.32158118404477131188E-11, +-3.14991652796324136454E-11, + 1.18891471078464383424E-11, + 4.94060238822496958910E-10, + 3.39623202570838634515E-9, + 2.26666899049817806459E-8, + 2.04891858946906374183E-7, + 2.89137052083475648297E-6, + 6.88975834691682398426E-5, + 3.36911647825569408990E-3, + 8.04490411014108831608E-1 +}; + + +extern double chbevl ( double, void *, int ); + +double i0(x) +double x; +{ +double y; + +if( x < 0 ) + x = -x; +if( x <= 8.0 ) + { + y = (x/2.0) - 2.0; + return( exp(x) * chbevl( y, A, 30 ) ); + } + +return( exp(x) * chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) ); + +} + + + + +double i0e( x ) +double x; +{ +double y; + +if( x < 0 ) + x = -x; +if( x <= 8.0 ) + { + y = (x/2.0) - 2.0; + return( chbevl( y, A, 30 ) ); + } + +return( chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) ); + +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index df68a26b..47c24210 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -242,7 +242,7 @@ target_sources(sfizz_internal PRIVATE ${SFIZZ_HEADERS} ${SFIZZ_SOURCES} ${FAUST_ target_include_directories(sfizz_internal PUBLIC "." "sfizz") target_link_libraries(sfizz_internal PUBLIC absl::strings absl::span sfizz::filesystem sfizz::atomic_queue - PRIVATE sfizz::parser sfizz::messaging absl::flat_hash_map Threads::Threads st_audiofile sfizz::pugixml sfizz::spline sfizz::tunings sfizz::hiir sfizz::kissfft sfizz::cpuid sfizz::threadpool sfizz::jsl sfizz::atomic) + PRIVATE sfizz::parser sfizz::messaging absl::flat_hash_map Threads::Threads st_audiofile sfizz::pugixml sfizz::spline sfizz::tunings sfizz::hiir sfizz::kissfft sfizz::cephes sfizz::cpuid sfizz::threadpool sfizz::jsl sfizz::atomic) if(SFIZZ_USE_SNDFILE) target_compile_definitions(sfizz_internal PUBLIC "SFIZZ_USE_SNDFILE=1") target_link_libraries(sfizz_internal PUBLIC st_audiofile) diff --git a/src/sfizz/MathHelpers.h b/src/sfizz/MathHelpers.h index 50b012bd..d04793fc 100644 --- a/src/sfizz/MathHelpers.h +++ b/src/sfizz/MathHelpers.h @@ -25,6 +25,13 @@ #include #endif +#if __cplusplus >= 201703L +static double i0(double x) { return std::cyl_bessel_i(0.0, x); } +#else +// external Bessel function from cephes +extern "C" double i0(double x); +#endif + template constexpr T max(T op1, T op2) { @@ -465,6 +472,54 @@ bool isReasonableAudio(absl::Span span) return true; } +/** + * @brief Compute the Kaiser window + * + * @param b Kaiser parameter beta + * @param window Span of real which receives the window + */ +template +void kaiserWindow(double b, absl::Span window) +{ + double i0b = i0(b); + for (size_t i = 0, n = window.size(); i < n; ++i) { + double x = i / static_cast(n - 1); + double t = x + x - 1.0; + window[i] = static_cast(i0(b * std::sqrt(1.0 - t * t)) / i0b); + } +} + +/** + * @brief Compute a single point of the Kaiser window + * This is less efficient than calculating the whole window at once. + * + * @param b Kaiser parameter beta + * @param x Point to evaluate, normalized in 0 to 1 + */ +inline double kaiserWindowSinglePoint(double b, double x) +{ + double t = x + x - 1.0; + return i0(b * std::sqrt(1.0 - t * t)) / i0(b); +} + +/** + * @brief Compute the cardinal sine + */ +template +T sinc(T x) +{ + return (x == T(0)) ? T(1) : (std::sin(x) / x); +} + +/** + * @brief Compute the normalized cardinal sine + */ +template +T normalizedSinc(T x) +{ + return sinc(pi() * x); +} + /** * @brief Finds the minimum size of 2 spans *