All Current Processors

All SFPL Processors As Of 27 November 2012

This is constant work in progress - I am working my way down the list as I get chance - and updating as stuff gets added/updated. Please note that some times the descriptions says 'change the volume...' for example. This is just to make reading easier; SFPL has 'immutable' types so it is not possible to change the volume of a signal or trim white space off a string. What actually happens is that a new one is made and the old one left unchanged.

December 2013 - this is quite out of data, I will work on it over the next few days...

• MakeSawTooth
  Convert a signal into a saw tooth using center crossing
• MakeSquare
  Convert a signal into a square wave using center crossing
• SubOctave
  Make a sin wave the period of which is half the rate of center crossing.
• MakeTriangle
  Make a signal into a triangle wave using centre crossing
• Saturate
  Distortion effect using 1-(1/x) law
• ShapedThreshold
  Convert a signal into a pulse wave using a threshold based on a shape
• Threshold
  Convert a signal into a pulse wave using a fixed threshold
• WaveShaper
  Distort a wave using a 6th order polynomial transfer function
• ClipToSafe
  Use DC removal and fadeout to attempt to get rid of clicks entering and leaving a audio signal
• Concatenate
  Join two or more audio signals end to end 
• Convolve
  Convolve one signal with another. Convolution is forwards and backwards, i.e. the convolving signals centre is placed at the current sample on the main signal.
• ConvolveDelay
  Mix a signal with a delayed version of its self where the delayed version is actually mixed with the original via convolution via a second convolving signal.
• ConvolveForward
  Convolve one signal with another. Convolution is forwards only (see Convolve).
• CrossFadeSplit
  Create a bunch of subsections of a signal where the sub sections are overlapping and are pre-faded.
• Cut
  Take a subsection of a signal.
• Divide
  Sample my sample division of one signal by another.
• Mix
  Mix two or more signals together using direct sample addition.
• MixAt
  Mix two or more signals together using direct sample addition and with defined offsets for each.
• Multiply
  Sample by sample multiplication of one signal by another.
• GaussianShape
  Generate a signal which follows a Gaussian shape.
• HammingShape
  Generate a signal which follows a Hamming shape.
• HannShape
  Generate a signal which follows a Hann shape.
• Note
  Take a string representation of a note in equal temperament and forward the frequency in Hz based on A4 being 440Hz. It uses a post fix of b or # for flat and sharp. E.g. "C5#" "E2b". 
• NumericShape
  Takes a bunch of time and magnitude values and creates a signal from the linear interpolation of the values.
• PhasedSinWave
  Create a sin wave of a given length and frequency and starting at a given phase.
• ShapedSinWave
  Create a sin wave the frequency of which is continuously set by a shape signal. This effectively is an infinitely variable sin wave generator and can create things like FM.
• Silence
  Generate a given length of silence.
• SimpleShape
  Generate a shape signal from a series of decible and time values.
• SinWave
  Generate a sin wave of a given length and frequency.
• SincShape
  Generate a Sinc shape signal with a given period and width.
• Slide
  Generate a sin wave which continuously in frequency varies based on a series of time and frequency values.
• WhiteNoise
  Generate pure white noise of a given length.
• ToJSON
  Convert a SFPL data type (signal, number, bunch etc) into a JSON string.
• CreateBesselBandPass
  Create a Bessel based band pass filter which can be passed to another processor.
• CreateButterworthBandPass
  Create a Butterworth based band pass filter which can be passed to another processor.
• CreateButterworthHighPass
  Create a Butterworth based high pass filter which can be passed to another processor.
• CreateButterworthLowPass
  Create a Butterworth based low pass filter which can be passed to another processor.
• DirectRelength
  Change the length of an audio signal by a fixed ratio via decimation and / or interpolation resulting in a new sample of different length and changing all the frequencies.
• DirectResample
  The same as DirectRelength except the resulting signal is cut or padded with silence to result in it being the same length as the original.
• FrequencyModulate
  Use an shape signal to continuously resample another signal so as to modulate its frequency. This is different from Resample in that the change in frequency is symmetric about zero for the modulating waveform. A positive value of the modulator will increase the frequency of the modulated and a negative modulator value decreases the frequency of the modulated. This method is centre frequency stable for stable modulation frequencies. The centre frequency of the modulate signal will shift if the frequency of the modulator signal shifts overall throughout the modulation period. For complex waveforms, things are more complex.
• ShapedLadderFilter
  A classic synthesiser low pass filter. This is loosely similar to something like the low pass on a  minimoog but it is not as fruity. It acts as a low pass with variable cut off and variable resonance.
• ButterwothBandPass
  Use a Butterworth filter for band pass filtering up to 5 poles.
• BesselBandPass
  Use a Bessel filter for band pass filtering up to 5 poles.
• ButterworthHighPass
  Use a Butterworth filter for high pass filtering up to 11 poles.
• BesselHighPass
  Use a Bessel filter for high pass filterring up to 11 poles.
• ButterworthLowPass
  Use a Butterworth filter for low pass filtering up to 11 poles.
• BesselLowPass
  Use a Bessel filter for low pass filtering up to 11 poles.
• RBJLowPass
  Use a Robert Bristow-Johnson filter for 2 pole low pass variable Q filtering.
• RBJHighPass
  Use a Robert Bristow-Johnson filter for 2 pole high pass variable Q filtering.
• RBJBandPass
  Use a Robert Bristow-Johnson filter for 2 pole band pass filtering with variable Q.
• RBJPeaking
  Use a Robert Bristow-Johnson filter for 2 pole filter as a peaking resonator with variable Q and gain.
• RBJBandReject
  Use a Robert Bristow-Johnson filter for 2 pole notch filter with variable Q.
• RBJLowShelf
  Use a Robert Bristow-Johnson filter for 2 pole low shelf filter with variable Q.
• RBJHighShelf
  Use a Robert Bristow-Johnson filter for 2 pole high shelf filter with variable Q.
• Resample
  Resample a signal at a rate determined by a second signal. Unlike FrequencyModulate, the next sample rate is simply the modulation signal value. 
• ShapedButterworthBandPass
  Use a Butterworth filter as a band pass filter up to 5 poles but make the low and high frequency shoulders controlled by two other signals allowing continuously variable filter width.
• Mixer
  Chooses a Java Audio mixer given a description.
• LineWait
  Waits for a Java Audio line to finish.
• Mixers
  Produces a bunch of description of the available Java Audio mixers.
• Monitor
  Plays a single audio signal on the default output device using Java Audio.
• PlayFile
  Plays the contents of a file on the default output device using Java Audio.
• ReadFile
  Reads a wav file into a bunch of signals. Forwards a bunch with one signal per audio channel in the file.
• Semitone
  Forwards a number representing the twelfth root of 2.
• SpeedOfSound
  Forwards a number giving the approximate speed of sound at sea level in meters per second.
• StereoMonitor
  Plays a stereo sound on the default audio device using Java Audio. Takes a bunch of two signals and plays them.
• WriteFile16
  Take a bunch of audio signals and write them out to a 16 bit wav file with one channel per element in the bunch.
• WriteFile32
  Take a bunch of audio signals and write them out to a 16 bit wav file with one channel per element in the bunch.
• WriteFileString
  Write a string to a file.
• DirectResonate
• FilterResonate
• ResonantFilter
  Filter via feedback resonance. Note that this is nothing at all like an IIR peaking resonator. 
• Resonate
• FilteredResonantFilter
  Filter via feedback resonance with an IIR filter in the feedback loop (see ResonantFilter). 
• Reverse
  Reverse the direction of an audio signal.
• ShapedResonantFilter
  Filter via feedback resonance where the feedback delay time is varied according to the value of a shaping signal. The shaping signal alters the delay in a simple multiplication way; e.g. a value of 2 will double the delay, a value of 0.5 will halve it (see ResonantFilter). 
• Frequency
  Convert a period for a single cycle in milliseconds into a frequency in Hz.
• Length
  Forward the length of the passed audio signal in milliseconds.
• Period
  Convert a frequency in Hz into a period for a single cycle in milliseconds.
• ValueAt
  Forward the numeric value of a signal at a point in time in milliseconds.
• Normalise
  Remove DC from an audio signal and make set the magnitude so the maximum excursion from zero is 1.
• NormaliseArea
  Scale an audio signal so that the area between the signal and zero sums to 1.
• RemoveDC
  Remove an DC component from an audio signal.
• ShapedPower
  Raise all the values of a signal to the power of an equivalent length shaping signal. E.g. if the shaping signal as a value of 2 at 10 milliseconds, then the value of the forwarded audio signal at 10 milliseconds will be the square of the input signal at that point.
• DoneAll
  Perform Done on a bunch of Do Tasks in order and forward a bunch holding the results.
• DoAll
  Do a bunch of code blocks, then accumulate the Done result in a bunch in order of the code blocks in to input bunch, Forward the resulting in a bunch of the same order.
• Max
  Forward the greatest (closest to positive infinity) numerical value from a bunch of numbers.
• Min
  Forward the lowest (closest to negative infinity) numeric value from a bunch of numbers.
• Prime
  Forward the closest prim number which is larger than the passed number.
• Random
  Forward a random number between 0 and 1.
• Truncate
  Forward a number with any decimal part removed. E.g. 1.57 will forward to 1.
• Clip
  Set any value in an audio signal which is above 1 to 1 and any below -1 to -1.
• DirectMix
  Add a constant value to every value in an audio signal and forward the result.
• Invert
  Invert an audio signal about zero.
• Power
  Raise every value in an audio signal to a index.
• Rectify
  Convert all negative values in an audio signal to positive values of the same magnitude.
• Volume
  Scale the volume of an audio signal by a value given in db.
• NumericVolume
  Scale the volume of an audio signal by a value given as a linear number (2 will double etc).
• FromDBs
  Convert a number passed in dbs to a linear number and forward that number.
• ToDBs
  The oposite of FromDBs.
• MaxValue
  Find the largest (closest to positive infinity) value of an audio signal.
• PrintLn
  Send a string to standard out and append a new line.
• |
  Forward the modulus of two numbers.
• -
  Forward the different of two numbers.
• /
  Forward one number divided by another.
• GT
  Forward TRUE if the first passed number is greater than the second, otherwise forward false.
• LT
  Forward TRUE if the first passed number is less than the second, otherwise forward false.
• EQ
  Forward TRUE if the first passed number is equal to the second, otherwise forward false.
• NOT
  If passed TRUE forward FALSE, if passed FALSE forward TRUE.
• OR
  Forward the Boolean OR of the passed bunch of Booleans.
• AND
  Forward the Boolean AND of the passed bunch of Booleans.
• XOR
  Forward the Boolean XOR of the bassed bunch of Booleans.
• StrCat
  Concatenate all the members of the passed bunch of string and forward the result.
• StrIndex
  Forward the index (zero based) of a string in another string. -1 if not found.
• String
  Convert any operand into a string representation of its self.
• StrLen
  Forward the length of the passed string.
• StrToUpper
  Forward the a string based on the passed string where all the characters are turned to upper case using the current locale.
• StrToLower
  Forward the a string based on the passed string where all the characters are turned to upper case using the current locale.
• StrTrim
  Forward a string with all the white space at the start and end of the passed string removed.
• Integer
  Round a number to the nearest integer using the Java default rounding rules.
• Follow
  Create an envelop signal from the magnitude of the passed audio signal with settable attack and release rates.
• Console
  Cause patch execution to pause the move over to interactive console working.
• Granulate
  Split an audio signal into tapered (fade in and fade out) grains of a given length with an optionally given random variation up from that length.
• AutoCorrelate
  Create a signal which is the direct auto-correlation of the passed signal.
• DBs-100 DBs-99 DBs-98 ... DBs+98 DBs+99 DBs-100
  Change the volume of the passed audio signal by an amount in dbs and forward the result.
• Pcnt-100 Pcnt-99 Pcnt-98 ... Pcnt+98 Pcnt+99 Pcnt+100
  Change the volume of the passed audio signal by an amount in linear numbers as a percent and forward the result.

Sonic Field 0.2 Released

Sonic Field 0.2 Has Been Released Under AGPL 3.0

The .zip of all the source and class files is here: http://www.sonicfield.co.uk/downloads/SonicField-0-2.zip

This download contains all the dependancy jars and the information of the licensing for them as well.

From the VERSION.txt file:


Version 0.2
===========

This version of Sonic Field has many enhancements including but not limited to:
1)  Granular synthesis
2)  RBJ IIR Filers (high,low,high-shelf,low-shelf,peaking)
3)  Cleaner implementation of resonance synthesis
4)  Better error handling - fix the line number bug
5)  System hang on error in paralleliser bug fixed
6)  Envelope follower.
7)  Fix to numerous bugs in shaped
8)  Addition of shaped butter worth and bessel bandpass filters
9)  Improved FM signal modulation and addition of fm signal generator
10) Addition of shaped power

Site Update

Just added two tracks to the examples page

Eat Bass and megaphone-1 are now available form this examples page.

Current:

Sonic Field Of Dreams: 
Dreams: A growing collection of perforances exploring the generation of digital sounds using inspiration from the analogue world.
Synthesis And Processing Examples: 
Eat Bass: A demonstration of granular synthesis for bass effects and pitch shifting.
Megaphone-1: A demonstration of granular synthesis, envelope following and eveloped filtering to make a stylophone sound like a church organ.
Early Stuff: 
Somme: 4 ab-initio synthesised saxophones with abinito rain and sampled drums. This piece also has fully spatialised reverberation.   
Romance Electric: 2 ab-initio clean tone electric guitars playing Romance.   
Romance Semi-Acoustic: 4 ab-initio semi-acoustic guitars playing Romance.   
Theme Of Four: An ab-initio rendering of hand played drums with variations around 4/4 time.   

Clips And Effects

Flute Dry A classical harp and flute clip unprocessed.
Flute Reverb' The flute clip with spatialised reverberation added.   

Granulate

Granulate - the gateway to granular synthesis

Granulate takes a signal and breaks it up into a series of sound grains. The over all length of the grain is passed into granulate:


(?signal,50)granulate !new

In the above, signal is broken up into 50ms grains. The grains are faded in and out so that the first 25 milliseconds (in the above) are fading in and the second are fading out.

The result of Granulate is a bunch of bunches. Each inner bunch consists of a sample as the grain and the time in milliseconds from which it came. Passing the output of Granulate into MixAt will exactly reproduce exactly the original signal appart from the first half and last half of a grain which are missing causing a short fade in and fade out section.

Here is a test case for Granulate:

[ Granulates a singal so that the middle section when reconstructed and mixed
  with the orignal should be within rounding error of silent
]
(10000,440)SinWave !signal
(?signal,50)granulate !new
>new MixAt Normalise !dog
(>dog,>signal Invert)Mix !dog
0 !zero
(
   1000,2000,
   {
      !x
      (>zero,((?dog,>x)ValueAt,2)**)+ !zero
   }
)Repeat
(?zero,0.1)LT Println

Implementing RBJ Filters In Java

A flare for simplicity.
Image From NASA Solar Observatory.

Infinite Impulse Response filters a super useful in many application but especially so for audio work where very large numbers of samples need to be filtered. 

IIR filters are sometimes called recursive filters. They work by feeding some of their signal back into the incoming signal. This means that even though for each sample of the incoming signal they only need to make a few calculations, they implicitly use the entire length of preceding signal. The limit to this model is only the numerical accuracy of the mathematics being used.

There is something of a down side to the seemingly magical properties of IIR filters - the mathematics used to create what are called the 'coefficients' is a bit tricky. I guess, given enough time, I could get my head around it, but I have other things to do and I doubt I would be that good at it anyhow! What we need is someone very much better at maths than I am who can make a simple cookbook of how to compute the coefficients for a set of general purpose filters which sit in the sweet spot between not filtering enough and going unstable.

Most pipe dreams do not come true but this one seems to have done so! Robert Bristow-Johnson created just such a cook book. It is called the Audio EQ Cookbook and it does just what we want. It shows how to compute the coefficients for the 'bi-quad' variants for IIR filters for all the handy filters we normally need.

package com.nerdscentral.audio.pitch.algorithm;

import com.nerdscentral.audio.SFConstants;

/*
 * RBJ Filters from C++ version by arguru[AT]smartelectronix[DOT]com based on eq filter cookbook by Robert Bristow-Johnson
 * <rbj@audioimagination.com> This code is believed to be public domain and license free after best efforts to establish its
 * licensing.
 */
public class SFRBJFilter
{

    // filter coeffs
    double b0a0, b1a0, b2a0, a1a0, a2a0;

    // in/out history
    double ou1, ou2, in1, in2;

    public SFRBJFilter()
    {
        // reset filter coeffs
        b0a0 = b1a0 = b2a0 = a1a0 = a2a0 = 0.0;

        // reset in/out history
        ou1 = ou2 = in1 = in2 = 0.0f;
    }

    public double filter(double in0)
    {
        // filter
        final double yn = b0a0 * in0 + b1a0 * in1 + b2a0 * in2 - a1a0 * ou1 - a2a0 * ou2;

        // push in/out buffers
        in2 = in1;
        in1 = in0;
        ou2 = ou1;
        ou1 = yn;

        // return output
        return yn;
    }

    public enum FilterType
    {
        LOWPASS, HIGHPASS, BANDPASS_SKIRT, BANDPASS_PEAK, NOTCH, ALLPASS, PEAK, LOWSHELF, HIGHSHELF
    }

    public void calc_filter_coeffs(final FilterType type, final double frequency, final double q, final double db_gain)
    {
        boolean q_is_bandwidth;
        final double sample_rate = SFConstants.SAMPLE_RATE;
        switch (type)
        {
        case ALLPASS:
        case HIGHPASS:
        case LOWPASS:
        case LOWSHELF:
        case HIGHSHELF:
            q_is_bandwidth = false;
            break;
        default:
            q_is_bandwidth = true;
            break;
        }
        // temp pi
        final double temp_pi = 3.1415926535897932384626433832795;

        // temp coef vars
        double alpha, a0 = 0, a1 = 0, a2 = 0, b0 = 0, b1 = 0, b2 = 0;

        // peaking, lowshelf and hishelf
        if (type == FilterType.PEAK || type == FilterType.HIGHSHELF || type == FilterType.LOWSHELF)
        {
            final double A = Math.pow(10.0, (db_gain / 40.0));
            final double omega = 2.0 * temp_pi * frequency / sample_rate;
            final double tsin = Math.sin(omega);
            final double tcos = Math.cos(omega);

            if (q_is_bandwidth) alpha = tsin * Math.sinh(Math.log(2.0) / 2.0 * q * omega / tsin);
            else
                alpha = tsin / (2.0 * q);

            final double beta = Math.sqrt(A) / q;

            // peaking
            if (type == FilterType.PEAK)
            {
                b0 = (1.0 + alpha * A);
                b1 = (-2.0 * tcos);
                b2 = (1.0 - alpha * A);
                a0 = (1.0 + alpha / A);
                a1 = (-2.0 * tcos);
                a2 = (1.0 - alpha / A);
            }

            // lowshelf
            if (type == FilterType.LOWSHELF)
            {
                b0 = (A * ((A + 1.0) - (A - 1.0) * tcos + beta * tsin));
                b1 = (2.0 * A * ((A - 1.0) - (A + 1.0) * tcos));
                b2 = (A * ((A + 1.0) - (A - 1.0) * tcos - beta * tsin));
                a0 = ((A + 1.0) + (A - 1.0) * tcos + beta * tsin);
                a1 = (-2.0 * ((A - 1.0) + (A + 1.0) * tcos));
                a2 = ((A + 1.0) + (A - 1.0) * tcos - beta * tsin);
            }

            // hishelf
            if (type == FilterType.HIGHSHELF)
            {
                b0 = (A * ((A + 1.0) + (A - 1.0) * tcos + beta * tsin));
                b1 = (-2.0 * A * ((A - 1.0) + (A + 1.0) * tcos));
                b2 = (A * ((A + 1.0) + (A - 1.0) * tcos - beta * tsin));
                a0 = ((A + 1.0) - (A - 1.0) * tcos + beta * tsin);
                a1 = (2.0 * ((A - 1.0) - (A + 1.0) * tcos));
                a2 = ((A + 1.0) - (A - 1.0) * tcos - beta * tsin);
            }
        }
        else
        {
            // other filters
            final double omega = 2.0 * temp_pi * frequency / sample_rate;
            final double tsin = Math.sin(omega);
            final double tcos = Math.cos(omega);

            if (q_is_bandwidth) alpha = tsin * Math.sinh(Math.log(2.0) / 2.0 * q * omega / tsin);
            else
                alpha = tsin / (2.0 * q);

            // lowpass
            if (type == FilterType.LOWPASS)
            {
                b0 = (1.0 - tcos) / 2.0;
                b1 = 1.0 - tcos;
                b2 = (1.0 - tcos) / 2.0;
                a0 = 1.0 + alpha;
                a1 = -2.0 * tcos;
                a2 = 1.0 - alpha;
            }

            // hipass
            if (type == FilterType.HIGHPASS)
            {
                b0 = (1.0 + tcos) / 2.0;
                b1 = -(1.0 + tcos);
                b2 = (1.0 + tcos) / 2.0;
                a0 = 1.0 + alpha;
                a1 = -2.0 * tcos;
                a2 = 1.0 - alpha;
            }

            // bandpass csg
            if (type == FilterType.BANDPASS_SKIRT)
            {
                b0 = tsin / 2.0;
                b1 = 0.0;
                b2 = -tsin / 2;
                a0 = 1.0 + alpha;
                a1 = -2.0 * tcos;
                a2 = 1.0 - alpha;
            }

            // bandpass czpg
            if (type == FilterType.BANDPASS_PEAK)
            {
                b0 = alpha;
                b1 = 0.0;
                b2 = -alpha;
                a0 = 1.0 + alpha;
                a1 = -2.0 * tcos;
                a2 = 1.0 - alpha;
            }

            // notch
            if (type == FilterType.NOTCH)
            {
                b0 = 1.0;
                b1 = -2.0 * tcos;
                b2 = 1.0;
                a0 = 1.0 + alpha;
                a1 = -2.0 * tcos;
                a2 = 1.0 - alpha;
            }

            // allpass
            if (type == FilterType.ALLPASS)
            {
                b0 = 1.0 - alpha;
                b1 = -2.0 * tcos;
                b2 = 1.0 + alpha;
                a0 = 1.0 + alpha;
                a1 = -2.0 * tcos;
                a2 = 1.0 - alpha;
            }
        }

        // set filter coeffs
        b0a0 = (b0 / a0);
        b1a0 = (b1 / a0);
        b2a0 = (b2 / a0);
        a1a0 = (a1 / a0);
        a2a0 = (a2 / a0);
    }
}

BiquadPeak

A Biquadratic Peaking Filter

(?signal,?f1,0.5,48)BiquadPeak  Normalise  !signal1

(?signal1,?f1,0.5,48)BiquadPeak Normalise  !signal1

(?signal1,?f1,0.2,48)BiquadPeak Normalise  !signal1

The above will make a tone like sounds form white noise.

(?signal,?centre-frequency,?q,?gain)BiquadPeak 

BiquadPeak is an IIR filter which filters just a peak. It is different form BandPass in that BandPass produces clean linearly slowed sides to the pass bands when graphed in DBs. Peak produces a sigmoidal side. The peak width can be very narrow or very broad, this is controlled by q. ?gain controls the amount the top/center of the peak is in increased amplitude compared to signal at a arbitrary distance from the that centre. The centre is naturally controlled bu centre-frequency. 

Amazingly Handy Document - RBJ EQ CookBook

A Huge thanks to Robert Bristow-Johnson for having written released this VERY useful document!

         Cookbook formulae for audio EQ biquad filter coefficients
----------------------------------------------------------------------------
           by Robert Bristow-Johnson  <rbj@audioimagination.com>


All filter transfer functions were derived from analog prototypes (that
are shown below for each EQ filter type) and had been digitized using the
Bilinear Transform.  BLT frequency warping has been taken into account for
both significant frequency relocation (this is the normal "prewarping" that
is necessary when using the BLT) and for bandwidth readjustment (since the
bandwidth is compressed when mapped from analog to digital using the BLT).

First, given a biquad transfer function defined as:

            b0 + b1*z^-1 + b2*z^-2
    H(z) = ------------------------                                  (Eq 1)
            a0 + a1*z^-1 + a2*z^-2

This shows 6 coefficients instead of 5 so, depending on your architechture,
you will likely normalize a0 to be 1 and perhaps also b0 to 1 (and collect
that into an overall gain coefficient).  Then your transfer function would
look like:

            (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
    H(z) = ---------------------------------------                   (Eq 2)
               1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

or

                      1 + (b1/b0)*z^-1 + (b2/b0)*z^-2
    H(z) = (b0/a0) * ---------------------------------               (Eq 3)
                      1 + (a1/a0)*z^-1 + (a2/a0)*z^-2


The most straight forward implementation would be the "Direct Form 1"
(Eq 2):

    y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]
                        - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]            (Eq 4)

This is probably both the best and the easiest method to implement in the
56K and other fixed-point or floating-point architechtures with a double
wide accumulator.



Begin with these user defined parameters:

    Fs (the sampling frequency)

    f0 ("wherever it's happenin', man."  Center Frequency or
        Corner Frequency, or shelf midpoint frequency, depending
        on which filter type.  The "significant frequency".)

    dBgain (used only for peaking and shelving filters)

    Q (the EE kind of definition, except for peakingEQ in which A*Q is
        the classic EE Q.  That adjustment in definition was made so that
        a boost of N dB followed by a cut of N dB for identical Q and
        f0/Fs results in a precisely flat unity gain filter or "wire".)

     _or_ BW, the bandwidth in octaves (between -3 dB frequencies for BPF
        and notch or between midpoint (dBgain/2) gain frequencies for
        peaking EQ)

     _or_ S, a "shelf slope" parameter (for shelving EQ only).  When S = 1,
        the shelf slope is as steep as it can be and remain monotonically
        increasing or decreasing gain with frequency.  The shelf slope, in
        dB/octave, remains proportional to S for all other values for a
        fixed f0/Fs and dBgain.



Then compute a few intermediate variables:

    A  = sqrt( 10^(dBgain/20) )
       =       10^(dBgain/40)     (for peaking and shelving EQ filters only)

    w0 = 2*pi*f0/Fs

    cos(w0)
    sin(w0)

    alpha = sin(w0)/(2*Q)                                       (case: Q)
          = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) )           (case: BW)
          = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 )         (case: S)

        FYI: The relationship between bandwidth and Q is
             1/Q = 2*sinh(ln(2)/2*BW*w0/sin(w0))     (digital filter w BLT)
        or   1/Q = 2*sinh(ln(2)/2*BW)             (analog filter prototype)

        The relationship between shelf slope and Q is
             1/Q = sqrt((A + 1/A)*(1/S - 1) + 2)

    2*sqrt(A)*alpha  =  sin(w0) * sqrt( (A^2 + 1)*(1/S - 1) + 2*A )
        is a handy intermediate variable for shelving EQ filters.


Finally, compute the coefficients for whichever filter type you want:
   (The analog prototypes, H(s), are shown for each filter
        type for normalized frequency.)


LPF:        H(s) = 1 / (s^2 + s/Q + 1)

            b0 =  (1 - cos(w0))/2
            b1 =   1 - cos(w0)
            b2 =  (1 - cos(w0))/2
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha



HPF:        H(s) = s^2 / (s^2 + s/Q + 1)

            b0 =  (1 + cos(w0))/2
            b1 = -(1 + cos(w0))
            b2 =  (1 + cos(w0))/2
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha



BPF:        H(s) = s / (s^2 + s/Q + 1)  (constant skirt gain, peak gain = Q)

            b0 =   sin(w0)/2  =   Q*alpha
            b1 =   0
            b2 =  -sin(w0)/2  =  -Q*alpha
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha


BPF:        H(s) = (s/Q) / (s^2 + s/Q + 1)      (constant 0 dB peak gain)

            b0 =   alpha
            b1 =   0
            b2 =  -alpha
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha



notch:      H(s) = (s^2 + 1) / (s^2 + s/Q + 1)

            b0 =   1
            b1 =  -2*cos(w0)
            b2 =   1
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha



APF:        H(s) = (s^2 - s/Q + 1) / (s^2 + s/Q + 1)

            b0 =   1 - alpha
            b1 =  -2*cos(w0)
            b2 =   1 + alpha
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha



peakingEQ:  H(s) = (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1)

            b0 =   1 + alpha*A
            b1 =  -2*cos(w0)
            b2 =   1 - alpha*A
            a0 =   1 + alpha/A
            a1 =  -2*cos(w0)
            a2 =   1 - alpha/A



lowShelf: H(s) = A * (s^2 + (sqrt(A)/Q)*s + A)/(A*s^2 + (sqrt(A)/Q)*s + 1)

            b0 =    A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
            b1 =  2*A*( (A-1) - (A+1)*cos(w0)                   )
            b2 =    A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
            a0 =        (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
            a1 =   -2*( (A-1) + (A+1)*cos(w0)                   )
            a2 =        (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha



highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)

            b0 =    A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha )
            b1 = -2*A*( (A-1) + (A+1)*cos(w0)                   )
            b2 =    A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha )
            a0 =        (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha
            a1 =    2*( (A-1) - (A+1)*cos(w0)                   )
            a2 =        (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha





FYI:   The bilinear transform (with compensation for frequency warping)
substitutes:

                                  1         1 - z^-1
      (normalized)   s  <--  ----------- * ----------
                              tan(w0/2)     1 + z^-1

   and makes use of these trig identities:

                     sin(w0)                               1 - cos(w0)
      tan(w0/2) = -------------           (tan(w0/2))^2 = -------------
                   1 + cos(w0)                             1 + cos(w0)


   resulting in these substitutions:


                 1 + cos(w0)     1 + 2*z^-1 + z^-2
      1    <--  ------------- * -------------------
                 1 + cos(w0)     1 + 2*z^-1 + z^-2


                 1 + cos(w0)     1 - z^-1
      s    <--  ------------- * ----------
                   sin(w0)       1 + z^-1

                                      1 + cos(w0)     1         -  z^-2
                                  =  ------------- * -------------------
                                        sin(w0)       1 + 2*z^-1 + z^-2


                 1 + cos(w0)     1 - 2*z^-1 + z^-2
      s^2  <--  ------------- * -------------------
                 1 - cos(w0)     1 + 2*z^-1 + z^-2


   The factor:

                    1 + cos(w0)
                -------------------
                 1 + 2*z^-1 + z^-2

   is common to all terms in both numerator and denominator, can be factored
   out, and thus be left out in the substitutions above resulting in:


                 1 + 2*z^-1 + z^-2
      1    <--  -------------------
                    1 + cos(w0)


                 1         -  z^-2
      s    <--  -------------------
                      sin(w0)


                 1 - 2*z^-1 + z^-2
      s^2  <--  -------------------
                    1 - cos(w0)


   In addition, all terms, numerator and denominator, can be multiplied by a
   common (sin(w0))^2 factor, finally resulting in these substitutions:


      1         <--   (1 + 2*z^-1 + z^-2) * (1 - cos(w0))

      s         <--   (1         -  z^-2) * sin(w0)

      s^2       <--   (1 - 2*z^-1 + z^-2) * (1 + cos(w0))

      1 + s^2   <--   2 * (1 - 2*cos(w0)*z^-1 + z^-2)


   The biquad coefficient formulae above come out after a little
   simplification.