Beat Frei
Institute for Computer Music and Sound
Technology (ICST)
Baslerstrasse 30,
CH-8048 Zurich, Switzerland
beat.frei@zhdk.ch, http://www.icst.net
Download the recent version from
http://www.icst.net/research/projects/dsp-library
The ICST
DSP Library is a package of C++ routines with focus on rapid development of audio
and measurement applications. It can be seen as a speed-optimized compact
technical computing tool with industrial grade audio algorithms that allows one
to develop, analyze and test math-centered real-time applications by directly
creating the final code without any intermediate steps. While being optimized
for modern x86 processors under Windows and Mac OS X, the library may be useful
on other platforms as well. Its scope is processing, so the author decided, with
regard to the large variety of hosts and protocols, to omit any real-time
interface functionality and leave this part to dedicated frameworks and
libraries.
·
Core
functions based on hand-optimized SSE2 code for maximum speed.
·
Efficient
single-precision floating point arithmetic whenever appropriate without
sacrificing signal quality required for professional audio.
·
Fast
block processing.
·
Zero
overhead support of streaming blocks with varying sizes.
·
Thread-safe.
·
Multiple
platform support.
·
Small
memory footprint.
·
No
dependencies on external libraries.
·
No
license restrictions. Uses only genuine and public domain code.
·
Patent-free
to the best of the author’s knowledge.
·
Make
the library folder visible to your C++ project or copy it into the project
folder.
·
Add
all .cpp and .h files to the project.
·
Include
the desired headers, whereby “Common.h” must always be included. For the full range
of functions and easy access to the library namespace without scope operators,
copy the following lines to the top of your implementation files:
#include "Common.h" // always include this
#include "AudioAnalysis.h"
#include "AudioFile.h"
#include "AudioSynth.h"
#include "BlkDsp.h"
#include "Chart.h"
#include "Neuro.h"
#include "SpecMath.h"
using
namespace icstdsp;
A Windows platform
is detected by the existence of the _WIN32 macro. Otherwise, a POSIX platform
is assumed. Furthermore, it is assumed that including a C++ standard library header
of type <cAnything> puts its functions not only into the std but also the
global namespace, which appears to be valid with popular compilers.
Side
effects:
·
The
namespace “icstdsp” is created.
·
The
macro ICSTDSP_SSEALIGN is defined.
·
On
the Windows platform, NOMINMAX is defined if it does not already exist.
·
The
following standard headers are automatically included by including “Common.h”:
o
<cmath>
o
<cfloat>
o
<cstdlib>
o
<cstdio>
o
<emmintrin.h>
(unless ICSTLIB_NO_SSEOPT is defined)
o
<afxwin.h>
(only if ICSTLIB_ENABLE_MFC is defined)
The library
provides several options to customize it for compatibility and speed. With the
default settings, the target processor must support SSE2 and may be in any
rounding mode. Add the following preprocessor definitions to the compiler’s
project settings as required:
ICSTLIB_NO_SSEOPT
Mandatory if either the target
processor is not guaranteed to support SSE2 or your compiler can’t translate
SSE intrinsics. No SSE code is generated with the effect that many core
functions run slower, typically 2 to 5 times.
ICSTLIB_DEF_ROUND
If you can ensure that the floating-point rounding mode will be in the
default state (round to nearest) whenever a library function is called. Maximizes
the speed of conversions that round floating-point numbers to integers, the “fexp”
function and functions that depend on these.
ICSTLIB_ENABLE_MFC
To indicate
that your C++ development environment supports Microsoft Foundation Classes.
This makes the Chart class available.
ICSTLIB_USE_IPP
Tell the
library to use an installed Intel Performance Primitives package in order to
speed up certain transforms. Consider this option if
your application spends over 20 percent of its time doing transforms. For
maximum performance, first initialize the IPP CPU dispatcher according to version
and linkage (ippInit or ippStaticInit) and then call PrepareTransforms. Remember to
call UnPrepareTransforms upon
termination of the application. Be aware that the IPP package is neither
supplied with this library nor BSD compatible.
Many SSE-optimized
functions internally exist in two different flavours, a faster one with aligned
and a slower one with unaligned memory access, of which the appropriate is dynamically
selected. Although this process is entirely transparent, you are responsible
for providing aligned data as needed to achieve optimum speed. Make yourself
familiar with the aligned memory support of the library (see Auxiliary) and refer to the remark section of
a function to learn more about its specific requirements.
Calculations
involving denormal numbers take a very long time on some popular processors. An
IEEE 754 single precision floating-point number, which is commonly used to
represent the C++ “float” data type, is denormal if it is non-zero and has an
absolute value below ≈ 10-38. If denormals occur frequently, the
CPU load skyrockets and time-constraints of real-time applications are easily
violated. While there may be a flag to tell a particular CPU or one of its
execution units to bash denormals to zero, a host program that calls a plug-in
written with help of this library may still decide to set it the other way and
even worse, other libraries may fail as they rely on denormals being processed
as such. This issue has plagued signal processing on general-purpose CPUs for a
long time and still no satisfactory solution is in sight. Thus, there are many
different concepts around, all tradeoffs between ease of use and efficiency. Here
is the one adopted for this library:
·
Definition:
Safe numbers are either zero or very unlikely (p < 0.001) to have an
absolute value below 10-18.
·
Input
vectors that consist of safe numbers produce a negligible amount of denormals
with any library function except for the rare case of evaluating a polynomial
that has more than two consecutive zero coefficients.
·
It’s
up to you not to produce bunches of unsafe numbers. Guidelines:
o
Choose
-1 to 1 as the standard interval for calculations.
o
Functions
that use an anti-denormal scheme produce small but safe numbers. If this output
is not an immediate input to another function that uses an anti-denormal scheme,
check whether successive operations might result in unsafe numbers and if so, call
“prune” with a limit as high as your
application allows to set small numbers to zero.
·
No
precautions are taken for functions that always produce safe numbers when the input
switches between values that produce safe numbers.
·
If
a function could produce a sequence of denormals internally or unsafe numbers at
the output when the input switches between values that produce safe numbers in
the steady-state, an anti-denormal scheme is used as described in the remark
section of the function.
·
The
library does neither depend on specific CPU settings nor does it change any
flags associated with the handling of denormals.
License: BSD (2-clause
type a.k.a. FreeBSD)
Language: C++
Platforms: Windows (tested with VS 2008, VC 6
and ICC 10)
MacOS (tested with GCC 4.0 under Xcode)
Linux (not tested yet)
Processors: x86, x64
Namespace: icstdsp
Size of download: ≈ 2.9 MB
Source code
lines: ≈ 15000
(excl. Ooura’s math routines)
|
Files |
Classes |
Description |
Code in files depends
on classes… |
|
BlkDsp.cpp BlkDsp.h |
BlkDsp CircBuffer |
Standard block signal processing |
SpecMath |
|
SpecMath.cpp SpecMath.h SpecMathInline.h |
SpecMath |
Specialized math functions |
|
|
AudioAnalysis.cpp AudioAnalysis.h |
AudioAnalysis |
Audio analysis |
BlkDsp SpecMath |
|
AudioSynth.cpp AudioSynth.h |
WaveOsc, RingMod, Noise, Delay, Amp, VarDelay,
SampleOsc, Envelope, MoogFilter, ChambFilter, FMOsc, VAOsc, AudioToPM, Lowpass1, Highpass1, Hilbert |
Audio synthesis |
BlkDsp SpecMath |
|
Neuro.cpp Neuro.h |
Neurons Axons Kohonen PCA |
Neural networks, principal and independent
component analysis |
BlkDsp SpecMath |
|
ffftoourad.cpp ffftoouraf.cpp ffftooura.h |
|
FFT routines by Takuya Ooura |
|
|
AudioFile.cpp AufioFile.h |
AudioFile |
Load and store audio data |
|
|
Chart.cpp Chart.h |
Chart |
2D visualization, WINDOWS ONLY |
|
|
MathDefs.h |
|
Internal math constants and macros |
|
|
CritSect.h |
|
Internals for thread safety |
|
|
Common.h |
|
Standard includes and data types for
interfacing, may be used to create a precompiled header file |
|
Copyright ©
2008, 2009, 2010,
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
THIS
SOFTWARE IS PROVIDED BY THE
This
version 1.2 of the library has beta status. Each function has been tested in
its standard implementation (ICSTLIB_NO_SSEOPT
defined) on Windows with VC 6. Afterwards, the library has been adapted for
POSIX compatibility and SSE implementations were developed on Mac OS X with Xcode
using GCC 4.0. Finally, they have been cross-checked on Windows with VS 2008.
The library has also been compiled with ICC 10 under Windows.
Subject to
minor changes:
·
Documentation
·
Code
for multi-platform compatibility (Linux compatibility assumed but still
untested)
·
Function
implementation (e.g. fine tuning of audio synthesis filters)
FUNCTIONAL REFERENCE
“Size” always denotes the number of elements to
process. Sizes must be at least 1 unless otherwise noted.
The vast majority of non-zero elements should
have absolute values within 10-18…1018. Larger values may
cause internal overflow, smaller non-zero values may produce denormal numbers
that take a very long time to process. Refer to the section about denormals for details.
Align frequently used arrays for optimum SSE
support. See Auxiliary.
Time series are represented as 1D-arrays with
the last element being the most recent.
Column (row) vectors are represented as 1D-arrays
starting with the top (leftmost) element.
Matrices with m rows and n columns are
represented as 1D-arrays as follows:
a[m*n] = {a11,..,a1n, a21,..,a2n,…..,
am1,..,amn}
A single complex number is of type “cpx”
defined in “SpecMath.h”.
Complex data series with “size” elements have
the following format:
d[2*size] = {re0, im0, re1, im1,..,
resize-1, imsize-1}
Polynomials have, unless otherwise noted, the
form sum(i = 0..degree){c[i]*xi} with real single or double
precision coefficients c[i].
Parameters denoted as continuation data are
used to save internal states for block-wise processing of contiguous data. They
must be initialized when the function is called to process the first block and left
unaltered for subsequent blocks.
Input only data is never altered by a function
call.
Functions
that don’t require internal states to be saved except for processing a subsequent
block are collected in memberless classes of static methods - otherwise, they
are member functions of ordinary objects, often with externally supplied work
data.
Examples:
// Class: BlkDsp
// Declaration: static void linear(float*
d, int size, float start, float end)
float d[2048]; //
work data
BlkDsp::linear(d, 2048, 0, 1.0f); //
static method call, class must be specified
// Class: Noise
// Declaration: void
Update(float* out, int samples)
float out[100]; //
work data, not part of Noise object
Noise mynoise(); // create Noise object
mynoise.Update(out, 100); // member function call,
object must be specified
Elementary Real Array Operations
Elementary Complex Array Operations
Decompose to Arbitrary Functions
Audio Synthesis Filters + Amplifiers
Kohonen Feature Maps + Vector
Quantization
Principal + Independent Component
Analysis
Fast
forward (fft) and inverse (ifft) Fourier transform of complex data
Fast
Fourier transform of real data to lower half spectrum
Fast
inverse Fourier transform of lower half spectrum to real data
Fast
Fourier transform of symmetrical lower half real data to lower half spectrum
Fast inverse Fourier transform
of lower half spectrum to symmetrical lower half real data
Fast
forward (dct) and inverse (idct) type 2 cosine transform
Fast
forward (dst) and inverse (idst) type 2 sine transform
Get
k-th bin of the Fourier transform of arbitrary size real data
Forward
Haar wavelet transform
Inverse
Haar wavelet transform
Create
asymptotic exponential segment
Create
logarithmically spaced data points
Create
counter-clockwise rotating complex phasor
Create
sine wave with linearly changing frequency
Create
sine wave with exponentially changing frequency
Create
uniformly distributed noise in the interval [-1,1]
Create
exponentially distributed noise with variance 1
Create
standard Cauchy distributed noise
Create
Gaussian noise with variance 1
Create
Hann (a.k.a. raised cosine) window
Create
3-term Blackman-Harris window with -67 dB side lobe
Create
4-term Blackman-Harris window with -92 dB side lobe
Elementary
Real Array Operations
Get
index of element with largest difference to reference value
Get
index of element with smallest difference to reference value
Fast
natural logarithm of absolute value
Fast
cross correlation (FFT-based)
Fast
biased autocorrelation (FFT-based)
Elementary Complex Array
Operations
Evaluate
polynomial with complex argument
Inverse
of tabulated monotonically increasing function
Sort for descending order while
keeping track of correspondences
Create
subset while keeping track of correspondences
Map
interval linearly to another
Linear
interpolated cyclic table lookup
Elementary
Real Matrix Operations
Subtract
one matrix from another
Multiply
matrix by column vector
Multiply
transpose of matrix by column vector
Multiply
inverse of matrix by another matrix
Time-varying
first order section
Frequency response of
discrete-time LTI system
Frequency
response of continuous-time LTI system
Group
delay of discrete-time LTI system
Group
delay of continuous-time LTI system
Arithmetic
mean of a distribution
Unbiased
variance of a distribution
CDF of the normal distribution
CDF
complement of the normal distribution
Confidence interval of normally
distributed data
Save data to file as raw series
of floats
Load
data from file as raw series of floats
Get
number of elements in file
Return
closest power of two greater or equal to i
Allocate array aligned for SSE
Declare variable or array aligned
for SSE (Macro)
Prepare transforms for maximum
performance
Reset
pointers and fill buffer with zeros
Read
at an offset from the current read pointer without advancing it
Get
maximum number of elements that can be read without overflow
Get maximum number of elements
that can be written without overflow
Quick
and dirty hyperbolic tangent
Quick
and dirty medium range exponential
Natural
logarithm of the gamma function
Bessel
function of the first kind
Fast
rounding double to int conversion
Fast
truncating double to int conversion
Fast truncating split into integer
and fractional part
Find complex roots of polynomial
with real coefficients
Construct polynomial from complex
roots
Symbolic
multiplication by constant
Symbolic
multiplication of two polynomials
Convert
Chebyshev series to polynomial
Convert
polynomial to Chebyshev series
Fit parabola and return extremum
Find
coefficients of Chebyshev approximation
Solve using forth order
Runge-Kutta method
Design
discrete time first order filter
Design
discrete time audio equalizer filter
Convert
continuous to discrete time transfer function
Design
continuous time first order filter
Design
continuous time second order filter
Get
normalized Butterworth low pass filter parameters
Get
normalized Bessel low pass filter parameters
Get
normalized Chebyshev type 1 low pass filter parameters
Convert
normalized to general low pass filter
Convert
normalized low pass to general high pass filter
Convert normalized low pass to
general band pass filter
High
resolution spectral analysis
Prepare
parameters for high resolution spectral analysis
Analysis
by Decomposition to Arbitrary Functions
Convert magnitude spectrum to real
cepstrum
Convert
real cepstrum to magnitude spectrum
Get
Mel frequency cepstral coefficients (MFCC)
Calculate LPC residual from
original signal
Restore original signal from LPC
residual using a direct form IIR filter
Restore original signal from LPC
residual using a lattice IIR filter
Calculate line spectral
frequencies (LSFs) from even order LPC coefficients
Restore original signal from even
order LPC residual using LSFs with an IIR filter
Fundamental
frequency detector
Verify
fundamental frequency against high resolution spectrum
Extract
harmonic spectrum and inharmonicity
Spectral
flatness within a defined frequency band
Amplitude-based
attack transient detector
Spectrum-based
general transient detector
Partial
tracker after McAulay/Quatieri
Add
time-varying cosine wave to existing data matching frequencies and phases
Add time-varying cosine wave to
existing data matching frequencies
Modulator
audio and control update
Carrier
audio and control update
Audio
and control update (linear version)
Audio
and control update (virtual analog version)
Virtual Analog Moog Low Pass
Filter
Convert Audio to Phase Modulation
Signal
Create
axons with randomized weights
Disable
batch learning and update weights
Update
weights by generalized Hebbian algorithm
Create
Kohonen, vector quantization or k-means map
Reset
network to state after construction
Get
output node corresponding to input vector
Decompose
vector into sum of components
Compose
vector as sum of components
Convert
principal to independent components
Load
raw file and specify format
Get
safe pointer to audio data
Plot
vertical bulk interval diagram
Plot
horizontal bulk interval diagram
Get
x value from mouse message
Get
y value from mouse message
Get
right button activity from mouse message
Get
left button activity from mouse message
Get
from mouse message: Has button just been pressed?
Get from mouse message: Has
button just been released?
Get
safe pointer to graphic device context
Standard
block signal processing.
// class: BlkDsp
// input: d[2*size]
= {re0, im0, re1, im1,.., resize-1,
imsize-1}
// output: d[2*size]
= {re0, im0, re1, im1,.., resize-1,
imsize-1}
// remarks: size must be a power of two greater than
1
// common speedup factor when
using an IPP package: > 4 (float), > 2 (double)
static void
fft(float* d, int size)
static void
fft(double* d, int size)
static void
ifft(float* d, int size)
static void
ifft(double* d, int size)
// class: BlkDsp
// input: d[size] = {re0, re1,..,
resize-1}
// output: d[size] = {re0, resize/2, re1, im1,.., resize/2-1,
imsize/2-1}
// remarks: size must be a power of two greater than
1
// common speedup factor when
using an IPP package: > 4 (float), > 2 (double)
static void
realfft(float* d, int size)
static void
realfft(double* d, int size)
// class: BlkDsp
// input: d[size] = {re0, resize/2, re1, im1,.., resize/2-1,
imsize/2-1}
// output: d[size] = {re0, re1,..,
resize-1}
// remarks: size must be a power of two greater than
1
// common speedup factor when
using an IPP package: > 4 (float), > 2 (double)
static void
realifft(float* d, int size)
static void
realifft(double* d, int size)
// class: BlkDsp
// input: d[size+1] = {re0, re1,..,
resize}
// output: d[size+1] = {re0, re1,..,
resize}
// remarks: size must be a power of two greater than
1
static void
realsymfft(float* d, int size)
static void
realsymfft(double* d, int size)
// class: BlkDsp
// input: d[size+1] = {re0, re1,..,
resize}
// output: d[size+1] = {re0, re1,..,
resize}
// remarks: size must be a power of two greater than
1
static void
realsymifft(float* d, int size)
static void
realsymifft(double* d, int size)
// class: BlkDsp
// input: d[size] = {re0, re1,..,
resize-1}
// output: d[size] = {re0, re1,..,
resize-1}
// remarks: size must be a power of two greater than
1
// idct = type 3 cosine
transform
// common speedup factor when
using an IPP package: > 2 (float), > 1 (double)
static void
dct(float* d, int size)
static void
dct(double* d, int size)
static void
idct(float* d, int size)
static void
idct(double* d, int size)
// class: BlkDsp
// input: d[size] = {re0, re1,..,
resize-1}
// output: d[size] = {re0, re1,..,
resize-1}
// remarks: size must be a power of two greater than
1
static void
dst(float* d, int size)
static void
dst(double* d, int size)
static void
idst(float* d, int size)
static void
idst(double* d, int size)
// class: BlkDsp
// input: d[size] = {re0, re1,..,
resize-1}
// output: return value = {rek, imk}
static cpx
goertzel(float* d, int size, int k)
// class: BlkDsp
// input: d[size] = {re0, re1,..,
resize-1}
// output: d[size] = {{re00},{re10},{re20,
re21},{re30, re31, re32, re33},..,{reN0,..,reNsize/2-1}}
// remarks: size must be the N-th power of two and greater
than 1
static void
hwt(float* d, int size)
// class: BlkDsp
// input: d[size] = {{re00},{re10},{re20,
re21},{re30, re31, re32, re33},..,{reN0,..,reNsize/2-1 }}
// output: d[size] = {re0, re1,..,
resize-1}
// remarks: size must be the N-th power of two and
greater than 1
static void
ihwt(float* d, int size)
All windows are symmetric: d[0] = d[size-1]
All signals include the end point.
Periodic versions are obtained by creating
vectors of size + 1 and discarding the last element.
// class: BlkDsp
// input: start = value of d[0]
// end =
value of d[size-1]
// output: d[size]
// remarks: much faster with d aligned
static void
linear(float* d, int size, float start,
float end)
// class: BlkDsp
// input: start = value of d[0]
// end =
value of d[size-1]
// time = time constant as a fraction of size
// output: d[size]
static void
exponential(float* d, int size, float start, float end, float time)
// class: BlkDsp
// input: start = value of d[0]
// end =
value of d[size-1]
// output: d[size]
static void
logspace(float* d, int size, float
start, float end)
// class: BlkDsp
// input: periods = number of periods
// phase = phase shift as fraction of a period
// center = true: set
phase relative to center
// false: set phase relative to beginning
// output: d[size]
// return value = size
> 1: instantaneous phase of d[size-1]
// size
= 1: instantaneous phase of hypothetical d[1]
static float
sine(float* d, int size, float periods, float phase = 0, bool center = false)
// class: BlkDsp
// input: periods = number of
periods
// phase = phase
shift as fraction of a period
// center = true: set phase relative to center
// false: set phase relative to beginning
// output: d[2*size]
= {re0, im0, re1, im1,..,
resize-1, imsize-1}
// return value = size
> 1: instantaneous phase of d[size-1]
// size
= 1: instantaneous phase of hypothetical d[1]
static float
cpxphasor(float* d, int size, float periods, float phase = 0, bool center =
false)
// class: BlkDsp
// input: startpd = number of
periods at d[0]
// endpd = number
of periods at d[size-1]
// phase = phase
shift as fraction of a period
// output: d[size]
// return value = size
> 1: instantaneous phase of d[size-1]
// size
= 1: instantaneous phase of hypothetical d[1]
static float
chirp(float* d, int size, float startpd, float endpd, float phase = 0)
// class: BlkDsp
// input: startpd = number of
periods at d[0]
// endpd = number
of periods at d[size-1]
// phase = phase
shift as fraction of a period
// output: d[size]
// return value = size
> 1: instantaneous phase of d[size-1]
// size
= 1: instantaneous phase of hypothetical d[1]
// remarks: absolute error < 1.1e-6
// the error is a continuous
function that adds low-order harmonics but no noise
static float
expchirp(float* d, int size, float startpd, float endpd, float phase = 0)
// class: BlkDsp
// input: periods = number of
periods
// symmetry = shape: 0 |\|\, 0.5 /\/\, 1 /|/|
// phase = phase
shift as fraction of a period
// center = true: set phase relative to center
// false: set phase relative to beginning
// output: d[size]
// return value = size
> 1: instantaneous phase of d[size-1]
// size
= 1: instantaneous phase of hypothetical d[1]
// remarks: phase = 0 => d[phase
reference point] = 0
static float
saw( float* d, int size, float
periods,
float symmetry = 1.0f, float phase =
0, bool center = false )
// class: BlkDsp
// output: d[size]
// remarks: SSE version uses MWC algorithm with an
approximate period of 2^63
// non-SSE version uses
32-bit congruential generator
static void
unoise(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
enoise(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
cnoise(float* d, int size)
// class: BlkDsp
// input: apxorder = order of
approximation
// output: d[size]
// remarks: set apxorder = 2 to create triangular
noise as used for dithering
static void
gnoise(float* d, int size, int apxorder = 10)
// class: BlkDsp
// input: periods = number of
periods
// output: d[size]
static void
sinc(float* d, int size, double periods)
// class: BlkDsp
// output: d[size]
static void
triangle(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
hann(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
hamming(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
blackman(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
bhw3(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
bhw4(float* d, int size)
// class: BlkDsp
// output: d[size]
static void
flattop(float* d, int size)
// class: BlkDsp
// input: sigma = width
parameter
// output: d[size]
static void
gauss(float* d, int size, double sigma)
// class: BlkDsp
// input: alpha = width
parameter
// output: d[size]
static void
kaiser(float* d, int size, double alpha)
// class: BlkDsp
// input: d[size]
// output: return value = sum of elements
// remarks: faster with d aligned
static float
sum(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = <d,d>
// remarks: faster with d aligned
static float
energy(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = <d,d>/size
// remarks: faster with d aligned
static float
power(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = sqrt(<d,d>/size)
// remarks: faster with d aligned
static float
rms(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = sqrt(<d,d>)
// remarks: faster with d aligned
static float
norm(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = index of most
positive element of d
static int
maxi(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = index of most
negative element of d
static int
mini(float* d, int size)
// class: BlkDsp
// input: d[size]
// r = reference
value
// output: return value = index of element
with largest |d-r|
static int
farthesti(float* d, float r, int size)
// class: BlkDsp
// input: d[size]
// r = reference
value
// output: return value = index of element
with smallest |d-r|
// remarks: “mtabinv” is much faster for
arrays that contain a monotonic function
static int
nearesti(float* d, float r, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// return value = scale
factor
// remarks: if the signal has no energy, it remains
unchanged and 1.0f is returned
// faster with d aligned
static float
normalize(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: faster with d aligned
static void
abs(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: output = 1 for input ≥ 0, output =
-1 for input ≤ 0, output is never 0
// faster with d aligned
static void
sgn(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
static void
finv(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
static void
fsqrt(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: absolute error < 1e-6 for |d| <
2pi
// the error consists of both
a continuous function due to the approximation and a
// noise term that grows linearly
with |d| due to finite argument precision
// slightly faster with d
aligned
static void
fcos(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: absolute error < 1.1e-6 for |d| <
2pi
// the error consists of both
a continuous function due to the approximation and a
// noise term that grows
linearly with |d| due to finite argument precision
// slightly slower than
“fcos”, slightly faster with d aligned
static void
fsin(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: absolute approximation error < 4e-7
// the approximation error is
a continuous function that adds no noise, it does not
// include errors caused by finite
argument and result precision
// slightly faster with d
aligned
static void
logabs(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: relative approximation error < 3e-7
// the approximation error is
a continuous function that adds no noise, it does not
// include errors caused by
finite argument and result precision
// slightly faster with d
aligned
static void
fexp(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
static void
reverse(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// designed for single block
processing
// for streaming blocks, use a
filter with transfer function H(z)
= z/(z-1), but keep
// an eye on the output as it
may grow excessively
static void
cumsum(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: d[size]
// remarks: normalized so that a linear input going
from 0 to 1 results in an output of 1
// designed for single block
processing
// for streaming blocks, a first order high pass filter may be
used, the simplest of
// which has the transfer
function H(z) = (z-1)/T with T = sampling interval
static void
diff(float* d, int size)
// class: BlkDsp
// input: d[size]
// lf = true: extra precision for low frequency signals
// dprev = d[-1], required if lf = true
// dnext = d[size],
required if lf = true
// output: d[size]
// remarks: normalized so that an input of 1 results
in a linear output going from 0 to 1
// designed for single block
processing
// for streaming blocks, a filter with transfer function H(z) =
0.5T(z+1)/(z-c)
// may serve as a leaky integrator,
where c is slightly less than 1 and T denotes
// the sampling interval
static void
integrate(float* d, int size, bool lf = false, float dprev = 0, float dnext = 0)
// class: BlkDsp
// input: c = reference value
// output: d[size] = all elements
set to the reference value
// remarks: much faster with d aligned
static void
set(float* d, float c, int size)
// class: BlkDsp
// input: d[size]
// lim = limit
// rep = value
replaced
// output: d[size] = for all
elements:
// if |d[n]| <= lim then
sign(d[n])*rep -> d[n]
// else
leave d[n] unchanged
// remarks: much faster with d aligned
static void
prune(float* d, float lim, float rep, int size)
// class: BlkDsp
// input: d[size]
// hi = upper
limit
// lo = lower
limit
// output: d[size]
// remarks: faster with d aligned
static void
limit(float* d, int size, float hi, float lo)
// class: BlkDsp
// input: r[size]
// output: d[size]
// remarks: arrays r and d may overlap
static void
copy(float* d, float* r, int size)
// class: BlkDsp
// input: d[size], r[size]
// output: r[size], d[size]
// remarks: faster with both d and r aligned
static void
swap(float* d, float* r, int size)
// class: BlkDsp
// input: d[size], r[size]
// output: d[size]
// remarks: faster with both d and r aligned
static void
max(float* d, float* r, int size)
// class: BlkDsp
// input: d[size], r[size]
// output: d[size]
// remarks: faster with both d and r aligned
static void
min(float* d, float* r, int size)
// .. array
r or constant c element-wise to d
// class: BlkDsp
// input: d[size]
// r[size], c
// output: d[size]
// remarks: faster with both d and r aligned
static void
add(float* d, float c, int size)
static void
add(float* d, float* r, int size)
// .. array
r or constant c element-wise from d
// class: BlkDsp
// input: d[size]
// r[size], c
// output: d[size]
// remarks: faster with both d and r aligned
static void
sub(float* d, float c, int size)
static void
sub(float* d, float* r, int size)
// .. array
r or constant c element-wise by d
// class: BlkDsp
// input: d[size]
// r[size], c
// output: d[size]
// remarks: faster with both d and r aligned
static void
mul(float* d, float c, int size)
static void
mul(float* d, float* r, int size)
// add c*r
element-wise to d
// class: BlkDsp
// input: d[size]
// r[size]
// c
// output: d[size]
// remarks: faster with both d and r aligned
static void
mac(float* d, float* r, float c, int size)
// class: BlkDsp
// input: d[size]
// r[size]
// output: return value = dot product
<d,r>
// remarks: faster with both d and r aligned
static float
dotp(float* d, float* r, int size)
// class: BlkDsp
// input: d[size]
// r[size]
// output: return value = squared
Euclidean distance
// remarks: faster with both d and r aligned
static
float sdist(float* d, float* r, int size)
// class: BlkDsp
// input: d[dsize]
// r[rsize]
// output: d[dsize + rsize - 1]
// remarks: faster with r aligned and dsize >
rsize, alignment of d
irrelevant
static void
conv(float* d, float* r, int dsize, int rsize)
// class: BlkDsp
// input: d[dsize]
// r[rsize]
// output: d[N] = d[0…dsize + rsize - 2]: convolution of d and r
// d[dsize
+ rsize – 1…N-1]: undefined
// r[N] = undefined
// remarks: N = next higher power of two larger or equal to dsize +
rsize, see nexthipow2
// slightly faster with both
d and r aligned
static void
fconv(float* d, float* r, int dsize, int rsize)
// class: BlkDsp
// input: d[dsize]
// r[rsize]
// output: d[pts] = cross correlation d x r
// remarks: pts ≤ dsize
// faster with both d and r
aligned
static void
ccorr(float* d, float* r, int dsize, int rsize, int pts)
// class: BlkDsp
// input: d[dsize]
// r[rsize]
// output: d[N] = d[0…dsize - 1]: cross correlation d x r
// d[dsize…N-1]: undefined
// r[N] = undefined
// remarks: N = next higher power of two larger or equal to dsize +
rsize, see nexthipow2
// slightly faster with both
d and r aligned
static void
fccorr(float* d, float* r, int dsize, int rsize)
// class: BlkDsp
// input: r[rsize]
// output: d[dsize]
// remarks: dsize ≤ rsize
// faster with r aligned, alignment of d irrelevant
static void
uacorr(float* d, float* r, int dsize,
int rsize)
// class: BlkDsp
// input: r[rsize]
// output: d[dsize]
// remarks: dsize ≤ rsize
// faster with r aligned, alignment of d irrelevant
static void
bacorr(float* d, float* r, int dsize, int rsize)
// class: BlkDsp
// input: d[size]
// output: d[N] = d[0…size - 1]: autocorrelation of d
// d[size…N-1]: undefined
// remarks: N = next higher power of two larger or equal to 2*size,
see nexthipow2
static void
facorr(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: return value = sum of elements
// remarks: faster with d aligned
static cpx
cpxsum(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: return value = <d,d*>
// remarks: faster with d aligned
static float
cpxenergy(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: return value = <d,d*>/size
// remarks: faster with d aligned
static float
cpxpower(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: return value = sqrt(<d,d*>/size)
// remarks: faster with d aligned
static float
cpxrms(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: return value = sqrt(<d,d*>)
// remarks: faster with d aligned
static float
cpxnorm(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: d[2*size]
// return value = scale
factor
// remarks: if the signal has no energy, it remains
unchanged and 1.0f is returned
// faster with d aligned
static float
cpxnormalize(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: d[2*size]
// remarks: faster with d aligned
static void
cpxconj(float* d, int size)
// class: BlkDsp
// input: d[2*size]
// fullrange = true: |d| < 2.95e-39 yields signed infinity
// false: much faster, |d| < 1.09e-19 yields signed
infinity
// output: d[2*size]
// remarks: divide by zero yields signed infinity
// slightly faster with d
aligned
static void
cpxinv(float* d, int size, bool fullrange=false)
// class: BlkDsp
// input: r[2*size]
// output: d[size]
// remarks: d = r is ok
// slightly faster with both d
and r aligned
static void
cpxmag(float* d, float* r, int size)
// class: BlkDsp
// input: r[2*size]
// output: d[size]
// remarks: d = r is ok
// faster with both d and r
aligned
static void
cpxpow(float* d, float* r, int size)
// class: BlkDsp
// input: r[2*size]
// output: d[size]
// remarks: d = r is ok
// equivalent to d[i] =
atan2(r[2*i + 1], r[2*i])
static void
cpxarg(float* d, float* r, int size)
// class: BlkDsp
// input: d[2*size]
// lim = limit
// rep = value
replaced
// output: d[2*size] = for all
elements (n = 0,2,4,..):
// if
|d[n]| ≤ lim then d[n] = rep, d[n+1] = 0
// else
leave d[n] and d[n+1] unchanged
// remarks: slightly faster with d aligned
static void
cpxprune(float* d, float lim, float rep, int size)
// class: BlkDsp
// input: c = reference value
// output: d[2*size] = set all
elements to the reference value
// remarks: much faster with d aligned
static void
cpxset(float* d, cpx c, int size)
// class: BlkDsp
// input: r[2*size]
// output: d[2*size]
// remarks: arrays r and d may overlap
static void
cpxcopy(float* d, float* r, int size)
// class: BlkDsp
// input: d[2*size], r[2*size]
// output: r[2*size], d[2*size]
// remarks: faster with both d and r aligned
static void
cpxswap(float* d, float* r, int size)
// .. array
r or constant c element-wise to d
// class: BlkDsp
// input: d[2*size]
// r[2*size], c
// output: d[2*size]
// remarks: faster with both d and r aligned
static void
cpxadd(float* d, cpx c, int size)
static void
cpxadd(float* d, float* r, int size)
// .. array
r or constant c element-wise from d
// class: BlkDsp
// input: d[2*size]
// r[2*size], c
// output: d[2*size]
// remarks: faster with both d and r aligned
static void
cpxsub(float* d, cpx c, int size)
static void
cpxsub(float* d, float* r, int size)
// .. array
r or constant c element-wise by d
// class: BlkDsp
// input: d[2*size]
// r[2*size], c
// output: d[2*size]
// remarks: faster with both d and r aligned
static void
cpxmul(float* d, cpx c, int size)
static void
cpxmul(float* d, float* r, int size)
// add c*r
element-wise to d
// class: BlkDsp
// input: d[2*size]
// r[2*size]
// c
// output: d[2*size]
// remarks: faster with both d and r aligned
static void
cpxmac(float* d, float* r, cpx c, int size)
// class: BlkDsp
// input: d[2*size]
// r[2*size]
// output: return value = dot product
<d,r*>
// remarks: faster with both d and r aligned
static cpx
cpxdotp(float* d, float* r, int size)
// class: BlkDsp
// input: d[2*size]
// output: re[size]
// remarks: re = d is ok
static void
cpxre(float* re, float* d, int size)
// class: BlkDsp
// input: d[2*size]
// output: im[size]
// remarks: im = d is ok
static void
cpxim(float* im, float* d, int size)
// class: BlkDsp
// input: re[size] = real data
(can be NULL)
// im[size] = imaginary
data (can be NULL)
// output: d [2*size]
// remarks: if re = NULL then Re{d} = 0, if im =
NULL then Im{d} = 0
static void
realtocpx(float* d, float* re, float* im, int size)
// class: BlkDsp
// input: mag[size] = magnitude
// arg[size] = argument
in radians
// output: d[2*size]
// remarks: relative error < 1.1e-6 for |arg|
< 2pi
// the error consists of both
a continuous function due to the approximation and a
// noise term that grows
linearly with |arg| due to finite argument precision
// slightly faster with d
aligned
static void
cpxptc(float* d, float* mag, float* arg, int size)
// class: BlkDsp
// input: d[2*size]
// mag[size] = new
magnitude
// output: d[2*size]
// remarks: an input (0, 0) results in an output (0,
0)
// slightly faster with both d
and mag aligned
static void
cpxcombine(float* d, float* mag, int size)
// class: BlkDsp
// input: r[size]
// output: d[size]
static void
ftod(double* d, float* r, int size)
// class: BlkDsp
// input: r[size]
// output: d[size]
static void
dtof(float* d, double* r, int size)
// class: BlkDsp
// input: r[size]
// output: d[size]
// remarks: clips output to [-231, 231-1]
// always rounds to nearest
independent of the actual FPU/SSE rounding mode
// faster with both d and r
aligned
static void
ftoi(int* d, float* r, int size)
// class: BlkDsp
// input: r[size]
// output: d[size]
// faster with both d and r
aligned
static void
itof(float* d, int* r, int size)
// r[n]
-> d[n*interval + offset], n = 0..rsize-1
// class: BlkDsp
// input: r[rsize]
// output: d[dsize]
// remarks: dsize > (rsize-1)*interval + offset ≥ 0
// offset ≥ 0
// skipped elements of d remain
unchanged
static void
interleave(float* d, float* r, int rsize, int interval, int offset)
// r[n*interval
+ offset] -> d[n], n = 0..dsize-1
// class: BlkDsp
// input: r[rsize]
// output: d[dsize]
// remarks: rsize > (dsize-1)*interval + offset ≥ 0
// offset ≥ 0
// d = r is ok if offset ≥
(dsize – 1)*(1 – interval)
static void
deinterleave(float* d, float* r, int dsize, int interval, int offset)
// sum(i = 0..order){c[i]*d[n]^i}->
d[n], n = 0..dsize-1
// class: BlkDsp
// input: d[dsize]
// c[order + 1] = real
coefficients of polynomial
// order = order
of polynomial
// output: d[dsize]
// very fast all-float
version, double version is much slower
static void
polyval(float* d, int dsize, float* c, int order)
static void
polyval(float* d, int dsize, double* c, int order)
// calculate
element-wise for complex numbers d: sum(i = 0..order){c[i]*d^i}-> d
// class: BlkDsp
// input: d[2*dsize]
// c[order + 1] = real
coefficients of polynomial
// order = order
of polynomial
// output: d[2*dsize]
static void
cpxpolyval(float* d, int dsize, float* c, int order)
static void
cpxpolyval(float* d, int dsize, double* c, int order)
// sum(i =
0..order){c[i]*Ti(d[n])} -> d[n], n = 0..dsize-1
// Ti(x)
= i-th order Chebyshev polynomial of the first kind
// class: BlkDsp
// input: d[dsize]
// c[order + 1] = coefficients
of series
// order = order
of series
// output: d[dsize]
// very fast all-float
version, double version is much slower
static void
chebyval(float* d, int dsize, float* c, int order)
static void
chebyval(float* d, int dsize, double* c, int order)
// class: BlkDsp
// input: d[size]
// output: idx[size] = idx[0..number
of peaks-1]: indices of peaks
// idx[number
of peaks..size-1]: undefined
// example:
d[idx[0]] is the value of the first peak
// return value = number
of peaks
// remarks: the first and last elements of d are
never classified as peaks
static int
findpeaks(float* d, int* idx, int size)
// class: BlkDsp
// input: d[size]
// output: idx[size] = idx[0..number
of dips-1]: indices of dips
// idx[number
of dips..size-1]: undefined
// example:
d[idx[0]] is the value of the first dip
// return value = number
of dips
// remarks: the first and last elements of d are
never classified as dips
static int
finddips(float* d, int* idx, int size)
// class: BlkDsp
// input: d[size] = tabulated
function, must be
monotonically increasing
// x = function
value
// output: return value = i of d[i]
closest to but not greater than x
// special
case: 0 if d[0] > x
// remarks: performs a bisection
// use “nearesti” for non-monotonic
functions
static int
mtabinv(float* d, float x, int size)
// class: BlkDsp
// input: d[size]
// output: d[size] = data
sorted so that d[n] ≥ d[n + 1]
// idx[size] = index
at which this element was in the original array
// remarks: uses quicksort algorithm
static void
isort(float* d, int* idx, int size)
// select a
single or two complementary subsets from the full data set
// class: BlkDsp
// input: d[size] = full
set
// sel[size] = selectors
// output: a[size] = a[0..asize-1]: subset with all d[i] with sel[i]
> 0
// a[asize..size-1]: undefined
// aidx[size] = aidx[0..asize-1]: full set indices of subset a elements
// aidx[asize..size-1]: undefined
// example: a[3] = d[aidx[3]]
// b[size] = b[0..bsize-1]: subset with all d[i] with sel[i] ≤ 0
// b[bsize..size-1]: undefined
// bidx[size] = bidx[0..bsize-1]: full set indices of subset b elements
// bidx[bsize..size-1]: undefined
// return value = asize
(note: bsize = size – asize)
// remarks: either a or b = d is ok
// output options: a = NULL: compute only aidx and bidx
// aidx
= NULL: compute
only a and b
// additional option: b = NULL: compute neither b nor bidx
static int
select(float* d, float* sel, float* a, float* b, int* aidx, int* bidx, int
size)
// class: BlkDsp
// input: d[size] = phase
data in radians wrapping around within [-pi, pi]
// output: d[size] = continuous
phase data in radians
static void
unwrap(float* d, int size)
// class: BlkDsp
// input: r[size] = original
data
// rhi = bound
1 of original data
// rlo = bound
2 of original data
// dhi = bound
1 is mapped to this value
// dlo = bound
2 is mapped to this value
// output: d[size] = linearly
mapped data
// remarks: no limiting, out-of-bound input results
in correctly mapped out-of-bound output
// faster with d aligned
static void
maplin(float* d, float* r, int size, float dhi, float dlo, float rhi, float rlo)
// the output
linearly fades from d to r with increasing index
// class: BlkDsp
// input: d[size]
// r[size]
// output: d[size]
static void
xfade(float* d, float* r, int size)
// (1 - w[n])*d[n]
+ w[n]*r[n] -> d[n], n = 0..size-1
// class: BlkDsp
// input: d[size]
// r[size]
// w[size]
// output: d[size]
static void
xfade(float* d, float* r, float* w, int size)
// the
output fades from d to r with increasing index while maintaining constant power
for
// uncorrelated
equal power input signals
// class: BlkDsp
// input: d[size]
// r[size]
// output: d[size]
static void
pxfade(float* d, float* r, int size)
// this
function creates a periodic signal from a single cycle given as a table of N
elements
// with N a
power of two
// class: BlkDsp
// input: t[N + 1] = lookup table
// fill
t[0.. N-1] with single cycle and set t[N] = t[0]
// ltsize = base-2 logarithm of N, hence N = 2^ltsize
// start = start position
// step = index increment
// output: d[dsize] = generated
signal
// return value = start
position to generate the next output value
// remarks: N > start ≥ 0, step ≥ 0
static float
linlookup(float* d, int dsize, float* t, int ltsize, float start = 0, float
step = 1.0f)
// class: BlkDsp
// input: a[n*n]
// output: return value = determinant of a
static float
mdet(float* a, int n)
// class: BlkDsp
// input: a[n*n]
// output: return value = trace of a
static float
mtrace(float* a, int n)
// class: BlkDsp
// input: a[m*n]
// output: a[n*m]
static void
mxpose(float* a, int m, int n)
// class: BlkDsp
// output: a[n*n]
static void
mident(float* a, int n)
// class: BlkDsp
// input: a[m*n]
// c
// output: a[m*n]
// remarks: faster with a aligned
static void
mmul(float* a, float c, int m, int n)
// class: BlkDsp
// input: a[m*n]
// b[m*n]
// output: a[m*n]
// remarks: faster with both a and b aligned
static void
madd(float* a, float* b, int m, int n)
// class: BlkDsp
// input: a[m*n]
// b[m*n]
// output: a[m*n]
// remarks: faster with both a and b aligned
static void
msub(float* a, float* b, int m, int n)
// class: BlkDsp
// input: a[m*n]
// r[n]
// output: d[m]
// remarks: faster with both a and r aligned and n a
multiple of 4, alignment
of d irrelevant
static void
mmulv(float* d, float* a, float* r, int m, int n)
// class: BlkDsp
// input: a[m*n]
// r[m]
// output: d[n]
// remarks: about half as fast as “mmulv”
static void
mtmulv(float* d, float* a, float* r, int m, int n)
// class: BlkDsp
// input: a[m*n]
// b[n*p]
// output: c[m*p]
// remarks: uses “mmulv”, not particularly fast for
small matrices
static void
mmulm(float* c, float* a, float* b, int m, int n, int p)
// class: BlkDsp
// input: a[m*n]
// b[m*m]
// output: a[m*n] = resulting matrix
// b[m*m] = undefined
// return value = determinant
of matrix b
// remarks: if the return value is zero, matrix a
remains unchanged
static float
minvmulm(float* a, float* b, int m, int n)
// class: BlkDsp
// input: r[size] = original
signal
// n = delay
in sampling intervals
// c[n], cp = continuation
data, init with 0
// output: d[size] = original
signal delayed by n sampling intervals
// c[n], cp = continuation
data
static void
delay(float* d, float* r, int size, float* c, int& cp, int n)
// transfer
function: H(z) = sum(i
= 0..order){b[i]*z^-i}
// class: BlkDsp
// input: d[size] = original
signal
// b[order + 1] = filter
coefficients
// c[order] = continuation
data, init with 0
// output: d[size] = filtered
signal
// c[order] = continuation
data
static void
fir(float* d, int size, float* b, int order, float* c)
static void
fir(float* d, int size, double* b, int order, float* c)
// transfer
function: H(z) = 1
/ (1 + sum(i = 1..order){a[i]*z^-i})
// class: BlkDsp
// input: d[size] = original
signal
// a[order + 1] = filter
coefficients
// c[order] = continuation
data, init with 0
// output: d[size] = filtered
signal
// c[order] = continuation
data
// remarks: coefficient a[0] is ignored
// anti-denormal scheme sets output elements with
absolute value < 1e-15 to zero
static void
iir(float* d, int size, double* a, int order, double* c)
// transfer
function: H(z) = (b[0]
+ b[1]z^-1) / (1 + a[1]z^-1)
// class: BlkDsp
// input: d[size] = original
signal
// a[2], b[2] = filter
coefficients
// c[1] = continuation data, init with 0
// output: d[size] = filtered
signal
// c[1] = continuation data
// remarks: coefficient a[0] is ignored
// anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
static void
iir1(float* d, int size, float* a, float* b, float* c)
// transfer
function: H(z) = (bs[0]..be[0]
+ bs[1]..be[1]*z^-1) / (1 + as[1]..ae[1]*z^-1)
//
coefficients evolve linearly from start to end values with increasing index of
d
// class: BlkDsp
// input: d[size] = original
signal
// as[2], bs[2] = start
filter coefficients
// ae[2], be[2] = end
filter coefficients
// c[1] = continuation data, init with 0
// output: d[size] = filtered
signal
// c[1] = continuation data
// remarks: coefficients as[0] and ae[0] are ignored
// anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
static void
viir1(float* d, int size, float* as, float* ae, float* bs, float* be, float* c)
// transfer
function: H(z) = (b[0]
+ b[1]z^-1 + b[2]z^-2) / (1 + a[1]z^-1 + a[2]z^-2)
// class: BlkDsp
// input: d[size] = original
signal
// a[3], b[3] = filter
coefficients
// c[2] = continuation data, init with 0
// output: d[size] = filtered
signal
// c[2] = continuation data
// remarks: coefficient a[0] is ignored
// anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
static void
biquad(float* d, int size, double* a, double* b, double* c)
// class: BlkDsp
// input: d[size] = original
signal
// c[2] = continuation data, init with 0
// output: d[size] = filtered
signal
// c[2] = continuation data
static void
med3(float* d, int size, float* c)
// class: BlkDsp
// input: d[size] = original
signal
// c[4] = continuation data, init with 0
// output: d[size] = filtered
signal
// c[4] = continuation data
static void
med5(float* d, int size, float* c)
// transfer
function: H(z) = sum(i
= 0..border){b[i]z^-i} / sum(i = 0..aorder){a[i]z^-i}
// class: BlkDsp
// input: f[size] = frequencies as a fraction of the
sample rate at which the
// response
is calculated
// a[aorder + 1] = transfer
function polynomial coefficients (denominator)
// b[border + 1] = transfer
function polynomial coefficients (numerator)
// aorder = order
of denominator polynomial
// border = order
of numerator polynomial
// output: mag[size] = magnitude
response, linear absolute value
// phase[size] = phase
response, radians wrapped to [-pi, pi]
static void
freqz( float* mag, float* phase,
float* f, int size,
double* a, int aorder, double* b,
int border )
// transfer
function: G(s) = sum(i
= 0..border){b[i]s^i} / sum(i = 0..aorder){a[i]s^i}
// class: BlkDsp
// input: f[size] = frequencies
in Hz at which the response is calculated
// a[aorder + 1] = transfer
function polynomial coefficients (denominator)
// b[border + 1] = transfer
function polynomial coefficients (numerator)
// aorder = order
of denominator polynomial
// border = order
of numerator polynomial
// output: mag[size] = magnitude
response, linear absolute value
// phase[size] = phase
response, radians wrapped to [-pi, pi]
static void
freqs( float* mag, float* phase,
float* f, int size,
float* a, int aorder, float* b, int
border )
// transfer
function: H(z) = sum(i
= 0..border){b[i]z^-i} / sum(i = 0..aorder){a[i]z^-i}
// class: BlkDsp
// input: f[size] = frequencies
as a fraction of the sample rate at which the
// group
delay is calculated
// a[aorder + 1] = transfer
function polynomial coefficients (denominator)
// b[border + 1] = transfer
function polynomial coefficients (numerator)
// aorder = order
of denominator polynomial
// border = order
of numerator polynomial
// output: gdelay[size] = group delay in
sampling intervals
static void
gdelz(float* gdelay, float* f, int size, double* a, int aorder, double* b, int
border)
// transfer
function: G(s) = sum(i
= 0..border){b[i]s^i} / sum(i = 0..aorder){a[i]s^i}
// class: BlkDsp
// input: f[size] = frequencies
in Hz at which the response is calculated
// a[aorder + 1] = transfer
function polynomial coefficients (denominator)
// b[border + 1] = transfer
function polynomial coefficients (numerator)
// aorder = order
of denominator polynomial
// border = order
of numerator polynomial
// output: gdelay[size] = group delay in
seconds
static void
gdels(float* gdelay, float* f, int size, float* a, int aorder, float* b, int
border)
// class: BlkDsp
// input: d[size]
// output: return value = arithmetic mean
of d
// remarks: faster with d aligned
static float
mean(float* d, int size)
// operates
on a number of bins specified by a value and its probability of occurrence
// class: BlkDsp
// input: d[size] = probability
(need not be normalized)
// r[size] = value
// output: return value = arithmetic mean
of r weighted by d
// remarks: much faster with both d and r aligned
// equivalent to calculating
the barycenter
static float
pdmean(float* d, float* r, int size)
// class: BlkDsp
// input: d[size]
// output: return value = unbiased
variance of d
// remarks: faster with d aligned
static float
var(float* d, int size)
// operates
on a number of bins specified by a value and its probability of occurrence
// class: BlkDsp
// input: d[size] = probability
(need not be normalized)
// r[size] = value
// output: return value = unbiased
variance of r weighted by d
static float
pdvar(float* d, float* r, int size)
// class: BlkDsp
// input: d[size]
// m = measure
of central tendency, e.g. mean or median of d
// output: return value = average absolute
deviation of d
static float
aad(float* d, int size, float m)
// class: BlkDsp
// input: d[size]
// output: return value = median of d
static float
median(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = median absolute
deviation of d
// remarks: very robust to outliers but also much
slower than aad with mean
static
float mad(float* d, int size)
// class: BlkDsp
// input: d[size] = data
// p = p-value
// output: return value = p-quantile
// remarks: uses Hazen’s rule
static float
quantile(float* d, float p, int size)
// class: BlkDsp
// input: d[size]
// output: return value = geometric mean
of d
static float
gmean(float* d, int size)
// class: BlkDsp
// input: d[size]
// output: return value = harmonic mean
of d
static float
hmean(float* d, int size)
// class: BlkDsp
// input: x[size], y[size]
// output: return value = covariance of x
and y
// remarks: faster with both x and y aligned
static float
cov(float* x, float* y, int size)
// class: BlkDsp
// input: x[size], y[size]
// output: return value = Pearson product-moment
correlation of x and y
// remarks: faster with both x and y aligned
static float
cpears(float* x, float* y, int size)
// class: BlkDsp
// input: x[size], y[size]
// output: return value = Spearman’s rank
correlation of x and y
static float
cspear(float* x, float* y, int size)
// class: BlkDsp
// input: x[size], y[size]
// output: return value =
// remarks: consider cspear as a faster alternative for
large size
static float
ckend(float* x, float* y, int size)
// calculate
weights c so that data points y[j] are estimated in a least-square sense as a linear
// combination
of data points x1..n[j] and a constant: yest[j] = c[0] +
sum(i=1..n){c[i]*xi[j]};
// equivalent
to the least-squares approximation of data vector y as a linear combination of
// basis vectors
x1..n and the vector (1,..,1): yest = c[0]*(1,..,1) + sum(i=1..n){c[i]*xi}.
// class: BlkDsp
// input: x[n*size] = data points
of regressors, basis vectors
// format:
{x1[size], .., xn[size]}
// y[size] = data
points of dependent variable, data vector
// n = number
of regressors or basis vectors, n < size
// cov[N] = if
not NULL and cov[0] > 0: use c covariance
matrix
// if
not NULL and cov[0] ≤ 0: save c covariance
matrix
// N
= (n+1)*(n+1)
// output: c[n+1] = regression coefficients, weights of basis
vectors
// y[size] = residual:
y - yest
// return value = re: coefficient of determination, fraction
of y
// variance
explained by linear model, if y has zero
// mean:
1 – residual/original y energy
// im: model efficiency based on the Akaike criterion,
// a
higher value indicates a better model
// p[2] = if not NULL: p[0] =
p-value of hypothesis c[1..n] = 0
// p[1]
= f-value of the underlying f-statistic
// with
degrees of freedom (n, size-n-1)
// sdc[n+1] = if
not NULL: standard deviation of c[0..n]
// note: c[i] errors follow a
t(size-n-1)-distribution under
// quite
general assumptions, so approximate c[i]
// confidence
intervals are easily found through civn
// cov[N] = if
not NULL: saved c covariance matrix, cov[0] > 0
// remarks: do one-way ANOVA by encoding each factor
as a regressor and watching p[0]
// repeated regressions with
the same regressors are much faster if the covariance
// matrix of the regression
coefficients is reused by setting cov[0] = 0 for the first
// call and passing cov unchanged
for subsequent calls
// consider “matching pursuit” if n is fixed and lower than
the number of potential
// basis vectors
static cpx
linreg( float* c, float* x,
float* y, int size, int n,
float* p=NULL, float* sdc=NULL,
float* cov=NULL )
// class: BlkDsp
// input: d = degree of polynomial, d < size
// x[size] = x-value
of value pair
// y[size] = y-value
of value pair
// w[size] = weight
of value pair (w = NULL: no
weighting)
// output: c[d+1] = coefficients of weighted least squares
polynomial:
// yest(x)
= sum(i = 0..d){c[i]*xi}
// y[size] = residual:
y - yest
// remarks: use “chebyapprox” to convert a
tabulated function to a high-order polynomial
static void
polyfit(float* c, int d, int size, float* x, float* y, float* w=NULL)
// create
histogram by classifying data into bins of equal width
// class: BlkDsp
// input: d[dsize] = data
// bins = number of bins
// dmin = minimum value to be classified, goes to bin[0]
// dmax = maximum
value to be classified, goes to bin[bins-1]
// output: bin[bins] = number of
data elements in bin
// return value = number
of values outside of [dmin, dmax]
static int histogram(float* d, int* bin, int dsize, int bins, float dmin, float dmax)
// class: BlkDsp
// input: x[size] = x data
// y[size] = y
data
// norm = false: use x data for x quantiles
// = true: use
standard normal distribution for x quantiles
// output: x[size] = x
quantiles in ascending order
// y[size] = y
quantiles in ascending order
static void
qqpairs(float* x, float* y, int size, bool norm=false)
// class: BlkDsp
// input: d[size] = data
// th = detection
threshold relative to general purpose setting
// output: score[size] = distance/critical
distance - 1.0f, range: -1.0f..inf
// >
0: suspected outlier
// return value = number
of outliers
// remarks: insensitive to the masking effect due to
median-based algorithm
// intended to be tuned with “th” to
application-specific distributions and sizes
// use sign to dichotomize scores if required
static int
outliers(float* score, float* d, int size, float th=1.0f)
// test the
hypothesis that there are no outliers in a normally distributed data set
// class: BlkDsp
// input: d[size] = data
// mv[2] = unused:
NULL
// faster:
{arithmetic mean of d, unbiased variance of d}
// output: return value = p-value
// idx[1] = if
not NULL: idx[0] = index i of farthest
d[i]
static
float grubbs(float* d, int size, int* idx=NULL, float* mv=NULL)
// test the
hypothesis that two means differ by a predetermined value
// class: BlkDsp
// input: ref = reference value, usually 0
// xsize = number of x data points
// xmean = arithmetic
mean of x, s. mean
// xvar = unbiased variance of x, s. var
// onsided = false: test for mean(x) – mean(y) = ref
// true: test for mean(x) – mean(y) <= ref
// ysize = number of y data points
// >
0: two-sample test of unpaired data
// ≤
0: one-sample or paired two-sample test
// ymean = arithmetic
mean of y
// yvar = unbiased variance of y
// output: return value = re: p-value
// im: Cohen’s d
// t[1] = if not NULL: t[0] =
t-value of the underlying t-statistic
// with
degrees of freedom xsize - 1 (one-sample and paired
// test)
or xsize + ysize – 2 (unpaired two-sample test)
// remarks: a two-sample test of paired data (x,y) is
equivalent to a one-sample test of x – y
static cpx
ttest( float ref, int xsize, float
xmean, float xvar, bool onesided=false,
int ysize=0, float ymean=0, float
yvar=0, float* t=NULL )
// test the
hypothesis that two normally distributed data sets have equal variance
// class: BlkDsp
// input: xvar = unbiased variance of x, s. var
// xsize = number of x data points, < 1e6
// yvar = unbiased variance of y
// ysize = y > 0: number of
y data points, < 1e6
// y
≤ 0: use predetermined exact y
variance
// onsided = false: test for var(x) = var(y)
// true: test for var(x) <= var(y)
// output: return value = p-value
// remarks: use the Levene test if you suspect that x
and y are not normally distributed or
// sizes are too big
static
float ftest(float xvar, int xsize, float yvar, int ysize, bool onesided=false)
// test the
hypothesis that two arbitrary distributed data sets have equal variance
// class: BlkDsp
// input: x[xsize] = x data
// y[ysize] = y
data
// onsided = false: test for var(x) = var(y)
// true: test for var(x) <= var(y)
// bf = false: Levene test, for symmetrical distributions
// true: Brown-Forsythe test, for skewed
distributions
// output: return value = p-value
// remarks: use the F test
if there is clear evidence that x and y are normally distributed
static
float levtest(float* x, int xsize, float* y, int ysize, bool onesided=false,
bool bf=false)
// test the
hypothesis that an observed distribution equals a theoretical one
// class: BlkDsp
// input: fo[size] = observed frequencies
(numbers of values per bin)
// fe[size] = expected
frequencies or bin probabilities
// all elements must be > 0
// size = number of bins, < 1e6
// typical
range: sqrt(N)..N/5, N = number of
samples
// rdf = reduction
of degrees of freedom: degrees = size – rdf
// typical
values: 1 for testing distributions, (columns +
// rows
-1) for testing contingency tables, rdf < size
// yates = false: no effect
// true: enable Yates correction, recommended if you
// have
both less than 100 samples and only one
// degree
of freedom, typically applied to 2*2
// contingency
tables
// output: return value = re: p-value
// im: chi squared
// remarks: as a rule of thumb 80% of the values in both
fo and fe should be at least 5 to
// get a meaningful result
(combining bins may help)
// Fisher’s exact test (not
provided by this library due to realtime constraints) is
// recommended for 2*2 tables
with less than 20 samples
// prefer the binomial test
for dichotomous data
static cpx
chi2test(float* fo, float* fe, int size, int rdf=1, bool yates=false)
// test the
hypothesis H0 that observed data falls in two categories with a
given probability
// class: BlkDsp
// input: pa = predetermined probability of category
A
// asize = number of observed category A elements
// tsize = sample size, < 1e6
// onesided = false: conservative p-value estimate of H0,
// good
if both pa*tsize and (1-pa)*tsize > 5,
// exact
for pa = 0.5
// true: p-value = exact chance of observing at
least
// asize
elements
// output: return value = p-value
// remarks: swap categories if both asize is near
zero and the difference from one matters
static float
binomtest(float pa, int asize, int tsize, bool onesided=false)
// class: BlkDsp
// input: x = upper limit
// mean = arithmetic mean
// std = standard
deviation
// output: return value = Phi(x)
// remarks: use cdfcn for good resolution far
out in the upper tail
static
float cdfn(float x, float mean, float std)
// class: BlkDsp
// input: x = upper limit
// mean = arithmetic mean
// std = standard
deviation
// output: return value = 1 - Phi(x)
// remarks: use cdfn for good resolution far
out in the lower tail
static
float cdfcn(float x, float mean, float std)
// class: BlkDsp
// input: mean = arithmetic mean
// std = standard
deviation
// a = confidence
level
// output: return value = re: lower limit
// im: upper limit
// remarks: data following a t-distribution with degrees
of freedom > 25 leads to a
// 95%-confidence interval that
is just 5% larger than the normal approximation,
// so this function matches a
broad range of applications, e.g. C.I. of the mean
// of N values from a standard
normal distribution ≈ civn(mean, 1.0f/sqrtf(N))
static cpx
civn(float mean, float std, float a=0.95f)
// class: BlkDsp
// input: d[size] = data
// size = number of elements to write
// filename = full
pathname including extension
// append = true: append data to existing file
// false: overwrite existing file
// output: return
value = number of elements written, -1 on failure
static int
save(float* d, int size, char* filename, bool append = false)
// class: BlkDsp
// input: filename = full
pathname including extension
// size = number of elements to read
// offset = first element to read
// output: d[size] = data
// return value = number
of elements read, -1 on failure
static int
load(float* d, int size, char* filename, int offset = 0)
// class: BlkDsp
// input: filename = full
pathname including extension
// output: return value = number of
elements in file, -1 on failure
static int
getsize(char* filename)
// class: BlkDsp
static int
nexthipow2(int i)
// class: BlkDsp
// input: size = number of elements
// output: return value = pointer to
first element
// remarks: if ICSTLIB_NO_SSEOPT
is defined then:
// m128, m128d, m128i versions are not
available
// an unaligned standard “malloc” is
used
static float*
sseallocf(int size)
static double*
sseallocd(int size)
static int*
ssealloci(int size)
static __m128*
sseallocm128(int size)
static __m128d*
sseallocm128d(int size)
static __m128i*
sseallocm128i(int size)
// class: BlkDsp
// input: p = pointer to first element
// remarks: use only with arrays allocated by “ssealloc{f, d, i,
m128, m128d, m128i}”
// uses an
unaligned standard “free” if ICSTLIB_NO_SSEOPT is defined
static void
ssefree(void* p)
// usage: place the macro before the usual declaration
// don’t declare multiple instances separated by
commas
// example: ICSTDSP_SSEALIGN float mydata[1000];
// remarks: not as portable as “ssealloc”
// has no effect if ICSTLIB_NO_SSEOPT is defined
ICSTDSP_SSEALIGN
// class: BlkDsp
// remarks: strongly recommended if an IPP package
is used, no effect otherwise
// call it during the
initialization phase of the application
// call UnPrepareTransforms
before the application terminates
static void
PrepareTransforms()
// class: BlkDsp
// remarks: must be called before the application
terminates if PrepareTransforms
has ever
// been called, has no effect
if PrepareTransforms was never called
static void
UnPrepareTransforms()
Fixed
length circular buffer.
A method call that would cause an overflow
condition has no effect.
// class: CircBuffer
// input: size = buffer size in elements
CircBuffer(int
size = 1)
// class: CircBuffer
void
Reset()
// class: CircBuffer
// input: d[n] = data to write
// n = number
of elements to write
// output: return value = true on
success, false on overflow
bool
Write(float* d, int n)
// class: CircBuffer
// input: n = number of elements to read
// output: d[n] = data read
// return value = true
on success, false on overflow
bool
Read(float* d, int n)
// class: CircBuffer
// input: n = number of elements to read
// offset = start
reading at the current location - offset
// output: d[n] = data read
// return value = true
on success, false on overflow
bool
ReadPassive(float* d, int n, int offset)
// class: CircBuffer
int
GetReadSize()
// class: CircBuffer
int
GetWriteSize()
Signal
processing specific mathematical functions and routines.
// if the
general shape matters more than exact values, suitable as a nonlinearity where
an
// unlimited
input must be limited by a function that is differentiable everywhere within a
// half-plane
// class: SpecMath
// input: x, y (single
version)
// x[4], y[4] (quad
version)
// output: return value = approximate
atan2(y, x) (single
version)
// x[4] = approximate atan2(y[4], x[4]) (quad
version)
// remarks: absolute error < 0.072
static float
qdatan2f(float y, float x)
static void
qdatan2f(float* y, float* x)
// a fair
compromise of general shape, shape of derivative, and good asymptotic
behaviour,
// suitable
as a nonlinearity where an unlimited input must be limited by a function that
is
// differentiable
everywhere
// class: SpecMath
// input: x (single
version)
// x[4] (quad
version)
// output: return value = approximate tanh(x) (single version)
// x[4] = approximate tanh(x[4])
(quad version)
// remarks: absolute error < 0.0034 (SSE
version), 0.029 (basic version)
static float
qdtanh(float x)
static void
qdtanh(float* x)
// class: SpecMath
// input: x (single
version)
// x[4] (quad
version)
// output: return value = approximate
exp(x) (single
version)
// x[4] = approximate exp(x[4])
(quad version)
// remarks: relative error < 1.2e-4 within
recommended input range 0 ≤ x ≤10
static float
qdexp(float x)
static void
qdexp(float* x)
// class: SpecMath
// input: x (single
version)
// x[4] (quad
version)
// output: return value = approximate
exp(x) (single
version)
// x[4] = approximate exp(x[4])
(quad version)
// remarks: relative error < 2.3e-6 within
recommended input range 0 ≤ x ≤10
static float
fexp(float x)
// class: SpecMath
// input: x
// output: return value = arcsinh(x)
static float
asinhf(float x)
static double
asinh(double x)
// class: SpecMath
// input: x
// output: return value = arccosh(x)
static float
acoshf(float x)
static double
acosh(double x)
// class: SpecMath
// input: x
// output: return value = arctanh(x)
static float
atanhf(float x)
static double
atanh(double x)
// class: SpecMath
// input: x
// output: return value = erf(x)
// remarks: consider erfc if you need the precise distance
to the asymptote at the tails
// use probit for the inverse
static float
erff(float x)
static double
erf(double x)
// class: SpecMath
// input: x
// output: return value = erfc(x)
// remarks: use erfc(-x) instead of 2 – erfc(x) for
x << 0 to prevent cancellation
// if high resolution at the
tails is not required, use erf which is about
twice as fast
// use probit for the inverse
static
float erfcf(float x)
static
double ercf(double x)
// inverse cumulative
distribution function of the standard normal distribution
// class: SpecMath
// input: 0 ≤ x ≤ 1
// output: return value = probit(x) = sqrt(2)*erf-1(2x
- 1)
static
float probit(float x)
// class: SpecMath
// input: x > 0
// output: return value = ln(gamma(x))
// remarks: float version relative + absolute error
< 3e-7
static
float gammalnf(float x)
static
double gammaln(double x)
// class: SpecMath
// input: 0 < x < 1e6
// y ≥ 0
// output: return value = on success: P(x,y)
// out
of range: -1.0f
// remarks: this function uses an iteration and is
potentially slow for large x very close to y
static
float rgamma(float x, float y)
// class: SpecMath
// input: 0 ≤ x ≤ 1
// 0 < a < 1e6
// 0 < b < 1e6
// output: return value = on success: I(x,a,b)
// out
of range: -1.0f
// remarks: this function uses an iteration and is
potentially slow if a or b are large and the
// output is close to 0.5
static
float rbeta(float x, float a, float b)
// class: SpecMath
// input: x
// n
// output: return value = Jn(x)
static
float bessj(float x, int n)
// class: SpecMath
// input: -231 ≤ x ≤ 231-1
// output: return value = integer closest
to x
// remarks: always rounds to nearest independent of
the actual FPU/SSE rounding mode
static int
fdtoi(double x)
// class: SpecMath
// input: -231 ≤ x ≤ 231-1
// output: return value = integer between
0 and x that is closest to x
static int
ffdtoi(double x)
// class: SpecMath
// input: 0 ≤ x ≤ 231-1
// output: x = fractional part
// return value = largest
integer ≤ input
static int
fsplit(double &x)
// class: SpecMath
// input: f = pointer to the function
"func" which has the form
// "float
myfunc(float x, float* p, int n)" with p
// independent
of x
// p = parameters
passed to “func”
// n = number of parameters
// xmax = maximum x included in the search
// xmin = minimum x included in the search
// maxerr = maximum
absolute error of root
// 0:
highest precision, >0: faster
// output: return value = root: x for which func(x) = 0
// remarks: xmax and xmin must produce non-zero function values of
opposite sign
// uses bisection to
guarantee that the root is always found within a fixed upper
// limit of function calls
static
float froot( float
(*f)(float,float*,int), float* p, int n, float xmax, float xmin,
float maxerr=0 )
See General Considerations for the representation of
Polynomials.
Format of Chebyshev series:
sum(i = 0..order){c[i]*Ti(x)}with Ti(x)
= i-th order Chebyshev polynomial of the first kind
// this
routine is especially tailored for real-time signal processing applications, it
is not
// intended
to be a general-purpose replacement for standard solvers
// class: SpecMath
// input: c[d + 1] = real coefficients
of polynomial
// d = degree
of polynomial
// output: r[2*d] = complex roots, format: {re0,
im0,..,red-1, imd-1}
// return value = number
of roots or
// -1:
pathologic case, rare in physics-based applications
// remarks: very fast (s.p. version 20 x faster than
Matlab 6.1 roots for x9 + 1, PM, VC6)
// finite runtime guaranteed
by adaptive precision algorithm
// double precision version is
the recommended default
// factors supporting the single
version:
// no double
precision coefficients available
// speed matters more than precision
// roots do not
aggregate around themselves
// few roots with Im(z)/Re(z)
< 0.05
// extensive
testing
static int
roots(float* c, float* r, int d)
static int
roots(double* c, float* r, int d)
// class: SpecMath
// input: r[2*d] = complex
roots, format: {re0, im0,..,red-1, imd-1}
// d = degree
of polynomial
// output: cre[d + 1] = real part of
polynomial coefficients
// cim[d + 1] = imaginary
part of polynomial coefficients
// remarks: the polynomial is normalized so that the
coefficient of the highest-degree term
// is 1, i.e. cre[d] = 1.0f,
cim[d] = 0
static void
roottops(double* cre, double* cim, float* r, int d)
// class: SpecMath
// input: a[da + 1] = coefficients
of polynomial a
// da = degree
of polynomial a
// b[db + 1] = coefficients
of polynomial b
// db = degree
of polynomial b
// output: a[N] = coefficients of polynomial a + b, N = max(da,db)
+ 1
// remarks: ensure that a provides enough space for the result if db
> da
static void
addpoly(float* a, int da, float* b, int db)
static void
addpoly(double* a, int da, double* b, int db)
// class: SpecMath
// input: a[da + 1] = coefficients
of polynomial a
// da = degree
of polynomial a
// b[db + 1] = coefficients
of polynomial b
// db = degree
of polynomial b
// output: a[N] = coefficients of polynomial a - b, N = max(da,db)
+ 1
// remarks: ensure that a provides enough space for the result if db
> da
static void
subpoly(float* a, int da, float* b, int db)
static void
subpoly(double* a, int da, double* b, int db)
// class: SpecMath
// input: c[d + 1] = coefficients
of polynomial c
// d = degree
of polynomial c
// a = constant
// output: c[d + 1] = coefficients
of polynomial a*c
static void
mulpoly(float* c, float a, int d)
static void
mulpoly(double* c, double a, int d)
// class: SpecMath
// input: a[da + 1] = coefficients
of polynomial a
// da = degree
of polynomial a
// b[db + 1] = coefficients
of polynomial b
// db = degree
of polynomial b
// output: a[N] = coefficients of polynomial a*b, N = da
+ db + 1
// remarks: ensure that a provides enough space for the result
static void
pmulpoly(float* a, int da, float* b, int db)
static void
pmulpoly(double* a, int da, double* b, int db)
// class: SpecMath
// perform
the following operation for each coefficient: a^n*c[n] -> c[n]
// input: c[d + 1] = coefficients
of polynomial c
// d = degree
of polynomial c
// a = constant
// output: c[d + 1] = coefficients
of scaled polynomial
static void
scalepoly(float* c, float a, int d)
static void
scalepoly(double* c, double a, int d)
// class: SpecMath
// replace
argument of polynomial a with polynomial b
// input: a[da + 1] = coefficients
of polynomial a
// da = degree
of polynomial a
// b[db + 1] = coefficients
of polynomial b
// db = degree
of polynomial b
// output: a[N] = coefficients of polynomial a(b), N =
da*db + 1
// remarks: ensure that a provides enough space for the result
static void
sspoly(float* a, int da, float* b, int db)
static void sspoly(double* a, int da, double* b, int db)
// class: SpecMath
// input: c[d + 1] = coefficients
of polynomial
// d = degree
of polynomial
// output: c[d + 1] = coefficients
of reversed polynomial
static void
revpoly(float* c, int d)
static void
revpoly(double* c, int d)
// class: SpecMath
// input: c[d + 1] = coefficients
of polynomial
// d = degree
of input polynomial
// output: c[d + 1] = coefficients
of differentiated polynomial, c[d] = 0
static void
diffpoly(float* c, int d)
static void
diffpoly(double* c, int d)
// class: SpecMath
// input: c[d] = coefficients of polynomial
// d = degree
of output
polynomial
// output: c[d + 1] = coefficients
of integrated polynomial
static void
integratepoly(float* c, int d)
static void
integratepoly(double* c, int d)
// class: SpecMath
// input: d = degree of polynomial
// output: c[d + 1] = coefficients
of Chebyshev polynomial
static void
chebypoly(double* c, int d)
// class: SpecMath
// input: c[d + 1] = coefficients
of Chebyshev series
// d = order
of Chebyshev series
// output: c[d + 1] = coefficients
of polynomial
static void
chebytops(double* c, int d)
// class: SpecMath
// input: c[d + 1] = coefficients
of polynomial
// d = order
of Chebyshev series
// output: c[d + 1] = coefficients
of Chebyshev series
static void
pstocheby(double* c, int d)
// class: SpecMath
// input: 0 ≤ x ≤ 1.0f = x
value
// y[2] = y values: y[0] = y(0), y[1] = y(1.0f)
// output: return value = linear interpolated
value y(x)
// remarks: use “fsplit” to convert a
fractional index to x and the base index of y
static float
fitlin(float* y, float x)
// class: SpecMath
// input: ym1 = y(-1.0f)
// y0 = y(0)
// y1 = y(1.0f)
// output: return value = {re, im}: {x, y} of extremum
// {0,
y0} if no extremum in x = [-1, 1]
static cpx
paraext(float ym1, float y0, float y1)
// the
interpolation function goes through (x0, y0) and (x1, y1) with specified gradients
// class: SpecMath
// input: x0, y0 = data
point 1
// m0 = gradient
at point 1
// x1, y1 = data
point 2
// m1 = gradient
at point 2
// output: c[4] = coefficients of interpolation
polynomial:
// yest
= c[0] + c[1]*x + c[2]*x2 + c[3]*x3
// remarks: best performance is obtained in the
interval x = [x0, x1]
static void
fitcubic(double* c, float x0, float y0, float m0, float x1, float y1, float m1)
// the
interpolation function goes through (x0, y0) and (x1, y1) with continuous gradients
// input: x[4] = data
points x values: {xm1, x0, x1, x2}
// y[4] = data points y values: {ym1, y0, y1, y2}
// output: c[4] = coefficients of interpolation
polynomial:
// yest
= c[0] + c[1]*x + c[2]*x2 + c[3]*x3
// remarks: best performance is obtained in the
interval x = [x0, x1]
static void
fitcubic(double* c, float* x, float* y)
// the
interpolation function goes through (0, y0) and (1.0f, y1) with continuous
gradients
// fast
version with fixed data point x values: x[4] = {-1.0f, 0, 1.0f, 2.0f}
// input: y[4] = data points y values: {ym1, y0, y1, y2}
// output: c[4] = coefficients of interpolation
polynomial:
// yest
= c[0] + c[1]*x + c[2]*x2 + c[3]*x3
// remarks: best performance is obtained in the
interval x = [0, 1.0f]
static void
fitcubic(float* c, float* y)
// the interpolation
function goes through (0, y0) and (1.0f, y1) with continuous gradients
// fast
version with fixed data point x values: x[4] = {-1.0f, 0, 1.0f, 2.0f}
// class: SpecMath
// input: y[4] = data points y values: {ym1, y0, y1, y2}
// 0
≤ x ≤ 1.0f = x at which y(x) is interpolated
// output: return value = interpolated value
// remarks: use “fsplit” to convert a
fractional index to x and the base index of y
static float
fitcubic(float* y, float x)
// class: SpecMath
// input: y[size] = y
values taken at constant
intervals from x = -1.0f to1.0f
// size = number of tabulated data points
// output: c[d+1] = coefficients of Chebyshev series
approximation:
// yest(x)
= sum(i = 0..d){c[i]*Ti(x)}
// return value = error
relative to maximum absolute value of y
// remarks: use “maplin” to handle functions
with any x range
// use “polyfit” if y is noisy or taken at irregular
intervals
// you may use “chebytops” to convert the
Chebyshev series to a polynomial
static float
chebyapprox(double* c, int d, float* y, int size)
//
approximate deltay using the classic 4th order Runge-Kutta algorithm so that
// y(x +
deltax) = y(x) + deltay given the equation dy/dx = func(y(x), p[0..n-1])
// with p
independent of x
// class: SpecMath
// input: y = current value of y(x)
// deltax = integration
step
// f = pointer
to the function "func" which has the form
// "float
myfunc(float y, float* p, int n)"
// p[n] = parameters passed to "func"
// n = number
of parameters
// output: return value = deltay
// remarks: handle systems of differential equations
as follows:
// 1. call rk4 for each equation passing all
system variables y0,y1,..
// in
p excluding the one that is the y of the current call,
// save
the deltay without updating the system variables
// 2. update all system variables at once by
adding the deltay
// if "func" has to
be time-dependent, create a new system variable t
// with equation dt/dx = 1.0f
and pass it to "func" in p
// usage: 1. define
function "func" as described above
// 2. set y to initial value
// 3. repeat y += rk4(myfunc,p,n,y,deltax) to get consecutive y
values
static float
rk4(float (*f)(float, float*, int), float* p, int n, float y, float deltax)
General purpose discrete time second order
filters are obtained by designing a continuous time version followed by “stoz”.
Design filters of order > 2 as a series of
first and second order sections.
// find
coefficients of transfer function (b[0] + b[1]z^-1) / (1 + a[1]z^-1)
// class: SpecMath
// input: fc = -3db corner or shelving midpoint
frequency as a
// fraction
of the sample rate
// type = filter type: 0
low pass, 2 high pass, 4 low shelving
// 5
high shelving, 7 allpass
// dbgain = shelving
gain in dB
// output: a[2], b[2] = filter coefficients,
a[0] = 1
static void
dzbilin(float* a, float* b, float fc, int type, float dbgain = 0)
// find
coefficients of transfer function (b[0] + b[1]z^-1 + b[2]z^-2) / (1 + a[1]z^-1
+ a[2]z^-2)
// of an
audio equalizer filter based on EQ cookbook formulae by Robert Bristow-Johnson
// class: SpecMath
// input: fc = -3db corner corner / center / shelving
midpoint frequency
// as
a fraction of the sample rate
// qbw = if
> 0: quality factor, if ≤ 0: -bandwidth in octaves
// type = filter type: 0
low pass, 1 band pass, 2 high pass,
// 3
notch, 4 low shelving, 5 high shelving,
// 6
peaking, 7 allpass
// dbgain = peaking
/ shelving gain in dB
// output: a[3], b[3] = filter
coefficients, a[0] = 1
static void
eqzbiquad(double* a, double* b, float fc, float qbw, int type, float dbgain = 0)
// by
applying the bilinear z-transform
// G(s) = (d[0]
+..+ d[order]s^order) / (c[0] +..+ c[order]s^order)
// H(z) =
(b[0] +...+ b[order]z^-order) / (1 + a[1]z^-1 +...+ a[order]z^-order)
// class: SpecMath
// input: c[order + 1] = coefficients of continuous
time transfer function
// d[order + 1] = coefficients
of continuous time transfer function
// order = order
of transfer function polynomials
// fs = sample
rate in Hz
// output: a[order + 1] = coefficients of discrete
time transfer function, a[0] = 1
// b[order + 1] = coefficients
of discrete time transfer function
static void
stoz(double* a, double* b, float* c, float* d, int order, float fs)
// find
coefficients of transfer function (b[0] + b[1]s) / (a[0] + s)
// class: SpecMath
// input: fc = -3db corner or shelving midpoint
frequency in Hz
// type = filter type: 0
low pass, 2 high pass, 4 low shelving
// 5
high shelving, 7 allpass
// dbgain = shelving
gain in dB
// output: a[2], b[2] = filter
coefficients, a[1] = 1
static void
dsbilin(float* a, float* b, float fc, int type, float dbgain = 0)
// find
coefficients of transfer function (b[0] + b[1]s + b[2]s^2) / (a[0] + a[1]s + s^2)
// class: SpecMath
// input: fc = -3db corner corner / center / shelving
midpoint frequency
// in
Hz
// qbw = if
> 0: quality factor, if ≤ 0: -bandwidth in octaves
// type = filter type: 0
low pass, 1 band pass, 2 high pass,
// 3
notch, 4 low shelving, 5 high shelving,
// 6
peaking, 7 allpass
// dbgain = peaking
/ shelving gain in dB
// output: a[2], b[2] = filter
coefficients, a[2] = 1
static void
dsbiquad(float* a, float* b, float fc, float qbw, int type, float dbgain = 0)
// calculate
corner frequencies and quality factors of a series of first and second order
lowpass
// filters
that realize a Butterworth lowpass filter with corner frequency = 1 Hz
// class: SpecMath
// input: d = filter order
// output: f[N] = corner frequencies
// q[N] = quality factors
// N
= d/2 for d even, (d+1)/2 for d odd
// remarks: all stages are of second order, except
for the last stage of odd order filters,
// which is of first order
with an assigned quality factor of 0
static void
fqbutter(float* f, float* q, int d)
// calculate
corner frequencies and quality factors of a series of first and second order
lowpass
// filters
that realize a Bessel lowpass filter with corner frequency = 1 Hz
// class: SpecMath
// input: d ≤ 10 = filter order
// output: f[N] = corner frequencies
// q[N] = quality factors
// N
= d/2 for d even, (d+1)/2 for d odd
// remarks: all stages are of second order, except
for the last stage of odd order filters,
// which is of first order
with an assigned quality factor of 0
static void
fqbessel(float* f, float* q, int d)
// calculate
corner frequencies and quality factors of a series of first and second order
lowpass
// filters
that realize a Chebyshev type 1 lowpass filter with corner frequency = 1 Hz
// class: SpecMath
// input: d = filter order
// rdb = pass
band ripple in db
// output: f[N] = corner frequencies
// q[N] = quality factors
// N
= d/2 for d even, (d+1)/2 for d odd
// remarks: all stages are of second order, except
for the last stage of odd order filters,
// which is of first order
with an assigned quality factor of 0
static void
fqcheby(float* f, float* q, float rdb, int d)
// class: SpecMath
// input: f[n] = corner frequencies of normalized filter
// n = number
of corner frequencies
// fc = new
corner frequency in Hz
// (fs = sample
rate in Hz, specify to map frequencies
// for
use with the bilinear z-transform, see „stoz“)
// output: f[n] = corner frequencies of modified filter
static void
fqlowpass(float* f, float fc, int n, float fs)
// class: SpecMath
// input: f[n] = corner frequencies of normalized
filter
// n = number
of corner frequencies
// fc = new
corner frequency in Hz
// (fs = sample
rate in Hz, specify to map frequencies
// for
use with the bilinear z-transform, see „stoz“)
// output: f[n] = corner frequencies of modified filter
static void
fqhighpass(float* f, float fc, int n, float fs)
// make a series
of 2n second order bandpass
filters out of n normalized low pass filters
// class: SpecMath
// input: f[n] = corner frequencies of normalized
filter
// q[n] = associated quality factors
// if
q < 0.5 then a single second order filter is created
// n = number
of corner frequencies
// fc = band
pass center frequency in Hz
// bw = absolute
bandwidth in Hz
// (fs = sample
rate in Hz, specify to map frequencies
// for
use with the bilinear z-transform, see „stoz“)
// output: f[2*n] = band pass center frequencies
// 0
indicates “no filter here” (happens for q < 0.5)
// q[2*n] = associated
quality factors
// g[2*n] = associated
gain factors
// remarks: check band pass center frequencies to determine whether a
filter exists
static void
fqbandpass(float* f, float* q, float* g, float fc, float bw, int n, float fs)
Signal
analysis functions with focus on audio.
// class: AudioAnalysis
// input: d[size] = data to
analyze
// size = length of data block to analyze,
// must be a power of two larger
than 1
// w[size + 1] = symmetric
window function, as used with “prespec”
// dw[size + 1] = differentiated
window as obtained from “prespec”
// rw[size + 1] = time-ramped
window as obtained from “prespec”
// temp[size] = just
workspace (initialization of content not required)
// aic[4] = interpolation
coefficients as obtained from “prespec”
// output: d[size] = complex
half spectrum of windowed data,
// format:
{re0, im0, re1, im1,.., resize/2-1,
imsize/2-1}
// freq[size] = freq[0..size/2-1]: reassigned frequency of bin as
// a
fraction of the sample rate
// freq[size/2..size-1]: undefined
// amp[size] = amp[0..size/2-1]: reassigned amplitude of peak bins,
// 0:
bin does not represent a peak
// amp[size/2..size-1] undefined
// time[size] = time[0..size/2-1]: reassigned time position of bin
// relative
to center as a fraction of
// “size”,
range: [-0.5, 0.5]
// time[size/2..size-1]: undefined
// remarks: freq / amp / time can be NULL, in which
case the vector is not calculated
// time reassignment yields
the energy barycenter of the windowed
input, i.e.
// a wide pulse is located well only if it falls
into a flat region of the window
static void
spec( float* d, float* freq, float*
amp, float* time, int size,
float* w, float* dw, float* rw,
float* temp, float* aic )
// class: AudioAnalysis
// input: w[size + 1] = symmetric
window function
// size = length of data block to analyze with “spec”
// output: dw[size + 1] = differentiated
window
// rw[size + 1] = time-ramped
window
// aic[4] = coefficients
for amplitude interpolation
static void
prespec(float* w, float* dw, float* rw, float* aic, int size)
//
approximate data as a weighted sum of atoms from a dictionary:
// d ~ sum(i
= 0..elements-1){weight[i]*atom(idx[i])}
// class: AudioAnalysis
// input: data[size] = data to be decomposed
// dict[atoms*size] = dictionary
of atoms with L2 norm = 1,
// format:
{atom0[size], .., atomatoms-1[size]}
// atoms = number of atoms in the dictionary
// elements = number of summands of the decomposition
// output: weight[elements] = weight
of atom
// idx[elements] = index
of atom
// data[size] = residual: original data – approximated data
// remarks: faster with both data and dict aligned
and size a multiple of 4
static void
matchingpursuit( float* weight, int*
idx, int elements, float* dict,
int atoms, float* data, int size )
// class: AudioAnalysis
// input: d[size + 1] = magnitude
lower half spectrum
// output: d[size + 1] = real lower
half cepstrum
// remarks: size must be a power of two larger than 1
static void
magspectorceps(float* d, int size)
// class: AudioAnalysis
// input: d[size + 1] = real lower
half cepstrum
// output: d[size + 1] = magnitude
lower half spectrum
// remarks: size must be a power of two larger than 1
static void
rcepstomagspec(float* d, int size)
// from
power (mathematically more conclusive) or magnitude spectrum (ETSI ES 201108)
// class: AudioAnalysis
// input: d[size + 1] = lower half power
or magnitude spectrum
// bands = number
of filter bands
// cofs = number of MFC coefficients
// fs = sample
rate in Hz
// fmax = highest analysis frequency in Hz
// fmin = lowest analysis frequency in Hz
// output: c[cofs] = MFC
coefficients
// remarks: typical frame specs: duration = 25 ms,
overlap = 50%, hamming windowed
static void
spectomfcc( float* d, float* c, int size,
int bands = 23, int cofs = 13,
float fs = 44100.0f, float fmax =
8000.0f, float fmin = 64.0f )
// *** Example code ***
// const
int order = 10;
// double
a[order + 1], k[order];
// float
rm[order + 1], c[order];
// BlkDsp::set(c,
0, order);
// BlkDsp::bacorr(rm,
inputdata, order + 1, 1000);
// result =
lpdurbin(a, k, rm, order);
// if
(result.re == 1.0f) { handle
unstable solution, never
// occurs with biased autocorrelation }
// lpanalyze(inputdata,
c, a, 1000, order);
// solve
rm[] * a[1 ..] = [err 0..0] for a[] and err
// class: AudioAnalysis
// input: rm[order + 1] = biased
autocorrelation (see "bacorr")
// order = predictor
order
// output: a[order + 1] = linear
prediction coefficients
// k[order] = reflection
coefficients
// return value = re: maximum absolute reflection coefficient,
// 1
if predictor is unstable
// im: residual to input energy ratio
// remarks: using biased autocorrelation guarantees
a stable predictor, which furthermore
// performs similar to the
burg and covariance method if the input is longer than
// 20*order and hamming or
hann windowed
static cpx
lpdurbin(double* a, double* k, float* rm, int order)
// class: AudioAnalysis
// input: d[size] = original
signal
// a[order + 1] = linear
prediction coefficients
// c[order] = continuation
data, init with zero
// order = predictor
order
// output: d[size] = residual
signal
// c[order] = continuation
data
static void
lpanalyze(float* d, float* c, double* a, int size, int order)
// class: AudioAnalysis
// input: d[size] = residual
signal
// a[order + 1] = linear
prediction coefficients
// c[order] = continuation
data, init with zero
// order = predictor
order
// output: d[size] = reconstructed
original signal
// c[order] = continuation
data
// remarks: anti-denormal scheme sets output elements with absolute
value < 1e-15 to zero
static void
lpdsynth(float* d, double* c, double* a, int size, int order)
// the
second version allows gliding linearly from one k-set to another
// class: AudioAnalysis
// input: d[size] = residual
signal
// k[order] = reflection
coefficients (static version)
// kstart[order] = reflection
coefficients at d[0] (dynamic version)
// kend[order] = reflection
coefficients at d[size-1] (dynamic versio
// c[order] = continuation
data, init with zero
// order = predictor
order
// output: d[size] = reconstructed
original signal
// c[order] = continuation
data
// remarks: anti-denormal scheme sets output elements with absolute
value < 1e-15 to zero
static void
lplsynth(float* d, double* c, double* k, int size, int order)
static void
lplsynth(float* d, double* c, double* kstart, double* kend, int size, int order)
// class: AudioAnalysis
// input: a[order + 1] = linear
prediction coefficients
// order = predictor
order, must be even
// grid = number of bins to look for frequencies, default: 1024
// must be a power of two larger
than 4*order + 4
// output: f[order] = ascending
line spectral frequencies as a fraction of the
// sample
rate
// return value = true: success
// false: frequencies not separable on given grid
// remarks: frequencies can be separated if they are
at least 2*fs/grid apart
static bool
lpctolsf(float* f, double* a, int order, int grid)
// the
second version allows gliding linearly from one f-set to another
// class: AudioAnalysis
// input: d[size] = residual
signal
// order = predictor
order, must be even
// f[order] = line
spectral frequencies (static version)
// fstart[order] = line
spectral frequencies at d[0] (dynamic version)
// fend[order] = line
spectral frequencies at d[size-1] (dynamic version)
// c[2*order +1] = continuation
data, init with zero
// output: d[size] = reconstructed
original signal
// c[2*order +1] = continuation
data
// remarks: anti-denormal scheme sets output elements with absolute
value < 1e-15 to zero
static void
lpssynth(float* d, float* f, double* c, int size, int order)
static void
lpssynth(float* d, float* fstart, float* fend, double* c, int size, int order)
// adaptively
estimate y based on x minimizing the error e = yest(x) - y
// class: AudioAnalysis
// input: x[size] = source
data
// y[size] = reference
data
// order = filter
order
// mu = adaptation
rate: 0 (no adaption) …1.0f (fastest stable)
// h[order + 1] = current
filter coefficients, init
with zero for first call
// c[order] = continuation
data, init with zero
// output: e[size] = error
signal
// h[order + 1] = updated
filter coefficients
// c[order] = continuation
data
static void
nlms( float* e, float* x, float*
y, float* h, float* c,
int size, int order, float mu = 1.0f )
// based on
a Costas loop, frequencies and bandwidths given as a fraction of the sample
rate
// class: AudioAnalysis
// input: d[size] = modulated
signal
// fstart = start frequency
// pstart = start
phase, range: [0, 2π]
// tbw = tracking
bandwidth, determines how quickly the carrier
// oscillator
follows frequency changes of the input signal
// dbw = demodulation
bandwidth, determines how quickly the
// demodulated
signal tracks phase changes of the input,
// dbw
is internally forced to be at least 2*tbw
// c[2] = continuation data, init with zero
// output: r[2*size] = demodulated
quadrature signal,
// format:{i0,
q0,…, isize-1, qsize-1}
// fstart = updated start frequency
// pstart = updated
start phase
// c[2] = continuation data
// remarks: anti-denormal scheme enforces minimum output amplitude of
≈ 1e-15
static void
costas( float* d, float* r, int
size, float tbw, float dbw,
float &fstart, float
&pstart, float* c )
// time
constants specified in sampling intervals
// class: AudioAnalysis
// input: d[size] = audio
signal
// atime = attack time constant
// rtime = release time constant
// type = detector: 0 absolute value, 1 power, 2 RMS
// c = continuation
data, init with 0
// output: r[size] = envelope
signal
// c = continuation
data
// remarks: anti-denormal scheme imposes a lower limit of 1e-15 on
output values
static void
envelope( float* d, float* r, float
&c, int size,
float atime = 200.0f, float rtime =
2000.0f, int type = 0)
// based on
normalized autocorrelation
// class: AudioAnalysis
// input: d[size] = audio
signal
// type = normalization scheme:
// 0 McLeod/Wyvill
(intended for musical applications)
// 1 Cauchy-Schwarz
inequality (mathematical F0)
// 2 biased C-S inequality (tunable
compromise)
// 3 YIN (originally for
speech)
// output: return value = re: fundamental
frequency as fraction of the sample rate
// im:
tonality measure 0 (untuned)...1.0f (purely harmonic)
// remarks: precision of tonality measure does not depend
on a correct fundamental
// THE INFORMATION BELOW IS
SUPPLIED WITHOUT LIABILITY:
// type 3 has been subject to
a patent application (WO02097793) that appears to
// be withdrawn in
// the author learns that nevertheless
there is an active patent, he will remove the
// algorithm and make type 3
fall back to type 2
static cpx
fundamental(float* d, int size, int type = 2)
// class: AudioAnalysis
// input: amp[size] = reassigned amplitudes
obtained from "spec"
// fo = fundamental
frequency as fraction of the sample rate
// output: return value = index pointing to spectral bin of
verified fundamental
// frequency
static int
fundverify(float* amp, int size, float fo)
// class: AudioAnalysis
// input: hms = number of harmonics to extract
// amp[size] = reassigned
amplitudes obtained from "spec"
// freq[size] = reassigned
frequencies obtained from "spec"
// fo = fundamental
frequency as fraction of the sample rate
// output: idx[hms] = index
pointing to bin of harmonic in amp[] and freq[],
// example: freq[idx[0]] = frequency of
fundamental
// return value = inharmonicity
J that approximates the n-th partial
// frequency as fn =
n*fo*(1 + J*(n^2 - 1))
static float
harmonics(int* idx, int hms, float* amp, float* freq, int size, float fo)
// class: AudioAnalysis
// input: d[size] = audio
signal
// c = continuation
data, init with zero
// output: return value = estimated number
of zero crossings per sample,
// 0:
no estimate made
// c = continuation
data
static float
zerocross(float* d, float &c, int size)
// class: AudioAnalysis
// input: d[size] = power
lower half spectrum
// f1, f2 = band limit frequencies as fractions of the sample rate
// output: return value = spectral
flatness: 0 (pure oscillations) .. 1 (white noise)
static float
spectralflatness(float* d, int size, float f1, float f2)
// class: AudioAnalysis
// input: d[size] = power
lower half spectrum
// pth = accumulated
power threshold:
// 0.5f
for cutoff frequency of first order lowpass spectrum
// 0.95f
for perceived upper band limit
// output: return value = cutoff
frequency as fraction of the sample rate
static float
ufc(float* d, int size, float pth = 0.95f)
// class: AudioAnalysis
// input: d[size] = power
lower half spectrum
// pth = accumulated
power threshold:
// 0.5f
for cutoff frequency of first order highpass spectrum
// 0.95f
for perceived lower band limit
// output: return value = cutoff
frequency as fraction of the sample rate
static float
lfc(float* d, int size, float pth = 0.95f)
// compare
two spectra independent of their power
// class: AudioAnalysis
// input: d[size] = magnitude
half spectrum:
// well-balanced
inclusion of weaker components
// power
half spectrum:
// changes of strong
components dominate
// c[size] = continuation
data, init with zero
// output: return value = spectral
similarity: 0 (identical) .. 1 (maximally different)
// c[size] = continuation
data
static float
spectralflux(float* d, float* c, int size)
// reacts
to a quickly rising amplitude, all times in sampling intervals
// class: AudioAnalysis
// input: d[size] = audio signal
// mintime = reduced
sensitivity for shorter rise times
// maxtime = reduced
sensitivity for longer rise times
// rtime = dead
time of detection after a transient
// th = inhibit
detection of signals with amplitude < th
// c[3] = continuation data, init with 0
// output: r[size] = transientness,
positive, indicates transient if > 1
// c[3] = continuation data
static void
transamp( float* d, float* r, float* c,
int size, float mintime = 0,
float maxtime = 50.0f, float rtime =
1000.0f, float th = 0.001f )
// reacts
to both amplitude and spectral changes
// class: AudioAnalysis
// input: d[size] = magnitude
half spectrum
// c[size] = continuation
data ,init with 0
// output: return value = transientness,
positive, indicates transient if > 1
// c[size] = continuation
data
static float
transspec(float* d, float* c, int size)
// track
ordering: higher index corresponds to higher frequency
// class: AudioAnalysis
// input: freq [size] = reassigned
frequencies obtained from "spec"
// amp[size] = reassigned
amplitudes obtained from "spec"
// tracks = number
of currently existing tracks, init with 0
// c[size + 1] = continuation
data, init with 0
// output: trkidx[tracks] = indices pointing
to freq/amp bins
// usage: freq[trkidx[2]] = frequency of track 2
// trklink[tracks] = track
links, -1 indicates new track
// usage: trklink[2] = old index of track that is now track 2
// tracks = updated
number of tracks
// c[size + 1] = continuation
data
// remarks: max(tracks) = size/2, be sure to provide enough space for
trkidx and trklink
static void
mqtracks( int* trkidx, int* trklink,
int& tracks,
float* freq, float* amp, float* c,
int size )
// .. at
both ends using cubic phase interpolation
//
amplitude changes linearly
// end point not included for
seamless concatenation of subsequent data segments
//
frequencies are specified as a fraction of the sample rate
// phases are specified in
radians and needn't be
unwrapped
// class: AudioAnalysis
// input: d[size] = existing
audio data
// fs = frequency
at d[0]
// fe = frequency
at hypothetical d[size]
// as = amplitude
at d[0]
// ae = amplitude
at hypothetical d[size]
// ps = phase
at d[0]
// pe = phase
at hypothetical d[size]
// output: d[size] = new audio
data
static void
mqsyn(float* d, int size, float fs, float fe, float as, float ae, float ps,
float pe)
// .. at
both ends and phase at the starting point
// end
point phase is the integral of the linearly interpolated frequency
//
amplitude changes linearly
// end point not included for
seamless concatenation of subsequent data segments
//
frequencies are specified as a fraction of the sample rate
// phases are specified in
radians and needn't be
unwrapped
// class: AudioAnalysis
// input: d[size] = existing
audio data
// fs = frequency
at d[0]
// fe = frequency
at hypothetical d[size]
// as = amplitude
at d[0]
// ae = amplitude
at hypothetical d[size]
// ps = phase
at d[0]
// output: d[size] = new audio
data
static float
mqsyn(float* d, int size, float fs, float fe, float as, float ae, float ps)
The sample rate must be greater than 0.
Nominal audio and control data range: -1.0f … 1.0f
All functions clip out-of-range control data to
their individual input range.
“Update”
interpolates all supplied control data in a way that it glides from the current
to the target settings over the frame. Other functions create smooth
transitions too if it makes sense from a musical perspective.
Multithreading capability for x86/AMD64/ECMA
memory models: “Update” runs in parallel with any other function called from a
separate thread.
Key
features: Crossfading between different tables, pitch-dependent table switching
for minimal aliasing.
// class: WaveOsc
// input: tablesize = size of a
single wave,
// must be a power of two and
larger than 2
// tables = waves
per wavetable, range: 1..65536
// maxpitch = maximum
oscillator frequency in Hz
// minpitch = minimum
oscillator frequency in Hz
// smprate = sample
rate in Hz
WaveOsc(int
tablesize, int tables, float maxpitch, float minpitch, float smprate)
// class: WaveOsc
// input: d[tablesize] = complex half
spectrum of wavetable data, format:
// {A1,
B1, ..., Atablesize/2, Btablesize/2} with Ax
and Bx =
// amplitude
of cosine and sine part of x-th harmonic
// idx = wave
index, range: 0..tables-1
// remarks: spectrum is normalized internally to RMS = 0.5
void
LoadTable(float* d, int idx)
// class: WaveOsc
// input: phase = set wave
pointer position, range: 0 (start) .. 1.0f (end)
void
SetPhase(float phase)
// control
characteristics: pitch / pitchmod exponential, wave linear
// class: WaveOsc
// input: samples = number of
samples to process
// invsmp = 1.0f
/ (float)samples
// pitch = static pitch control, range: 0 (minpitch) .. 1.0f
(maxpitch)
// pitchmod = pitch
modulation control, range: -1.0f .. 1.0f
// wave = wave selection control,
// range:
0 (wave[0]) .. 1.0f (wave[tables-1])
// pmod[samples] = audio
rate phase modulation input,
// range:
-2^31 (-pi) .. 2^31-1 (pi)
// output: d[samples] = audio output
void
Update( float* d, int samples, float
invsmp, float pitch, float pitchmod,
float wave, int* pmod )
Includes
prefilters (zero @ fs/2), correction postfilter, and DC trap with fc = 5 Hz.
// class: RingMod
// input: smprate = sample
rate in Hz
RingMod(float
smprate)
// class: RingMod
// input: in1[samples] = audio input 1
// in2[samples] = audio
input 2
// samples = number
of samples to process
// output: out[samples] = audio output
// remarks: anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
void
Update(float* in1, float* in2, float* out, int samples)
uniformly
distributed white noise, range: -1..1, RMS amplitude: 0.577
pink noise,
accuracy: +/- 0.3db (0.00045..0.45fs), RMS amplitude: 0.577
// class: Noise
Noise()
// class: Noise
// input: type = noise
type: 0 (white), 1 (pink)
void SetType(int
type)
// class: Noise
// input: samples = number of
samples to process
// output: out[samples] = audio output
void
Update(float* out, int samples)
Implements
clickless delay time change by this sequence:
// class: Delay
// input: maxlen = maximum
delay time in samples
// smprate = sample
rate in Hz
// zlevel = audio
thru gain during delay change
Delay(int
maxlen, float smprate, float zlevel = 1.0f)
// class: Delay
// input: dlen = delay time in samples
void
SetType(int dlen)
// class: Delay
// input: in[samples] = audio input
// samples = number
of samples to process
// output: out[samples] = audio output
void
Update(float* in, float* out, int samples)
Uses 1st
order allpass interpolation.
// class: VarDelay
// input: maxlen = maximum
delay time in samples
VarDelay(float
maxlen)
// class: VarDelay
// input: in[samples] = audio input
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// length = delay
time in samples
// output: out[samples] = audio output
// remarks: anti-denormal scheme adds output noise spectrally concentrated
at 0 and fs/2,
// amplitude ≈ 1e-15
void
Update(float* in, float* out, int samples, float invsmp, float length)
Playback:
Loops:
Fractional and integer length
allowed, switching to unlooped playback ok.
Integer length allowed, specify loop
length in PreComp, reloading through
PreComp required for switching to unlooped
playback.
No loop support for reverse playback
Reverse
Playback: transpose < 0,
abs(transpose) = transpose factor
// class: SampleOsc
// remarks: constructing the first instance takes
some extra time as a shared table is built
SampleOsc()
// class: SampleOsc
// input: d[length] = original
sample
// length = sample
length in samples
// looplen = loop
length, specify if loop end > length - 8
// srate = specify
only if known: original sample rate in Hz
// output: d[length+16] = precompensated sample
void
PreComp(float* d, int length, int looplen=0, float srate=0)
// class: SampleOsc
// input: d[length+16] = precompensated sample
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// transpose = transpose
factor
// endpos = end
position (is also loop end), range: 0..length
// startpos = ≥
0: jump to specified start position
// <
0: continue
// looplen = loop
length in samples, 0: no loop
// output: out[samples] = audio output
// return value = playback
state: 0 normal
//
1 sample end reached
// -1
sample start reached in reverse direction
// remarks: if playback state = -1 or 1, the output
is zeroed with negligible CPU load in
// subsequent calls of Update
until smpstart >= 0 is specified or the playback
// direction has changed
int Update( float* out, float* d, int samples, float
invsmp,
float transpose, int endpos, int
startpos = -1, float looplen = 0 )
// class: Envelope
// input: tmax = maximum segment time constant in
samples
// tmin = minimum segment time constant in samples
// maxsegs = maximum
number of segments
Envelope(float
tmax, float tmin = 1.0f, int maxsegs = 16)
// class: Envelope
// input: pid = global parameter ID:
// ENVELOPE_LMOD:
level modulation intensity
// ENVELOPE_RELSEG: release segment index
// segment specific
parameter ID:
// ENVELOPE_TIME
: time constant
// ENVELOPE_TMOD:
time modulation intensity
// ENVELOPE_LEVEL:
end level
// ENVELOPE_NEXT: index of next segment
// ENVELOPE_STYPE:
segment type
// value = parameter value, range:
// ENVELOPE_LMOD:
0..1.0f, linear
// ENVELOPE_RELSEG: 0..maxsegs
// ENVELOPE_TIME: 0..1.0f, exponential
// ENVELOPE_TMOD:
-1.0f..1.0f, exponential
// ENVELOPE_LEVEL:
-1.0f..1.0f, linear
// ENVELOPE_NEXT: 0..maxsegs
// ENVELOPE_STYPE:
any segment type ID
// segment type IDs:
// ENVELOPE_ST_CONCAVE = finite
line segment with rounded start
// ENVELOPE_ST_CONVEX = finite
line segment with rounded end
// ENVELOPE_ST_CLASSICATK = attack
shape of classic ADSR generator
// ENVELOPE_ST_LINEAR = finite line segment
// ENVELOPE_ST_SIGMOID = finite
s-shaped segment
// ENVELOPE_ST_INFINITEEXP = infinite
exponential asymptote
//
// segment = envelope
segment, ignored for global parameters
void
SetParam(int pid, float value, int segment = 0)
// class: Envelope
// input: pid = any preset ID
// preset IDs:
// ENVELOPE_P_ADSR = classic
ADSR, segment indices:
// 0
attack, 1 decay + sustain, 2 release
// ENVELOPE_P_ABDSR = ADSR
with linear breakpoint segment between
// attack
and decay, segment indices:
// 0
attack, 1 breakpoint, 2 decay + sustain, 3 release
void
Preset(int pno)
// class: Envelope
// input: id = event: 0
(key has been pressed)
// 1
(key has been released)
// 2
(set envelope to idle state)
// remarks: Event 0 causes the envelope to switch to
segment 0
// Event 1 causes the
envelope to switch to the release segment (see SetParam)
// Event 2 causes the
envelope to switch to the idle segment (= maxsegs)
void
Event(int id)
// class: Envelope
// input: values = number
of output data points to generate
// tmod = time modulation input, range: -1.0f..1.0f, exponential
// lmod = level modulation input, range: 0..1.0f, linear
// output: out[values] = envelope data
points
// remarks: anti-denormal scheme truncates output absolute values below
5e-11 to zero
void
Update(float* out, int values, float tmod = 0, float lmod = 0)
// class: Amp
Amp()
// class: Amp
// input: curve = control
characteristic: 0 (linear), 1 (square), 2 (quartic)
void
SetType(int curve)
// class: Amp
// input: data[samples] = audio input
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// amp = amplitude control, range: 0..1.0f
// output: data[samples] = audio output
// remarks: anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
void
Update(float* data, int samples, float invsmp, float amp)
Filter
input by first order lowpass with additional zero @ fs/2 and convert it to a
wrapped integer.
// class: AudioToPM
// input: smprate = sample
rate in Hz
// fc = filter
cutoff frequency in Hz
AudioToPM(float
smprate, float fc = 400.0f)
// class: AudioToPM
// input: in[samples] = audio input,
nominal range: -1.0f..1.0f
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// modint = modulation
intensity, range: 0..1.0f
// output: out[samples] = phase modulation
output
// remarks: An input within the nominal range is
converted to an output in the range
// -2^34..2^34 which is
wrapped to an ordinary integer. As the phase
// modulation input of an
oscillator usually maps a value of 2^31 to a
// phase shift of pi, the
maximum modulation depth is +/- 4 periods.
void
Update(float* in, int* out, int samples, float invsmp, float modint)
Generates
two orthogonal output signals with a maximum phase error of 0.017 degrees
(0.0003 radians) in the frequency range 0.00015..0.49985 * sample rate.
// class: Hilbert
Hilbert()
// class: Hilbert
// input: in[samples] = original
signal
// samples = number
of samples to process
// output: out1[samples] = phase shifted signal
// out2[samples] = phase
shifted signal leading out1 by 90 degrees
// remarks: anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
void
Update(float* in, float* out1, float* out2, int samples)
// class: Lowpass1
// input: smprate = sample
rate in Hz
// fmin = minimum -3dB cutoff frequency (Hz)
Lowpass1(float
smprate, float fmin = 5.0f)
// class: Lowpass1
// input: data[samples] = original signal
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// freq = exponential
cutoff frequency control,
// warped
to infinity in the top range,
// range: 0 (fmin) .. 1.0f (infinity)
// output: data[samples] = filtered signal
// remarks: anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
void
Update(float* data, int samples, float invsmp, float freq)
// class: Highpass1
// input: smprate = sample
rate in Hz
// fmin = minimum -3dB cutoff frequency (Hz)
Highpass1(float
smprate, float fmin = 5.0f)
// class: Highpass1
// input: data[samples] = original signal
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// freq = exponential
cutoff frequency control,
// warped to infinity in the top range,
// range: 0 (fmin) .. 1.0f (infinity)
// output: data[samples] = filtered signal
// remarks: anti-denormal scheme adds output noise spectrally
concentrated at 0 and fs/2,
// amplitude ≈ 1e-15
void
Update(float* data, int samples, float invsmp, float freq)
// class: ChambFilter
// input: smprate = sample
rate in Hz
// fmin = minimum cutoff / center frequency in Hz
// remarks: constructing the first instance takes
some extra time as a shared table is built
ChambFilter(float
smprate, float fmin = 5.0f)
// class: ChambFilter
// input: type = characteristic: 0 (lowpass), 1 (bandpass), 2 (highpass),
// 3 (peaking), 4 (notch), other (bypass)
void
SetType(int type)
// class: ChambFilter
// input: data[samples] = original signal
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// freq = exponential
cutoff / center frequency control,
// range: 0 (fmin) .. 1.0f (0.42 * sample rate)
// res = resonance
control, reciprocal characteristic,
// range:
0 (minimum) .. 1.0f (maximum)
// output: data[samples] = filtered signal
// remarks: transparent with gain = 0 dB @ res = 0, freq
= 1, type = 0
// anti-denormal scheme adds ≈
3e-15 of pink noise referenced to the input
void
Update(float* data, int samples, float invsmp, float freq, float res)
// with amplitude
compression and resonance up to self oscillation
// class: ChambFilter
// input: data[samples] = original signal
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// freq = exponential
cutoff / center frequency control,
// range: 0 (fmin) .. 1.0f (0.42 * sample rate)
// res = resonance
control, reciprocal characteristic,
// range:
0 (minimum) .. 1.0f (self-oscillation)
// output: data[samples] = filtered signal
// remarks: approximate low-level gain = -6dB @ res =
0, freq = 1, type = 0
// cascade linear versions with a total Q limited to
about 5 and this version as
// last stage to create higher order virtual analog
filters with minimum aliasing
// anti-denormal scheme adds ≈
3e-15 of pink noise referenced to the input
void
UpdateCmp(float* data, int samples, float invsmp, float freq, float res)
Key
features: Low-alias wideband nonlinearity, amplitude compression, resonance up
to self-oscillation with precise tuning, improved gain evolution with
increasing resonance.
// class: MoogFilter
// input: smprate = sample
rate in Hz
// fmin = minimum cutoff frequency in Hz
// remarks: constructing the first instance takes
some extra time as a shared table is built
MoogFilter(float
smprate, float fmin = 5.0f)
// class: MoogFilter
// input: data[samples] = original signal,
nominal range: -1.0f..1.0f
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// freq = exponential
cutoff frequency control,
// range: 0 (fmin) .. 1.0f (0.42 * sample rate)
// res = resonance
control, reciprocal characteristic,
// range:
0 (minimum) .. 1.0f (self-oscillation)
// output: data[samples] = filtered signal
// remarks: nearly transparent @ res = 0, freq = 1
// approximate low-level gain
@ res = 0, freq = 1: -6dB
// anti-denormal scheme adds ≈
3e-15 of pink noise referenced to the input
void
Update(float* data, int samples, float invsmp, float freq, float res)
// class: FMOsc
// input: maxpitch = maximum
oscillator frequency in Hz
// minpitch = minimum
oscillator frequency in Hz
// smprate = sample
rate in Hz
// remarks: constructing the first instance takes
some extra time as a shared table is built
FMOsc(float
maxpitch, float minpitch, float smprate)
// class: FMOsc
// input: phase = set
wave pointer position, range: 0 (start) .. 1.0f (end)
void
SetPhase(float phase)
// class: FMOsc
// input: out[samples] = phase modulation
signal (if add = true)
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// pitch = exponential pitch control,
// range:
0 (minpitch) .. 1.0f (maxpitch)
// amp = linear
amplitude control, range: 0..1
// fbk = linear
feedback control,
// range:
0 (none) .. 1 (maximum stable without noise)
// pmod[samples] = audio
rate phase modulation input,
// range:
-2^31 (-pi) .. 2^31-1 (pi)
// add = output
mode: true: oscillator signal + input
// false: oscillator signal only
// output: out[samples] = phase modulation
output,
// range:
-2^35..2^35-1, wrapped to integer
void
UpdateMod( int* out, int samples, float
invsmp, float pitch,
float amp, float fbk, int* pmod, bool
add = false )
// class: FMOsc
// input: out[samples] = audio signal (if
add = true)
// samples = number
of samples to process
// invsmp = 1.0f
/ (float)samples
// pitch = exponential pitch control,
// range:
0 (minpitch) .. 1.0f (maxpitch)
// amp = linear
amplitude control, range: 0..1
// fbk = linear
feedback control,
// range:
0 (none) .. 1 (maximum stable without noise)
// pmod[samples] = audio
rate phase modulation input,
// range:
-2^31 (-pi) .. 2^31-1 (pi)
// add = output
mode: true: oscillator signal + input
// false: oscillator signal only
// output: out[samples] = audio output,
range: -1.0f..1.0f
void
UpdateCar( float* out, int
samples, float invsmp, float pitch,
float amp,
float fbk, int* pmod, bool add = false)
// class: VAOsc
// input: maxpitch = maximum
oscillator frequency in Hz
// minpitch = minimum
oscillator frequency in Hz
// smprate = sample
rate in Hz
// remarks: constructing the first instance takes
some extra time as a shared table is built
VAOsc(float
maxpitch, float minpitch, float smprate)
// class: VAOsc
// input: shape = wave
shape: 0 (sawtooth), 1 (pulse)
void
SetShape(int shape)
// class: VAOsc
// input: phase = set
wave pointer position, range: 0 (start) .. 1.0f (end)
void
SetPhase(float phase)
// class: VAOsc
// input: samples = number of
samples to process
// invsmp = 1.0f
/ (float)samples
// pitch = exponential pitch control,
// range:
0 (minpitch) .. 1.0f (maxpitch)
// pwidth = linear
pulse width control (pulse shape only),
// range:
0 (0%) .. 1 (100%)
// pmod[samples] = audio
rate phase modulation input,
// range:
-2^31 (-pi) .. 2^31-1 (pi)
// output: out[samples] = audio output
void
Update(float* out, int samples, float invsmp, float pitch, float pwidth, int*
pmod)
Generates a
naive sawtooth wave without any anti-aliasing measures and is intended as a building
block for more complex oscillators and LFOs. Enters a dedicated LFO mode for
maxpitch < 0.0038*smprate that allows a minimum minpitch of 1e-7*smprate.
// class: RawSawOsc
// input: maxpitch = maximum
oscillator frequency in Hz
// minpitch = minimum
oscillator frequency in Hz
// smprate = sample
rate in Hz
RawSawOsc(float
maxpitch, float minpitch, float smprate)
// class: RawSawOsc
// input: phase = set
wave pointer position, range: 0 (start) .. 1.0f (end)
void
SetPhase(float phase)
// class: RawSawOsc
// input: samples = number of
samples to process
// invsmp = 1.0f
/ (float)samples
// pitch = exponential pitch control,
// range:
0 (minpitch) .. 1.0f (maxpitch)
// pmod[samples] = audio
rate phase modulation input,
// range:
-2^31 (-pi) .. 2^31-1 (pi)
// output: out[samples] = audio output
void
Update(float* out, int samples, float invsmp, float pitch, int* pmod)
Neural
networks and other classifiers. Principal and independent component analysis.
Provides a
layer of neurons with the following structure:
k neurons {N1,N2,...,Nk} accessible by neuron
indices 1..k
a constant bias node at neuron index 0
Data
exchange between neuron layers: connect "Neurons" by
"Axons"
Data
exchange between neurons and the outside world: directly read/write from/to
n[1..k]
Neuron
types:
NT_SIG 0 = sigmoidal nonlinearity: tanh(input sum)
NT_LIN 1 = linear: input sum
NT_BIN 2 = binary: 1.0f (input sum ≥ 0), -1.0f
(input sum < 0)
// class: Neurons
// input: k = number of neurons, default: sigmoidal
nonlinearity
// type = neuron type
// remarks: bias node n[0] is set to 1
// choose k = 4n – 1 with n any integer for fastest
processing
Neurons(int
k, int type = NT_SIG)
// class: Neurons
//
output : return value = number
of neurons, failed: -1
int
GetNofNeurons()
// class: Neurons
//
output : return value = type
of neurons
int GetType()
// class: Neurons
// input: n[k] = neurons {N1, N2,..., Nk}: n[1..k]
//
output : n[k] = neurons {N1, N2,..., Nk}: n[1..k]
// remarks: n[0] is the bias node and should not be written directly
float* n
Connect two
neuron layers, perform forward and back propagation, provide learning
functionality with continuous (default) and batch weight update.
Connection: create "Neurons"
create "Axons" that point
to them
note: you
may also connect a "Neurons" object to itself ->
FeedForward updates neurons serially
in ascending order
Evaluation: write input to n[1..k] of topmost
"Neurons"
call FeedForward top-down for every
"Axons"
read output from n[1..k] of bottom
"Neurons"
BP learn: evaluate input training vector
(reading the output is optional)
call BackPropagate(learning
rate,output training vector) for bottom "Axons"
call BackPropagate(learning rate)
bottom-up for remaining"Axons" excluding
the topmost
call BackPropagate(learning
rate,NULL,true) for topmost "Axons"
note: if bottom=topmost "Axons"
call BackPropagate(learning rate,output
training vector,true)
Apply GHA: init weights by calling
Reset(WI_ZEROBIASRNDOTHER)
loop: write input vector to n[1..k] of input "Neurons"
call FeedForward
call GHA(learning rate) -> if output
neurons are linear, then:
·
the
weights of output neurons n[1..k] converge to the ordered principal components
of the input vector set
·
higher
neuron index corresponds to less significant component
·
the
current input vector is approximated by output neuron values as
sum(i=1..k){n[i] * weights of n[i]}
·
an
input vector set with zero empirical mean yields MSE-optimal components
// class: Axons
// input: in = pointer to input neurons
// out = pointer
to output neurons
Axons( Neurons* in, Neurons* out)
// class: Axons
// input: mode = initialization mode:
// WI_ZERO
(zero all),
// WI_ZEROBIASRNDOTHER
(zero bias weights, randomize other)
// WI_RND
(randomize all), default
void Reset(int
mode = WI_RND)
// class: Axons
void
FeedForward()
// class: Axons
void
EnableBatchLearning()
// class: Axons
void
DisableBatchLearning()
// class: Axons
// input: lrate = learning
rate
// out = if
specified: process output neurons as output layer with
// training
vector out[0..output_neurons-1] = {O1,O2,..,Ok}
// in = if
true: process input neurons as input layer
void
BackPropagate(float lrate = 0, float* out = NULL, bool in = false)
// based on
the current input / output neuron values, bias weights are ignored
// class: Axons
// input: lrate = learning
rate
void
Hebb(float lrate)
// class: Axons
// input: lrate = learning
rate
void
GHA(float lrate)
// class: Axons
// input: n = input neuron index, 0: bias node
// output: d[M] = weights from input neuron n to output
neurons d[i] with
// i
= output neuron index 1..M,
// M
= number of output neurons + 1
// remarks: d[0] = 0 as bias nodes have no input
port
void
ReadInputWeights(float* d, int n)
// class: Axons
// input: n = input neuron index, 0: bias node
// output: d[M] = weights from input neuron n to output
neurons d[i] with
// i
= output neuron index 1..M,
// M
= number of output neurons + 1
// remarks: d[0] ignored as bias nodes have no input
port
void
WriteInputWeights(float* d, int n)
// class: Axons
// input: n = output neuron index > 0
// output: d[M] = weights from output neuron n to input
neurons d[i] with
// i
= input neuron index 1..M, d[0] = weight of bias node,
// M
= number of input neurons + 1
void
ReadOutputWeights(float* d, int n)
// class: Axons
// input: n = output neuron index > 0
// output: d[M] = weights from output neuron n to input
neurons d[i] with
// i
= input neuron index 1..M, d[0] = weight of bias node,
// M
= number of input neurons + 1
void
WriteOutputWeights(float* d, int n)
Provides Kohonen
feature maps, vector quantization, and k-means clustering. Nodes are initialized
from node 0 up with the first processed input vectors. Decay: value drops to
value*exp(-decay) per call of "Learn".
// class: Kohonen
// input: nodes = number
of output nodes
// vsize = input vector size
// lrate = initial learning rate
// alpha = learning rate decay
// beta = neighborhood decay (for distance = 1.0f)
// dist[N] = output
node topology (if any):
// distance
between nodes, N = nodes*nodes
// remarks: if vsize is larger than 63, make it a multiple of 4 for
fastest processing
Kohonen(int
nodes, int vsize, float lrate, float
alpha = 0, float beta = 0, float* dist = NULL)
// class: Kohonen
void
Reset()
// class: Kohonen
// input: d[vsize] = vector to
classify
// output: return value = output node
index
// remark: faster with both d aligned and vsize a
multiple of 4 larger than 63
int
GetNode(float* d)
// class: Kohonen
// input: node = output node index
// output: d[vsize] = vector of
specified output node
void
GetVector(float* d, int node)
// class: Kohonen
// input: d[vsize] = vector to
learn
// remark: faster with both d aligned and vsize a
multiple of 4 larger than 63
void
Learn(float* d)
// class: Kohonen
// input: d[vecs*vsize] = set of input
vectors,
// format:
{vector0[vsize], .. ,vectorvecs-1[vsize]}
// vecs = number of input vectors
// output: kmnode[vecs] = node associated
with vector
// return value = number
of iterations
// remark: faster with both d aligned and vsize a
multiple of 4 larger than 63
int
KMeans(float* d, int vecs, int* kmnode)
Principal and
independent component analysis.
Component
properties: orthonormal, zero mean, sign undetermined
// class: PCA
// input: d[vecs*vsize] = input vectors,
// format:
{vector0[vsize], .. , vectorvecs-1[vsize]}
// vsize = input vector size
// vecs = number of input vectors
// n = number
of components to calculate
// output: pcs[n*vsize] = principal
components,
// format:
{component0[vsize], .. , componentn-1[vsize]},
// ordering:
higher index corresponds to lower significance
// mean[vsize] = arithmetic
mean of input vectors
// var[n] = variance
explained by component as a fraction of total
// variance
// d[vecs*vsize] = residual:
original data – reconstruction from components
// return value = number
of computed components
// remarks: faster with both d and pcs aligned and
vsize a multiple of 4
static int
pcomp(float* pcs, float* mean, float* var,
float* d, int vsize, int vecs, int n)
// … so
that original vector = sum(i = 0..n-1){a[i]*component(i)} + mean + residual
// class: PCA
// input: pcs[n*vsize] = principal
components,
// format:
{component0[vsize], .. , componentn-1[vsize]}
// d[vsize] = input
vector
// vsize = input vector size
// n = number
of components to sum up
// mean[vsize] = arithmetic
mean of input vectors, add to sum if specified
// output: a[n] = component weights
// d[vsize] = residual:
original vector – (sum of components + mean)
// remarks: faster with both pcs and d aligned and
vsize a multiple of 4
static void
decompose(float* a, float* pcs, float* d, int vsize, int n, float* mean = NULL
)
// … given
by sum(i = 0..n-1){a[i]*component(i)} + mean + residual
// class: PCA
// input: d[vsize] = residual
// pcs[n*vsize] = principal
components,
// format:
{component0[vsize], .. , componentn-1[vsize]}
// a[n] = component weights
// vsize = input vector size
// n = number
of components to sum up
// mean[vsize] = arithmetic
mean of input vectors, add to sum if specified
// output: d[vsize] = reconstructed
vector
// remarks: faster with both d and pcs aligned and
vsize a multiple of 4
static void
compose(float* d, float* pcs, float* a, int vsize, int n, float* mean = NULL )
// … with
FastICA algorithm
// class: PCA
// input: pcs[n*vsize] = principal
components,
// format:
{component0[vsize], .. , componentn-1[vsize]}
// vsize = input and output vector size
// m = number
of independent components, m ≤ n
// n = number
of principal components
// output: ics[m*vsize] = independent
components,
// format:
{component0[vsize], .. , componentm-1[vsize]}
// remarks: faster with pcs aligned and vsize a
multiple of 4
static void
Encapsulates
an audio data buffer of float -1.0f..1.0f and provides file IO.
Allows safe
buffer access from a second thread via GetSafePt (any channel, true).
If a file
with supported write format is loaded and the audio buffer remains unchanged, “SaveWave”
and “AppendWave” will create an exact copy of the original audio data.
If the
buffer values range exceeds 1.0f, the data is normalized before it is written
to file.
Supported
read formats:
WAVE
(PCM8/16/24/32, FLOAT32, each standard + WAVE_FORMAT_EXTENSIBLE)
AIFF
(PCM8/16/24/32)
AIFFC
(PCM8/16/24/32, FLOAT32)
RAW
(PCM8s+u/16/24/32, FLOAT32, big and little endian)
Supported
write formats:
WAVE
(standard: PCM16, WAVE_FORMAT_EXTENSIBLE: PCM16/24, FLOAT32)
Supported
append formats:
WAVE
(PCM16/24, FLOAT32, each standard + WAVE_FORMAT_EXTENSIBLE)
Supported
speaker positions of WAVE_FORMAT_EXTENSIBLE format. Combine any of these flags
to assign channels to speakers, use 0 for undefined speaker arrangements:
SPK_FRONT_LEFT,
SPK_FRONT_RIGHT, SPK_FRONT_CENTER, SPK_LOW_FREQUENCY, SPK_BACK_LEFT, SPK_BACK_RIGHT,
SPK_FRONT_LEFT_OF_CENTER, SPK_FRONT_RIGHT_OF_CENTER, SPK_BACK_CENTER, SPK_SIDE_LEFT,
SPK_SIDE_RIGHT, SPK_TOP_CENTER, SPK_TOP_FRONT_LEFT, SPK_TOP_FRONT_CENTER,
SPK_TOP_FRONT_RIGHT, SPK_TOP_BACK_LEFT, SPK_TOP_BACK_CENTER, SPK_TOP_BACK_RIGHT
Specific
error codes:
NOFILE //
file not found
CORRUPT //
file corrupt
ERRREAD //
read file error
UNSAVED //
could not write to file at all
INCOMPLETE //
could write only a part to file
FMTERR //
no WAVE/AIFF format
NOSUPPORT //
unsupported data format
MISMATCH //
properties do not match
NODATA //
audio buffer empty
NOMEMORY //
out of audio data memory
NOMEMDEL //
NOMEMORY, audio buffer deleted
// class: AudioFile
AudioFile()
// class: AudioFile
// input: nsize = size
in frames (samples / channels)
// nchannels = number
of channels
// nresolution = resolution
in bit
// nrate = sample rate in Hz
// nspkpos = speaker
positions (see above)
int Create( unsigned int nsize = 0, unsigned int
nchannels = 1,
unsigned int nresolution = 16, unsigned
int nrate = 44100,
unsigned int nspkpos = 0 )
// class: AudioFile
// input: filename = full
pathname including extension
// offset = frame
offset for reading audio data
// frames = number
of frames to read, 0: all
// nodata = true:
read properties only,
// an
existing audio buffer will be deleted
// output: return value = error code (see
above)
// remarks: format is determined automatically
// a frame comprises a single
sample of each channel
int Load(char
*filename, unsigned int offset = 0, unsigned int frames = 0, bool nodata = false)
// class: AudioFile
// input: filename = full
pathname including extension
// offset = frame
offset for reading audio data
// frames = number
of frames to read, 0: all
// nchannels = number
of channels
// nresolution = resolution:
// 1-8
(format = false): PCM8 signed
// 1-8
(format = true): PCM8 unsigned
// 9-16: PCM16
// 17-24: PCM24
// >
24 (format = false): PCM32
// >
24 (format = true): float
// nrate = sample rate
// nspkpos = speaker
positions (see above)
// format = (see
nresolution)
// bigendian = true:
data interpreted as big endian
// output: return value = error code (see
above)
// remarks: a frame comprises a single sample of
each channel
// THIS FUNCTION HAS NOT BEEN EXHAUSTIVELY TESTED YET!
int
LoadRaw( char *filename, unsigned int
offset = 0, unsigned int frames = 0,
unsigned int nchannels = 1,unsigned
int nresolution = 16,
unsigned int nrate = 44100, unsigned
int nspkpos = 0,
bool format = false, bool bigendian
= false )
// class: AudioFile
// 16 bit
audio is saved as WAVE_FORMAT_EXTENSIBLE if channels > 2
// input: filename = full
pathname including extension
// output: return value = error code (see
above)
int
SaveWave(char *filename)
// class: AudioFile
// input: filename = full
pathname including extension
// output: return value = error code (see
above)
int
AppendWave(char *filename)
// class: AudioFile
// output: return value = size in frames
unsigned
int GetSize()
// class: AudioFile
// output: return value = sample rate in
Hz
unsigned
int GetRate()
// class: AudioFile
// output: return value = number of
channels
unsigned
int GetChannels()
// class: AudioFile
// output: return value = resolution in
bits
unsigned
int GetResolution()
// class: AudioFile
// output: return value = speaker
position (see above)
unsigned
int GetSpkPos()
// class: AudioFile
// input: channel = selected
channel, range: 0..channels-1
// lock = if true and return value not NULL:
// pointer
and properties remain unchanged until the next
// call
to “GetSafePt” with false, Create / Load / LoadRaw
// are
locked in a non-blocking wait loop during this time!
// output: return value = valid pointer
or NULL
float*
GetSafePt(unsigned int channel = 0, bool lock = false)
XY chart
class with optional Z color.
Special
features: mouse-controlled selection, highlighting.
Send
message on mouse event inside diagram:
msg = msgbase, wParam = msgtag,
lParam = [x(31..18) y(17..4) release(3)
press(2) left(1) right(0)]
Get value pairs
corresponding to mouse position from lParam with MsgToX and MsgToY.
Requires Windows and a compiler with MFC
support (windows, graphics, mouse)!
// class: Chart
Chart()
// call
once after creation
// class: Chart
// input: pwnd = parent window
// wframe
= chart window frame
// labelsize = label
font size
// titlesize = title
font size
// font = font type
// title = true: chart has title
// xlabel = true:
chart has x axis label
// xtickl = true: chart has x axis tick labels
// ylabel = true:
chart has y axis label
// ytickl = true: chart has y axis tick labels
// zlabel = true:
chart has z color label
// frame = true:
diagram has frame
// dcenter = true:
center labels at diagram frame
// false: center labels at window frame
// border = minimum
chart border in pixels
// ytadj = adjust y tick space in +/- pixels
// wndstyle = window
style: 0..3 child, 4 free float, 5 popup
// msgtag = tag
to identify chart in mouse messages
// msgbase = chart
mouse event message ID
// output: return value = 0: success,
other: failed
int Init( CWnd* pwnd = NULL,
CRect wframe = CRect(0,0,500,300),
int labelsize = 10, int titlesize =
14,
CString font = "Arial",
bool title = true,
bool xlabel = true, bool xtickl =
true,
bool ylabel = true, bool ytickl =
true,
bool zlabel = false,
bool frame = true,
dcenter = false,
int border = 0,
int ytadj = 0,
int wndstyle = 0,
unsigned int msgtag = 0, unsigned
int msgbase = WM_APP )
// class: Chart
// input: title = chart
title
// xlabel = x
axis label
// ylabel = y axis label
// zlabel = z
color label, "": no rainbow
// xdivs = number of x axis
divisions
// xmin = x axis minimum value
// xmax = x axis maximum value
// xlog = true: x axis has
logarithmic scale
// xgrid = true: x axis has
grid
// ydivs = number of y axis
divisions
// ymin = y axis minimum value
// ymax = y axis maximum value
// ylog = true: y axis has
logarithmic scale
// ygrid = true: y axis has
grid
// zmin = z color minimum value
// zmax = z color maximum value
// zlog = true: z color has
logarithmic scale
// bgcolor = chart
background color
// txtcolor = label
and title text color
// frmcolor = diagram
frame color
// dgcolor = diagram
background color
// gdcolor = diagram
grid color
// selmask = selection
XOR color mask, 0: invert selection
// dxborder = diagram
passive x border in pixels
// dyborder = diagram passive
y border in pixels
void
DrawChart( CString title =
"Untitled",
CString xlabel
= "x",
CString ylabel
= "y",
CString zlabel
= "z",
int xdivs = 10,
float xmin = -1.0f,
float xmax = 1.0f,
bool xlog = false,
bool xgrid = true,
int ydivs = 10,
float ymin = -1.0f,
float ymax = 1.0f,
bool ylog = false,
bool ygrid = true,
float zmin = 0,
float zmax = 1.0f,
bool zlog = false,
COLORREF
bgcolor = RGB(80,80,80),
COLORREF txtcolor = RGB(200,200,200),
COLORREF frmcolor = RGB(140,140,140),
COLORREF dgcolor = RGB(180,180,180),
COLORREF
gdcolor = RGB(160,160,160),
COLORREF selmask = RGB(64,64,64),
int dxborder = 0,
int dyborder = 0 )
// class: Chart
// input: xdata[size] = x data array
// ydata[size] = y
data array
// zdata[size] = z
data array
// size = number of data points
// pltcolor = plot
color
// hidden = true: plot to memory, no chart update
// false: plot and update chart from memory
void
PlotLine( float* xdata, float* ydata, int
size, COLORREF pltcolor = RGB(0,0,255),
bool hidden = false )
void
PlotLine( float* xdata, float* ydata, float*
zdata, int size, bool hidden = false)
// class: Chart
// input: xdata[size] = x data array
// ydata[size] = y
data array
// zdata[size] = z
data array
// size = number of data points
// pltcolor = plot
color
// hidden = true: plot to memory, no chart update
// false: plot and update chart from memory
void
PlotPoint(float* xdata, float* ydata, int size, COLORREF pltcolor =
RGB(0,0,255),
bool hidden = false )
void
PlotPoint(float* xdata, float* ydata, float* zdata, int size, bool hidden =
false)
// class: Chart
// input: xdata[size] = x data array
// ydata[size] = y
data array
// zdata[size] = z
data array
// size = number of data points
// pltcolor = plot
color
// diameter = circle
diameter in pixels
// hidden = true: plot to memory, no chart update
// false: plot and update chart from memory
void
PlotCircle(float* xdata, float* ydata, int size, COLORREF pltcolor =
RGB(0,0,255),
unsigned int diameter = 5, bool hidden = false)
void
PlotCircle(float* xdata, float* ydata, float* zdata, int size,
unsigned int diameter = 5, bool
hidden = false)
// class: Chart
// input: xdata[size] = x data array
// ydata[size] = y
data array
// zdata[size] = z
data array
// size = number of data points
// pltcolor = plot
color
// thickness = bulk
width: 1 + 2*thickness pixels
// hidden = true: plot to memory, no chart update
// false: plot and update chart from memory
void
PlotVBulk(float* xdata, float* ydata, int size,
COLORREF pltcolor = RGB(0,0,255), int
thickness = 0, bool hidden = false)
void
PlotVBulk(float* xdata, float* ydata, float* zdata, int size, int thickness = 0,
bool hidden = false)
// class: Chart
// input: xdata[size] = x data array
// y1data[size] = y
data array, interval endpoint 1
// y2data[size] = y
data array, interval endpoint 2
// zdata[size] = z
data array
// size = number of data points
// pltcolor = plot
color
// thickness = bulk
width: 1 + 2*thickness pixels
// hidden = true: plot to memory, no chart update
// false: plot and update chart from memory
void
PlotVBulkInterval(float* xdata, float* y1data, float* y2data, int size,
COLORREF pltcolor = RGB(0,0,255), int
thickness = 0, bool hidden = false)
void
PlotVBulkInterval(float* xdata, float* y1data, float* y2data, float* zdata,
int size, int thickness = 0, bool hidden =
false)
// class: Chart
// input: x1data[size] = x data array,
interval endpoint 1
// x2data[size] = x
data array, interval endpoint 2
// ydata[size] = y
data array
// zdata[size] = z
data array
// size = number of data points
// pltcolor = plot
color
// thickness = bulk
width: 1 + 2*thickness pixels
// hidden = true: plot to memory, no chart update
// false: plot and update chart from memory
void
PlotHBulkInterval(float* x1data, float* x2data, float* ydata, int size,
COLORREF pltcolor = RGB(0,0,255), int
thickness = 0, bool hidden = false)
void
PlotHBulkInterval(float* x1data, float* x2data, float* ydata, float* zdata,
int size, int thickness = 0, bool
hidden = false)
// class: Chart
// input: hidden = true: clear memory, no chart update
// false: clear chart and memory
void
Clear(bool hidden = false)
// class: Chart
// input: color = frame
color
// hidden = true: draw to memory, no chart update
// false: draw and update chart from memory
void
DrawFrame(COLORREF color, bool hidden = false)
// reversibly
highlight a rectangular zone of the diagram, call again to unlight
// class: Chart
// input: x1, x2 = x
interval of zone
// y1, y2 = y interval of zone
// mask = color XOR mask, 0: invert zone
// hidden = true: draw to memory, no chart update
// false: draw and update chart from memory
void
Highlight(float x1, float x2, float y1, float y2, COLORREF mask, bool hidden =
false)
// select a
rectangular zone of the diagram
// class: Chart
// input: x1, x2 = x
interval of zone
// y1, y2 = y interval of zone
// hidden = true: draw to memory, no chart update
// false: draw and update chart from memory
void
Select(float x1, float x2, float y1, float y2, bool hidden = false)
// class: Chart
// input: hidden = true: draw to memory, no chart update
// false: draw and update chart from memory
void
Unselect( bool hidden = false)
// class: Chart
// output: r[4] = values of selected segment,
// format: {xmin, xmax, ymin, ymax}
// return value = true: valid selection exists
// false: currently unselected
bool GetSelection(float*
r)
// mouse
messages are sent on mouse moves and clicks inside the diagram
// class: Chart
void
EnableMouseResponse()
// class: Chart
void
DisableMouseResponse()
// mouse
actions inside the diagram will call Select() and Unselect() appropriately
// class: Chart
// input: x = true: enable
selection on x axis
// false: disable selection on x axis,
// GetSelection()
returns diagram x range
// y = true: enable selection on y axis
// false: disable selection on y axis,
// GetSelection()
returns diagram y range
void
EnableMouseSelection(bool x = true, bool y = true)
// class: Chart
void
DisableMouseSelection()
// class: Chart
// input: lParam = mouse
message
// output: return value = diagram x value
associated with mouse message
float
MsgToX(unsigned int lParam);
// class: Chart
// input: lParam = mouse
message
// output: return value = diagram y value
associated with mouse message
float
MsgToY(unsigned int lParam)
// class: Chart
// input: lParam = mouse
message
// output: return value = true: right button active
bool RightButton(unsigned
int lParam)
// class: Chart
// input: lParam = mouse
message
// output: return value = true: left button active
bool
LeftButton(unsigned int lParam)
// class: Chart
// input: lParam = mouse
message
// output: return value = true: button has just been pressed
bool
Pressed(unsigned int lParam)
// class: Chart
// input: lParam = mouse
message
// output: return value = true: button has just been released
bool Released(unsigned
int lParam)
// useful
to draw anything into the diagram, see also GetDiagramFrame()
// class: Chart
// output: return value = pointer to
device context (NULL if uninitialized)
CDC*
GetSafeDeviceContext()
// useful
to draw anything into the diagram, see also GetSafeDeviceContext()
// class: Chart
// output: return value = diagram frame
(0 if uninitialized)
CRect
GetDiagramFrame()