ICST DSP Library

Version 1.2

 

 

 

Beat Frei

Institute for Computer Music and Sound Technology (ICST)

Zurich University of the Arts

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

 

 


Abstract

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.

 

Key Features

·        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.

 

Usage

·        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)

 

Configuration

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.

 

Special Issue: Alignment

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.

 

Special Issue: Denormals

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.

 

Profile

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

 

 

License

Copyright © 2008, 2009, 2010, Zurich University of the Arts, Beat Frei. All rights reserved.

 

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

 

  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

 

THIS SOFTWARE IS PROVIDED BY THE ZURICH UNIVERSITY OF THE ARTS ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ZURICH UNIVERSITY OF THE ARTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 

Beta Release Note

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

 

General Considerations

 

“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

 

Function Groups

Transforms

Signals and Windows

Elementary Real Array Operations

Elementary Complex Array Operations

Special Array Operations

Elementary Matrix Operations

Filters and Delays

System Analysis

Statistics

Storage

Auxiliary

Circular Buffer

Functions

Polynomials

Interpolation

Differential Equations

Filter Design

Spectral Analysis

Decompose to Arbitrary Functions

Cepstral Analysis

Linear Prediction

Adaptation

Feature Extraction

Resynthesis

Audio Synthesis Oscillators

Audio Synthesis Filters + Amplifiers

Audio Synthesis Miscellaneous

Audio Synthesis Control

Neurons

Axons

Kohonen Feature Maps + Vector Quantization

Principal + Independent Component Analysis

Read + Write Audio Files

Visualization (2D)

 


Functions

Transforms

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

 

Signals and Windows

Create linear segment

Create asymptotic exponential segment

Create logarithmically spaced data points

Create sine wave

Create counter-clockwise rotating complex phasor

Create sine wave with linearly changing frequency

Create sine wave with exponentially changing frequency

Create sawtooth wave

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 Sinc window

Create triangular window

Create Hann (a.k.a. raised cosine) window

Create Hamming window

Create Blackman window

Create 3-term Blackman-Harris window with -67 dB side lobe

Create 4-term Blackman-Harris window with -92 dB side lobe

Create 5-term flat top window

Create Gaussian window

Create Kaiser window

 

Elementary Real Array Operations

            Sum of elements

Energy

Power

RMS value

L2 norm

Get index of maximum

Get index of minimum

Get index of element with largest difference to reference value

Get index of element with smallest difference to reference value

Normalize to L2 norm = 1

Absolute value

Fast sign

Fast reciprocal

Fast square root

Fast cosine

Fast sine

Fast inverse tangent

Fast natural logarithm of absolute value

Fast exponential

Reverse element order

Cumulative sum

Differentiate

Integrate incrementally

Set

Prune

Limit to interval

Copy

Exchange

Element-wise maximum

Element-wise minimum

Add

Subtract

Multiply

Multiply-Accumulate

Dot product

Squared distance

Convolution

Fast convolution (FFT-based)

Cross correlation

Fast cross correlation (FFT-based)

Unbiased autocorrelation

Biased autocorrelation

Fast biased autocorrelation (FFT-based)

 

Elementary Complex Array Operations

            Sum of elements

Energy

Power

RMS value

L2 norm

Normalize to L2 norm = 1

Conjugate

Reciprocal

Magnitude

Squared magnitude

Argument in radians

Prune

Set

Copy

Exchange

Add

Subtract

Multiply

Multiply-Accumulate

Dot product

Get real part

Get imaginary part

Real to complex

Polar to cartesian

Change magnitude

 

Special Array Operations

            Copy float to double

            Copy double to float

Copy and clip float to int

Copy int to float

Interleave data

Deinterleave data

Evaluate polynomial

Evaluate polynomial with complex argument

Evaluate Chebyshev series

Find peaks

Find dips

Inverse of tabulated monotonically increasing function

Sort for descending order while keeping track of correspondences

Create subset while keeping track of correspondences

Unwrap phase

Map interval linearly to another

Linear crossfade

Custom crossfade

Power crossfade

Linear interpolated cyclic table lookup

 

Elementary Real Matrix Operations

            Determinant

            Trace

            Transpose

            Create identity matrix

            Multiply by constant

Add two matrices

Subtract one matrix from another

Multiply matrix by column vector

Multiply transpose of matrix by column vector

Multiply two matrices

Multiply inverse of matrix by another matrix

 

Filters and Delays

            Static delay

Static FIR filter

Static IIR filter

Static first order section

Time-varying first order section

Static second order section

Median filter (3-point)

Median filter (5-point)

 

System Analysis

            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

 

Statistics

Arithmetic mean

Arithmetic mean of a distribution

Unbiased variance

Unbiased variance of a distribution

Average absolute deviation

Median

Median absolute deviation

Quantile

Geometric mean

Harmonic mean

Covariance

Pearson correlation

Spearman’s rank correlation

Kendall tau rank correlation

Linear regression

Fit polynomial

Histogram

Make quantile-quantile pairs

Find outliers

Grubbs test

T test

F test

Levene / Brown-Forsythe test

Chi square test

Binomial test

CDF of the normal distribution

CDF complement of the normal distribution

Confidence interval of normally distributed data

 

Storage

Save data to file as raw series of floats

Load data from file as raw series of floats

Get number of elements in file

 

Auxiliary

Return closest power of two greater or equal to i

            Allocate array aligned for SSE

            Free array aligned for SSE

            Declare variable or array aligned for SSE (Macro)

            Prepare transforms for maximum performance

            Unprepare transforms

 

Circular Buffer

Create buffer

Reset pointers and fill buffer with zeros

Write

Read

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

 

Functions

Quick and dirty atan2

Quick and dirty hyperbolic tangent

Quick and dirty medium range exponential

Fast medium range exponential

Inverse hyperbolic sine

Inverse hyperbolic cosine

Inverse hyperbolic tangent

Error function

Error function complement

Probit function

Natural logarithm of the gamma function

Regularized gamma function

Regularized beta 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 root of a function

 

Polynomials

            Find complex roots of polynomial with real coefficients

            Construct polynomial from complex roots

Symbolic addition

Symbolic subtraction

Symbolic multiplication by constant

Symbolic multiplication of two polynomials

Symbolic scaling

Symbolic substitution

Reverse

Symbolic differentiation

Symbolic integration

Create Chebyshev polynomial

Convert Chebyshev series to polynomial

Convert polynomial to Chebyshev series

 

Interpolation

Linear interpolation

Fit parabola and return extremum

Fit cubic polynomial

Cubic interpolation

Find coefficients of Chebyshev approximation

 

Differential Equations

            Solve using forth order Runge-Kutta method

 

Filter Design

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

 

Spectral Analysis

High resolution spectral analysis

Prepare parameters for high resolution spectral analysis

 

Analysis by Decomposition to Arbitrary Functions

            Matching pursuit

 

Cepstral Analysis

            Convert magnitude spectrum to real cepstrum

Convert real cepstrum to magnitude spectrum

Get Mel frequency cepstral coefficients (MFCC)

 

Linear Prediction

            Levinson-Durbin recursion

            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

 

Adaptation

Normalized LMS algorithm

Tracking demodulator

 

Feature Extraction

Envelope follower

Fundamental frequency detector

Verify fundamental frequency against high resolution spectrum

Extract harmonic spectrum and inharmonicity

Count zero crossings

Spectral flatness within a defined frequency band

Upper cutoff frequency

Lower cutoff frequency

Spectral flux

Amplitude-based attack transient detector

Spectrum-based general transient detector

Partial tracker after McAulay/Quatieri

 

Resynthesis

Add time-varying cosine wave to existing data matching frequencies and phases

Add time-varying cosine wave to existing data matching frequencies

 

Wavetable Oscillator

Create oscillator

Load wavetable with spectrum

Set phase synchronously

Audio and control update

 

FM Oscillator

Create oscillator

Set phase synchronously

Modulator audio and control update

Carrier audio and control update

 

Virtual Analog Oscillator

Create oscillator

Set wave shape

Set phase synchronously

Audio and control update

 

Sample Playback Oscillator

Create oscillator

Prepare sample

Audio and control update

 

Raw Sawtooth Oscillator

            Create oscillator

Set phase synchronously

Audio and control update

 

Noise Generator

Create noise generator

Set noise type

Audio and control update

 

First Order Lowpass Filter

Create filter

Audio and control update

 

First Order Highpass Filter

Create filter

Audio and control update

 

Second Order Multimode Filter

Create filter

Set filter characteristic

Audio and control update (linear version)

Audio and control update (virtual analog version)

 

Virtual Analog Moog Low Pass Filter

Create filter

Audio and control update

 

Controlled Amplifier

Create amplifier

Set control characteristic

            Audio and control update

 

Static Delay Line

Create delay

Set delay time

            Audio and control update

 

Time-varying Delay Line

Create delay

Audio and control update

 

Convert Audio to Phase Modulation Signal

Create converter

Audio and control update

 

Ring Modulator

Create ring modulator

Audio and control update

 

Hilbert Transformer

            Create transformer

Audio and control update

 

Envelope Generator

Create envelope

Set parameter

Set preset envelope

Create envelope event

Control update

 

Neurons

Create neurons

Get number of neurons

Get type of neurons

Content of neurons

 

Axons

Create axons with randomized weights

Inititalize weights

Update output neurons

Enable batch learning

Disable batch learning and update weights

Back propagation step

Simple Hebbian weight update

Update weights by generalized Hebbian algorithm

Read input neuron weights

Write input neuron weights

Read output neuron weights

Write output neuron weights

 

Kohonen

Create Kohonen, vector quantization or k-means map

Reset network to state after construction

Get output node corresponding to input vector

Get vector of output node

Learn input vector

Perform k-means clustering

 

Principal Component Analysis

Get principal components

Decompose vector into sum of components

Compose vector as sum of components

Convert principal to independent components

 

AudioFile

Create audio file object

Create new audio file

Load file

Load raw file and specify format

Save as WAVE file

Append to existing WAVE file

Get size

Get sample rate

Get number of channels

Get resolution

Get speaker positions

Get safe pointer to audio data

 

Chart

Create chart

Initialize chart

Draw chart

Plot line diagram

Plot point diagram

Plot filled circle diagram

Plot vertical bulk diagram

Plot vertical bulk interval diagram

Plot horizontal bulk interval diagram

Clear diagram

Draw diagram frame

Highlight diagram zone

Select diagram segment

Unelect diagram segment

Get selection

Enable mouse response

Disable mouse response

Enable selection by mouse

Disable selection by mouse

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

Get diagram frame


 

BlkDsp

 

Standard block signal processing.

Transforms

// Fast forward (fft) and inverse (ifft) Fourier transform of complex data

// 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)

 

// Fast Fourier transform of real data to lower half spectrum

// 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)

 

// Fast inverse Fourier transform of lower half spectrum to real data

// 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)

 

// Fast Fourier transform of symmetrical lower half real data to lower half spectrum

// 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)

 

// Fast inverse Fourier transform of lower half spectrum to symmetrical lower half real

// data

// 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)

 

// Fast forward (dct) and inverse (idct) type 2 cosine transform

// 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)

 

// Fast forward (dst) and inverse (idst) type 2 sine transform

// 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)

 

// Get k-th bin of the Fourier transform of arbitrary size real data

// class:             BlkDsp

// input:             d[size] = {re0, re1,.., resize-1}    

// output:          return value = {rek, imk}

static cpx goertzel(float* d, int size, int k)

 

// Forward Haar wavelet transform

// 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)

 

// Inverse Haar wavelet transform

// 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)

 


Signals and Windows

 

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.

 

 

// Create linear segment

// 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)

 

// Create asymptotic exponential segment

// 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)

 

// Create logarithmically spaced data points

// 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)

 

// Create sine wave

// 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)

 

// Create counter-clockwise rotating complex phasor

// 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)

                                                                                                                                 

// Create sine wave with linearly changing frequency

// 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)

 

// Create sine wave with exponentially changing frequency

// 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)

 

// Create sawtooth wave

// 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       )

 

// Create uniformly distributed noise in the interval [-1,1]            

// 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)

 

// Create exponentially distributed noise with variance 1             

// class:             BlkDsp

// output:          d[size]

static void enoise(float* d, int size)

 

// Create standard Cauchy distributed noise                     

// class:             BlkDsp

// output:          d[size]

static void cnoise(float* d, int size)

 

// Create Gaussian noise with variance 1

// 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)

 

// Create Sinc window

// class:             BlkDsp

// input:             periods            =          number of periods                   

// output:          d[size]

static void sinc(float* d, int size, double periods)

 

// Create triangular window 

// class:             BlkDsp

// output:          d[size]

static void triangle(float* d, int size)

 

// Create Hann (a.k.a. raised cosine) window        

// class:             BlkDsp

// output:          d[size]

static void hann(float* d, int size)

 

// Create Hamming window  

// class:             BlkDsp

// output:          d[size]

static void hamming(float* d, int size)

 

// Create Blackman window

// class:             BlkDsp

// output:          d[size]

static void blackman(float* d, int size)

 

// Create 3-term Blackman-Harris window with -67 dB side lobe

// class:             BlkDsp

// output:          d[size]

static void bhw3(float* d, int size)

 

// Create 4-term Blackman-Harris window with -92 dB side lobe            

// class:             BlkDsp

// output:          d[size]

static void bhw4(float* d, int size)

 

// Create 5-term flat top window

// class:             BlkDsp

// output:          d[size]

static void flattop(float* d, int size)

 

// Create Gaussian window

// class:             BlkDsp

// input:             sigma               =          width parameter

// output:          d[size]

static void gauss(float* d, int size, double sigma)

 

// Create Kaiser window

// class:             BlkDsp

// input:             alpha                =          width parameter

// output:          d[size]

static void kaiser(float* d, int size, double alpha)

 


Elementary Real Array Operations

// Sum of elements

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          sum of elements

// remarks:        faster with d aligned

static float sum(float* d, int size)

 

// Energy

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          <d,d>

// remarks:        faster with d aligned

static float energy(float* d, int size)

 

// Power

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          <d,d>/size

// remarks:        faster with d aligned

static float power(float* d, int size)

 

// RMS value

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          sqrt(<d,d>/size)

// remarks:        faster with d aligned

static float rms(float* d, int size)

 

// L2 norm

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          sqrt(<d,d>)

// remarks:        faster with d aligned

static float norm(float* d, int size)

 

// Get index of maximum

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          index of most positive element of d

static int maxi(float* d, int size)

 

// Get index of minimum

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          index of most negative element of d

static int mini(float* d, int size)

 

// Get index of element with largest difference to reference value

// 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)

 

// Get index of element with smallest difference to reference value

// 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)

 

// Normalize to L2 norm = 1

// 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)

 

// Absolute value

// class:             BlkDsp

// input:             d[size]

// output:          d[size]

// remarks:        faster with d aligned

static void abs(float* d, int size)

 

// Fast sign

// 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)

 

// Fast reciprocal

// class:             BlkDsp

// input:             d[size]

// output:          d[size]

static void finv(float* d, int size)

 

// Fast square root

// class:             BlkDsp

// input:             d[size]

// output:          d[size]

static void fsqrt(float* d, int size)

 

// Fast cosine

// 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)

 

// Fast sine

// 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)

 

// Fast natural logarithm of absolute value

// 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)

 

// Fast exponential

// 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)

 

// Reverse element order

// class:             BlkDsp

// input:             d[size]

// output:          d[size]

static void reverse(float* d, int size)

 

// Cumulative sum

// 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)

 

// Differentiate

// 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)

 

// Integrate incrementally

// 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)

 

// Set

// 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)

 

// Prune

// 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)

 

// Limit to interval

// 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)

 

// Copy

// class:             BlkDsp

// input:             r[size]

// output:          d[size]

// remarks:        arrays r and d may overlap

static void copy(float* d, float* r, int size)       

 

// Exchange

// 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)

 

// Element-wise maximum

// 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)

 

// Element-wise minimum

// 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)

 

// Add

// .. 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)

 

// Subtract

// .. 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)

 

// Multiply

// .. 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)

 

// Multiply-Accumulate

// 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)

 

// Dot product

// 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)

 

// Squared distance

// 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)

 

// Convolution

// 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)

 

// Fast convolution (FFT-based)

// 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)

 

// Cross correlation

// 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)                

 

// Fast cross correlation (FFT-based)

// 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)

 

// Unbiased autocorrelation

// 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)        

 

// Biased autocorrelation

// 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)

 

// Fast biased autocorrelation (FFT-based)

// 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)

 


Elementary Complex Array Operations

 

See General Considerations for the representation of complex numbers and vectors.

 

 

// Sum of elements

// class:             BlkDsp

// input:             d[2*size]

// output:          return value      =          sum of elements

// remarks:        faster with d aligned

static cpx cpxsum(float* d, int size)

 

// Energy

// class:             BlkDsp

// input:             d[2*size]

// output:          return value      =          <d,d*>

// remarks:        faster with d aligned

static float cpxenergy(float* d, int size)

 

// Power

// class:             BlkDsp

// input:             d[2*size]

// output:          return value      =          <d,d*>/size

// remarks:        faster with d aligned

static float cpxpower(float* d, int size)

 

// RMS value

// 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)

 

// L2 norm

// class:             BlkDsp

// input:             d[2*size]

// output:          return value      =          sqrt(<d,d*>)

// remarks:        faster with d aligned

static float cpxnorm(float* d, int size)

 

// Normalize to L2 norm = 1

// 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)

 

// Conjugate

// class:             BlkDsp

// input:             d[2*size]

// output:          d[2*size]

// remarks:        faster with d aligned

static void cpxconj(float* d, int size)

 

// Reciprocal

// 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)

 

// Magnitude

// 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)

 

// Squared magnitude

// 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)

 

// Argument in radians ( = two-argument inverse tangent)

// 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)

 

// Prune

// 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)

 

// Set

// 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)

 

// Copy

// 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)  

 

// Exchange

// 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)

 

// Add

// .. 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)

 

// Subtract

// .. 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)

 

// Multiply

// .. 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)

 

// Multiply-Accumulate

// 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)

           

// Dot product

// 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)

 

// Get real part

// class:             BlkDsp

// input:             d[2*size]

// output:          re[size]

// remarks:        re = d is ok

static void cpxre(float* re, float* d, int size)

 

// Get imaginary part

// class:             BlkDsp

// input:             d[2*size]

// output:          im[size]

// remarks:        im = d is ok

static void cpxim(float* im, float* d, int size)

 

// Real to complex

// 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)

 

// Polar to cartesian

// 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)

 

// Change magnitude

// 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)


Special Array Operations

// Copy float to double

// class:             BlkDsp

// input:             r[size]

// output:          d[size]

static void ftod(double* d, float* r, int size)

 

// Copy double to float

// class:             BlkDsp

// input:             r[size]

// output:          d[size]

static void dtof(float* d, double* r, int size)

 

// Copy and clip float to int

// 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)

 

// Copy int to float

// class:             BlkDsp

// input:             r[size]

// output:          d[size]

//                      faster with both d and r aligned

static void itof(float* d, int* r, int size)

 

// Interleave data

// 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)

 

// Deinterleave data

// 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)

 

// Evaluate polynomial

// 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)

 

// Evaluate polynomial with complex argument

// 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)

 

// Evaluate Chebyshev series

// 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)

 

// Find peaks

// 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)   

 

// Find dips

// 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)

 

// Inverse of tabulated monotonically increasing function

// 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)

 

// Sort for descending order while keeping track of correspondences

// 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)

 

// Create subset while keeping track of correspondences

// 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)

                                                                                                                                            

// Unwrap phase

// 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)

 

// Map interval linearly to another

// 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)

 

// Linear crossfade

// 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)

 

// Custom crossfade

// (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)

 

// Power crossfade

// 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)

 

// Linear interpolated cyclic table lookup

// 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)

 


Elementary Real Matrix Operations

 

See General Considerations for the representation of vectors and matrices.

 

 

// Determinant of square matrix a(n,n)

// class:             BlkDsp

// input:             a[n*n]             

// output:          return value      =          determinant of a            

static float mdet(float* a, int n)

 

// Trace of square matrix a(n,n)

// class:             BlkDsp

// input:             a[n*n]             

// output:          return value      =          trace of a          

static float mtrace(float* a, int n)

 

// Transpose matrix a(m,n)

// class:             BlkDsp

// input:             a[m*n]            

// output:          a[n*m]

static void mxpose(float* a, int m, int n)

 

// Create identity matrix a(n,n)

// class:             BlkDsp

// output:          a[n*n]

static void mident(float* a, int n)

 

// Multiply matrix a(m,n) by constant c

// 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)

 

// Add two matrices a(m,n) and b(m,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)

 

// Subtract matrix b(m,n) from matrix a(m,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)


// Multiply matrix a(m,n) by column vector r(n) to get column vector d(m)

// 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)

 

// Multiply transpose of matrix a(m,n) by column vector r(m) to get column vector d(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)

 

// Multiply matrix a(m,n) by matrix b(n,p) to get matrix c(m,p)

// 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)

 

// Multiply inverse of matrix b(m,m) by matrix a(m,n) to get matrix a(m,n)

// 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)


Filters and Delays

// Static delay

// 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)

 

// Static FIR filter

// 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)

 

// Static IIR filter

// 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)

 

// Static first order section

// 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)

 

// Time-varying first order section

// 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)

 

// Static second order section

// 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)

                       

// Median filter (3-point)

// 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)

 

// Median filter (5-point)

// 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)

 


System Analysis

// Frequency response of discrete-time LTI system

// 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     )          

 

// Frequency response of continuous-time LTI system

// 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 )

 

// Group delay of discrete-time LTI system

// 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)

 

// Group delay of continuous-time LTI system

// 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)


Statistics

// Arithmetic mean

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          arithmetic mean of d

// remarks:        faster with d aligned

static float mean(float* d, int size)

 

// Arithmetic mean of a distribution

// 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)

 

// Unbiased variance

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          unbiased variance of d

// remarks:        faster with d aligned

static float var(float* d, int size)

 

// Unbiased variance of a distribution

// 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)

 

// Average absolute deviation

// 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)

 

// Median

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          median of d

static float median(float* d, int size)

 

// Median absolute deviation

// 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)

 

// Quantile

// 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)

 

// Geometric mean

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          geometric mean of d

static float gmean(float* d, int size)

 

// Harmonic mean

// class:             BlkDsp

// input:             d[size]

// output:          return value      =          harmonic mean of d

static float hmean(float* d, int size)

 

// Covariance

// 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)

 

// Pearson correlation

// 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)

 

// Spearman’s rank correlation

// 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)

 

// Kendall tau rank correlation

// class:             BlkDsp

// input:             x[size], y[size]

// output:          return value      =          Kendall tau rank correlation of x and y

// remarks:        consider cspear as a faster alternative for large size

static float ckend(float* x, float* y, int size)

 

// Linear regression

// 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           )

 

// Fit polynomial

// 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)

 

// Histogram

// 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)

 

// Make quantile-quantile pairs

// 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)

 

// Find outliers

// 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)

 

// Grubbs test

// 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)

 

// T test

// 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                     )

 

// F test

// 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)

 

// Levene / Brown-Forsythe test

// 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)

 

// Chi square test

// 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)

 

// Binomial test

// 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)

 

// Cumulative distribution function of the normal distribution

// 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)

 

// Cumulative distribution function complement of the normal distribution

// 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)

 

// Confidence interval of normally distributed data

// 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)

 


Storage

// Save data to file as raw series of floats

// 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)

 

// Load data from file as raw series of floats

// 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)

                                              

// Get number of elements in file

// class:             BlkDsp

// input:             filename           =          full pathname including extension         

// output:          return value      =          number of elements in file, -1 on failure

static int getsize(char* filename)

 


Auxiliary

// Return closest power of two greater or equal to i

// class:             BlkDsp

static int nexthipow2(int i)

 

// Allocate array aligned for SSE

// 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)

 

// Free array aligned for SSE

// 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)

 

// Declare variable or array aligned for SSE (Macro)

// 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

 

// Prepare transforms for maximum performance

// 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()

 

// Unprepare transforms

// 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()


 

CircBuffer

 

Fixed length circular buffer.

 

A method call that would cause an overflow condition has no effect.

 

// Create buffer

// class:             CircBuffer

// input:             size                  =          buffer size in elements  

CircBuffer(int size = 1)

 

// Reset pointers and fill buffer with zeros

// class:             CircBuffer

void Reset()

 

// Write

// 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)

 

// Read

// 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)

 

// Read at an offset from the current read pointer without advancing it

// 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)

 

// Get maximum number of elements that can be read without causing an overflow

// class:             CircBuffer

int GetReadSize()

 

// Get maximum number of elements that can be written without causing an overflow

// class:             CircBuffer

int GetWriteSize()


 

SpecMath

 

Signal processing specific mathematical functions and routines.

Functions

// Quick and dirty atan2

// 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)

 

// Quick and dirty hyperbolic tangent

// 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)

 

// Quick and dirty medium range exponential

// 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)

 

// Fast medium range exponential

// 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)

static void fexp(float* x)


// Inverse hyperbolic sine

// class:             SpecMath

// input:             x

// output:          return value      =          arcsinh(x)

static float asinhf(float x)

static double asinh(double x)

 

// Inverse hyperbolic cosine

// class:             SpecMath

// input:             x

// output:          return value      =          arccosh(x)

static float acoshf(float x)

static double acosh(double x)

 

// Inverse hyperbolic tangent

// class:             SpecMath

// input:             x

// output:          return value      =          arctanh(x)

static float atanhf(float x)

static double atanh(double x)

 

// Error function

// 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)

 

// Error function complement

// 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)

 

// Probit function

// 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)

 

// Natural logarithm of the gamma function

// 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)

 

// Regularized gamma function

// 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)

 

// Regularized beta function

// 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)

 

// Bessel function of the first kind

// class:             SpecMath

// input:             x

//                      n

// output:          return value      =          Jn(x)

static float bessj(float x, int n)

 

// Fast rounding double to int conversion

// 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)

 

// Fast truncating double to int conversion

// 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)

 

// Fast truncating split into integer and fractional part

// class:             SpecMath

// input:             0 ≤ x ≤ 231-1

// output:          x                      =          fractional part

//                      return value      =          largest integer ≤ input

static int fsplit(double &x)

 

// Find root of a function

// 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                                                                        )

 


Polynomials

 

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

 

 

// Find complex roots of polynomial with real coefficients

// 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)

 

// Construct polynomial from complex roots

// 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)

 

// Symbolic addition

// 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)

 

// Symbolic subtraction

// 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)

 

// Symbolic multiplication by constant

// 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)

 

// Symbolic multiplication of two polynomials

// 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)

 

// Symbolic scaling

// 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)

 

// Symbolic substitution

// 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)

 

// Reverse

// 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)

 

// Symbolic differentiation

// 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)

 

// Symbolic integration

// 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)

 

// Create Chebyshev polynomial of the first kind

// class:             SpecMath

// input:             d                      =          degree of polynomial

// output:          c[d + 1]           =          coefficients of Chebyshev polynomial

static void chebypoly(double* c, int d)

 

// Convert Chebyshev series to polynomial

// 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)

 

// Convert polynomial to Chebyshev series

// 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)

 


Interpolation

// Linear interpolation

// 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)

 

// Fit parabola to y(x) and return extremum

// 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)

 

// Fit cubic polynomial to y(x)

// 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)

 

// Cubic interpolation

// 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)

 

// Find coefficients of Chebyshev approximation of y(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)

 


Differential Equations

// Solve using forth order Runge-Kutta method

// 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)

 


Filter Design

 

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.

 

 

// Design discrete time first order filter

// 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)

 

// Design discrete time audio equalizer filter

// 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)

 

// Convert continuous to discrete time transfer function

// 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)

 

// Design continuous time first order filter

// 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)

 

// Design continuous time second order filter

// 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)

 

// Get normalized Butterworth low pass filter parameters

// 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)

 

// Get normalized Bessel low pass filter parameters

// 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)

 

// Get normalized Chebyshev type 1 low pass filter parameters

// 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)

 

// Convert normalized to general low pass filter

// 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)

 

// Convert normalized low pass to general high pass filter

// 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)

 

// Convert normalized low pass to general band pass filter

// 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)


 

AudioAnalysis

 

Signal analysis functions with focus on audio.

Spectral Analysis

// High resolution spectral analysis

// 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   )

 

// Prepare parameters for high resolution spectral analysis

// 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)

 


Analysis by Decomposition to Arbitrary Functions

// Matching pursuit

// 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                           )

Cepstral Analysis

// Convert magnitude spectrum to real cepstrum

// 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)

 

// Convert real cepstrum to magnitude spectrum

// 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)

 

// Get Mel frequency cepstral coefficients (MFCC)

// 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  )


Linear Prediction

// ***   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);

                                                                                                                                            

// Levinson-Durbin recursion

// 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)

 

// Calculate LPC residual from original signal

// 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)

 

// Restore original signal from LPC residual using a direct form IIR filter

// 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)


// Restore original signal from LPC residual using a lattice IIR filter

// 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)

 

// Calculate line spectral frequencies (LSFs) from even order LPC coefficients

// 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)

 

// Restore original signal from even order LPC residual using LSFs with an IIR filter

// 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)

 


Adaptation

// Normalized LMS algorithm

// 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                     )

 

// Tracking demodulator

// 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                              )


Feature Extraction

// Envelope follower

// 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)

 

// Fundamental frequency detector

// 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 Europe and Japan, no corresponding US patent was found, if

//                      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)

 

// Verify fundamental frequency against high resolution spectrum

// 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)

 

// Extract harmonic spectrum and inharmonicity

// 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)

 

// Count zero crossings

// 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)

 

// Spectral flatness within a defined frequency band

// 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)

 

// Upper cutoff frequency

// 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)

 

// Lower cutoff frequency

// 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)

 

// Spectral flux

// 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)

 

// Amplitude-based attack transient detector

// 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 )

 

// Spectrum-based general transient detector

// 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)

 

// Partial tracker after McAulay/Quatieri

// 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            )


Resynthesis

// Add time-varying cosine wave to existing data matching frequencies and phases

// .. 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)

 

// Add time-varying cosine wave to existing data matching frequencies

// .. 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)


 

AudioSynth

 

 

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.

 

Wavetable Oscillator

Key features: Crossfading between different tables, pitch-dependent table switching for minimal aliasing.

 

// Create oscillator

// 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)

 

// Load wavetable with spectrum

// 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)

 

// Set phase synchronously

// class:             WaveOsc

// input:             phase               =          set wave pointer position, range: 0 (start) .. 1.0f (end)

void SetPhase(float phase)

 

// Audio and control update

// 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                                                                         )

 


Ring Modulator

Includes prefilters (zero @ fs/2), correction postfilter, and DC trap with fc = 5 Hz.

 

// Create ring modulator

// class:             RingMod

// input:             smprate            =          sample rate in Hz

RingMod(float smprate)

 

// Audio and control update

// 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)

 

Noise Generator

uniformly distributed white noise, range: -1..1, RMS amplitude: 0.577

pink noise, accuracy: +/- 0.3db (0.00045..0.45fs), RMS amplitude: 0.577

 

// Create noise generator

// class:             Noise

Noise()

 

// Set noise type

// class:             Noise

// input:             type                 =          noise type: 0 (white), 1 (pink)  

void SetType(int type)

 

 // Audio and control update

// class:             Noise

// input:             samples            =          number of samples to process

// output:          out[samples]    =          audio output

void Update(float* out, int samples)

 

Static Delay Line

Implements clickless delay time change by this sequence:

  1. crossfade delayed to original signal * zlevel within 20 ms
  2. wait new delay time
  3. crossfade back to delayed signal within 20 ms

 

// Create delay

// 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)

 

// Set delay time

// class:             Delay

// input:             dlen                 =          delay time in samples   

void SetType(int dlen)

 

// Audio and control update

// 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)

 

Time-varying Delay Line

Uses 1st order allpass interpolation.

 

// Create delay

// class:             VarDelay

// input:             maxlen             =          maximum delay time in samples

VarDelay(float maxlen)

 

// Audio and control update

// 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)

 


Sample Playback Oscillator

Playback:

  1. Call PreComp once to prepare a sample for playback.
  2. Call Update deliberately to play the sample with different parameters.

 

Loops:

Loop end position is at least 8 samples away from the sample end (recommended):

Fractional and integer length allowed, switching to unlooped playback ok.

Loop end is closer than 8 samples to the sample end (in case of given loop data):

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

 

// Create oscillator

// class:             SampleOsc

// remarks:        constructing the first instance takes some extra time as a shared table is built

SampleOsc()

 

// Prepare sample

// 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)

 

// Audio and control update

// 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  )

 


Envelope Generator

// Create envelope

// 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)

 

// Set parameter

// 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)

 

// Set preset envelope

// 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)


// Create envelope event

// 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)

 

// Control update

// 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)

 


Controlled Amplifier

// Create amplifier

// class:             Amp

Amp()

 

// Set control characteristic

// class:             Amp

// input:             curve               =          control characteristic: 0 (linear), 1 (square), 2 (quartic)

void SetType(int curve)

 

// Audio and control update

// 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)

 

Convert Audio to Phase Modulation Signal

Filter input by first order lowpass with additional zero @ fs/2 and convert it to a wrapped integer.

 

// Create converter

// class:             AudioToPM

// input:             smprate            =          sample rate in Hz

//                      fc                     =          filter cutoff frequency in Hz

AudioToPM(float smprate, float fc = 400.0f)

 

// Audio and control update

// 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)

 


Hilbert Transformer

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.

 

// Create transformer

// class:             Hilbert

Hilbert()

 

// Audio and control update

// 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)

 

First Order Lowpass Filter

// Create filter

// class:             Lowpass1

// input:             smprate            =          sample rate in Hz

//                      fmin                 =          minimum -3dB cutoff frequency (Hz)

Lowpass1(float smprate, float fmin = 5.0f)

 

// Audio and control update

// 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)

                                  


First Order Highpass Filter

// Create filter

// class:             Highpass1

// input:             smprate            =          sample rate in Hz

//                      fmin                 =          minimum -3dB cutoff frequency (Hz)

Highpass1(float smprate, float fmin = 5.0f)

 

// Audio and control update

// 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)

 


Second Order Multimode Filter                                    

// Create filter

// 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)

 

// Set filter characteristic

// class:             ChambFilter

// input:             type                 =          characteristic:   0 (lowpass), 1 (bandpass), 2 (highpass),

//                                                                                 3 (peaking), 4 (notch), other (bypass)

void SetType(int type)

 

// Audio and control update (linear version)

// 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)

 

// Audio and control update (virtual analog version)

// 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)

 


Virtual Analog Moog Low Pass Filter

Key features: Low-alias wideband nonlinearity, amplitude compression, resonance up to self-oscillation with precise tuning, improved gain evolution with increasing resonance.

 

// Create filter

// 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)

 

// Audio and control update

// 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)

 


FM Oscillator

// Create oscillator

// 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)

 

// Set phase synchronously

// class:             FMOsc

// input:             phase               =          set wave pointer position, range: 0 (start) .. 1.0f (end)

void SetPhase(float phase)

 

// Modulator audio and control update

// 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          )

 

// Carrier audio and control update

// 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)

 


Virtual Analog Oscillator

// Create oscillator

// 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)

 

// Set wave shape

// class:             VAOsc

// input:             shape               =          wave shape:     0 (sawtooth), 1 (pulse)

void SetShape(int shape)

 

// Set phase synchronously

// class:             VAOsc

// input:             phase               =          set wave pointer position, range: 0 (start) .. 1.0f (end)

void SetPhase(float phase)

 

// Audio and control update

// 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)

 


Raw Sawtooth Oscillator

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.

 

// Create oscillator

// 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)

 

// Set phase synchronously

// class:             RawSawOsc

// input:             phase               =          set wave pointer position, range: 0 (start) .. 1.0f (end)

void SetPhase(float phase)

 

// Audio and control update

// 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)


 

Neuro

 

Neural networks and other classifiers. Principal and independent component analysis.

Neurons

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)

 

// Create neurons

// 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)

 

// Get number of neurons

// class:             Neurons

// output :         return value      =          number of neurons, failed: -1

int GetNofNeurons()

 

// Get type of neurons

// class:             Neurons

// output :         return value      =          type of neurons

int GetType()

 

// Content of neurons

// 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

 


Axons

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

 

// Create axons with randomized weights

// class:             Axons

// input:             in                     =          pointer to input neurons

//                      out                   =          pointer to output neurons

Axons( Neurons* in, Neurons* out)

 

// Inititalize weights

// 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)

 

// Update output neurons

// class:             Axons

void FeedForward()                                       


// Enable batch learning

// class:             Axons

void EnableBatchLearning()

 

// Disable batch learning and update weights

// class:             Axons

void DisableBatchLearning()

 

// Back propagation step

// 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)

 

// Simple Hebbian weight update

// based on the current input / output neuron values, bias weights are ignored

// class:             Axons

// input:             lrate                 =          learning rate

void Hebb(float lrate)

 

// Update weights by generalized Hebbian algorithm

// class:             Axons

// input:             lrate                 =          learning rate

void GHA(float lrate)

 

// Read input neuron weights

// 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)

 

// Write input neuron weights

// 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)

 

// Read output neuron weights

// 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)

 

// Write output neuron weights

// 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)

 


Kohonen

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".

 

// Create Kohonen, vector quantization or k-means map

// 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)

 

// Reset network to state after construction

// class:             Kohonen

void Reset()                                       

 

// Get output node corresponding to input vector

// 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)             

 

// Get vector of output node

// class:             Kohonen

// input:             node                =          output node index

// output:          d[vsize]            =          vector of specified output node

void GetVector(float* d, int node)

 

// Learn input vector

// 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)

 

// Perform k-means clustering

// 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)

 


PCA

Principal and independent component analysis.

Component properties: orthonormal, zero mean, sign undetermined

 

// Get principal components

// 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)

 

// Decompose vector into sum of components

// … 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 )

 

// Compose vector as sum of components

// … 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 )

           


// Convert principal to independent components

// … 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 ica(float* ics, float* pcs, int vsize, int m, int n)


 

AudioFile

 

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


// Create audio file object

// class:             AudioFile

AudioFile()

 

// Create new audio file

// 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                                                         )

 

// Load file

// 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)

 

// Load raw file and specify format

// 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                                          )


// Save as WAVE file

// 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)

 

// Append to existing WAVE file

// class:             AudioFile

// input:             filename           =          full pathname including extension

// output:          return value      =          error code (see above)

int AppendWave(char *filename)

 

// Get size

// class:             AudioFile

// output:          return value      =          size in frames

unsigned int GetSize()

                                  

// Get sample rate

// class:             AudioFile

// output:          return value      =          sample rate in Hz

unsigned int GetRate()

 

// Get number of channels

// class:             AudioFile

// output:          return value      =          number of channels

unsigned int GetChannels()                             

 

// Get resolution

// class:             AudioFile

// output:          return value      =          resolution in bits

unsigned int GetResolution()                

 

// Get speaker positions

// class:             AudioFile

// output:          return value      =          speaker position (see above)

unsigned int GetSpkPos()

 

// Get safe pointer to audio data

// 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)


 

Chart

 

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)!

 

// Create chart

// class:             Chart

Chart()

 

// Initialize 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     )


// Draw chart

// 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          )

 

// Plot line diagram

// 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)

 

// Plot point diagram

// 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)

 

// Plot filled circle diagram

// 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)

 

// Plot vertical bulk diagram

// 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)

 

// Plot vertical bulk interval diagram

// 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)

 

// Plot horizontal bulk interval diagram

// 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)

 

// Clear diagram

// class:             Chart

// input:             hidden              =          true:     clear memory, no chart update

//                                                         false:    clear chart and memory

void Clear(bool hidden = false)

 

// Draw diagram frame

// 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)

 

// Highlight diagram zone

// 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 diagram segment

// 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)

 

// Unelect diagram segment

// class:             Chart

// input:             hidden              =          true:     draw to memory, no chart update

//                                                         false:    draw and update chart from memory

void Unselect(  bool hidden = false)

 

// Get selection

// 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)

 

// Enable mouse response

// mouse messages are sent on mouse moves and clicks inside the diagram

// class:             Chart

void EnableMouseResponse()

 

// Disable mouse response                          

// class:             Chart

void DisableMouseResponse()                       

 

// Enable selection by mouse

// 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)

 

// Disable selection by mouse                     

// class:             Chart

void DisableMouseSelection()

 

// Get x value from mouse message

// class:             Chart

// input:             lParam             =          mouse message

// output:          return value      =          diagram x value associated with mouse message

float MsgToX(unsigned int lParam);                

 

// Get y value from mouse message

// class:             Chart

// input:             lParam             =          mouse message

// output:          return value      =          diagram y value associated with mouse message

float MsgToY(unsigned int lParam)

 

// Get right button activity from mouse message

// class:             Chart

// input:             lParam             =          mouse message

// output:          return value      =          true:     right button active

bool RightButton(unsigned int lParam)

 

// Get left button activity from mouse message

// class:             Chart

// input:             lParam             =          mouse message

// output:          return value      =          true:     left button active

bool LeftButton(unsigned int lParam)

 

// Get from mouse message: Has button just been pressed?

// class:             Chart

// input:             lParam             =          mouse message

// output:          return value      =          true:     button has just been pressed

bool Pressed(unsigned int lParam)

 

// Get from mouse message: Has button just been released?

// class:             Chart

// input:             lParam             =          mouse message

// output:          return value      =          true:     button has just been released

bool Released(unsigned int lParam)

 

// Get safe pointer to graphic device context

// useful to draw anything into the diagram, see also GetDiagramFrame()

// class:             Chart

// output:          return value      =          pointer to device context (NULL if uninitialized)

CDC* GetSafeDeviceContext()

 

// Get diagram frame

// useful to draw anything into the diagram, see also GetSafeDeviceContext()

// class:             Chart

// output:          return value      =          diagram frame (0 if uninitialized)

CRect GetDiagramFrame()