Add the kaiser window helper
This commit is contained in:
Jean Pierre Cimalando
parent
ea0c122e0f
commit
75c169217d
6 changed files with 456 additions and 1 deletions
|
|
@ -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"
|
||||
|
|
|
|||
119
external/cephes/LICENSE.txt
vendored
Normal file
119
external/cephes/LICENSE.txt
vendored
Normal file
|
|
@ -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: <steve@moshier.net>
|
||||
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 <julien@cornebise.com>; 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 <julien@cornebise.com>; 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 <steve@moshier.net>
|
||||
X-X-Sender: steve@d510
|
||||
To: Julien Cornebise <julien@cornebise.com>
|
||||
Subject: Re: Cephes: permission to wrap+distribute for Lua
|
||||
In-Reply-To: <52653AD3.1010004@cornebise.com>
|
||||
Message-ID: <alpine.DEB.2.02.1310250827040.17646@d510>
|
||||
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 <COPYRIGHT HOLDER> 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.
|
||||
82
external/cephes/src/chbevl.c
vendored
Normal file
82
external/cephes/src/chbevl.c
vendored
Normal file
|
|
@ -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) );
|
||||
}
|
||||
193
external/cephes/src/i0.c
vendored
Normal file
193
external/cephes/src/i0.c
vendored
Normal file
|
|
@ -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 <math.h>
|
||||
|
||||
/* 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) );
|
||||
|
||||
}
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -25,6 +25,13 @@
|
|||
#include <xmmintrin.h>
|
||||
#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 <class T>
|
||||
constexpr T max(T op1, T op2)
|
||||
{
|
||||
|
|
@ -465,6 +472,54 @@ bool isReasonableAudio(absl::Span<Type> span)
|
|||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Compute the Kaiser window
|
||||
*
|
||||
* @param b Kaiser parameter beta
|
||||
* @param window Span of real which receives the window
|
||||
*/
|
||||
template <class T>
|
||||
void kaiserWindow(double b, absl::Span<T> window)
|
||||
{
|
||||
double i0b = i0(b);
|
||||
for (size_t i = 0, n = window.size(); i < n; ++i) {
|
||||
double x = i / static_cast<double>(n - 1);
|
||||
double t = x + x - 1.0;
|
||||
window[i] = static_cast<T>(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 <class T>
|
||||
T sinc(T x)
|
||||
{
|
||||
return (x == T(0)) ? T(1) : (std::sin(x) / x);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Compute the normalized cardinal sine
|
||||
*/
|
||||
template <class T>
|
||||
T normalizedSinc(T x)
|
||||
{
|
||||
return sinc(pi<T>() * x);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Finds the minimum size of 2 spans
|
||||
*
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue