Filters

Experienced DSP engineers may find also this document useful.

Contents:

General notes
Filter design and analysis
Filtering
Numerical problems
IIR filters
FIR filters
Adaptive filters
Miscellaneous
1. Median filter
2. Wiener filter
64-bit filters

General notes

In the most general sense, the filter (system), is basically anything that somehow processes an input signal and produces an output signal. Hence, the high-level interfaces for filters are very generic:

IFilter:
- ApplyTo(signal)
IOnlineFilter:
- Process(sample)
- Reset()
additional methods (IFilterExtensions):
- Process(inBuffer, outBuffer, length, startPos, endPos)
- Process(signal)

IFilter interface defines the contract for offline filter, whereas IOnlineFilter defines the contract for online filter. The term 'offline' means that the filtering routine accepts entire signal, processes it and returns the entire signal as well (ApplyTo() method). Online filters, on the other hand, allow processing sample after sample (returning each new sample immediately). Usually they maintain some state under the hood, and this state can be reset (Reset() method).

Ideally, all actual filters would implement both interfaces, and ApplyTo() method would simply reduce to calling / inlining the Process() method sample after sample. And this is how most filters are implemented in NWaves. Yet, there are situations when alternative implementation of ApplyTo() is better in terms of performance.

Filters broadly fall into two categories:

linear filters
non-linear filters

The former category is ubiquitous and theoretically well-studied. The latter is more difficult to analyze, and only some examples are implemented in NWaves: median filter, overdrive and distortion audioeffects.

Linearity of a filter means that it produces the linear combination of output signals in response to a linear combination of input signals, and the form of this combination doesn't change. If the useful property of time invariance also holds (i.e. time shift of input signal leads to the same time shift of the output signal), then linear filter is called an LTI filter.

LTI filters are represented with abstract class LtiFilter. It inherits both interfaces IFilter / IOnlineFilter and adds an abstract property for transfer function. The transfer function uniquely and completely describes any LTI filter.

Filter design and analysis (FDA)

Transfer function and difference equation

Transfer function is the representaion of an LTI filter in z-domain, whereas difference equation is analogous representation of the filter in time-domain.

Difference equations look like this:

Coefficients b[k] define the non-recursive part of the filter and a[m] are recursive coefficients.

This is the most general form of an LTI filter. If we feed the unit impulse to it, the resulting signal (i.e. impulse response) will be infinite. That's why such filters are called IIR filters (infinite impulse reponse).

If there is no recursive part in the formula above, then it reduces to:

This kind of LTI filters is called FIR filter (finite impulse response). I renamed coefficients b[k] to h[k] to emphasize the fact that now it's just a convolution kernel and the whole operation became the convolution of input signal x[n] and the finite filter kernel h[k] of length N.

Transfer function is another representation of an LTI system, where time delay coefficients now become the polynomial coefficients:

So the denominator is the recursive part, and numerator is combined from non-recursive coefficients. The denominator of FIR filters is simply array of one number {1}.

IIR and FIR filters have their own classes in NWaves - IirFilter and FirFilter respectively. They inherit from LtiFilter class. There's also the third powerful and useful subclass - ZiFilter, and we'll talk about it later.

There are two ways to construct them:

specify TransferFunction object
specify filter coefficients as arrays

The difference is important: in the first case filter object will aggregate TF object; in the second case the internal reference to TF will always be null and calling Tf property will create TransferFunction object on the fly from filter coefficients.

The reason for this design decision is (in short): FirFilter / IirFilter classes are used only for filtering. Everything related to filter design & analysis is concentrated in TransferFunction class. Prior to ver.0.9.2 FDA responsibilities were somewhat mixed between filters and TF (for example, the FrequencyResponse() method was part of filters and not TF, although this concept is related to FDA, not to filtering process per se).

Prior to ver.0.9.2 filter objects contained arrays of both 32-bit and 64-bit coefficients although filtering was (and still is) done with 32-bit floats and FDA was (and still is) done with 64-bit doubles.

As of ver.0.9.2 filters contain only 32-bit coefficients for filtering and the reference to TransferFunction object (not null - only if needed). TF objects contain 64-bit numerator / denominator coefficients.

To sum up, the general simple rule is:

Use LtiFilter subclasses for filtering; use TransferFunction class for FDA.

var b = new[] { 1, 0.2, -0.3 };
var a = new[] { 1, 0.9 };

var filter = new IirFilter(new TransferFunction(b, a));
var tf = filter.Tf;
// here the filter stores the reference to TF with 64-bit precision
// and Tf property returns this reference.
// If we use this filter only for filtering, then this constructor is redundant
// and memory inefficient.

var filter = new IirFilter(b, a);
var tf = filter.Tf;
// in this case TransferFunction will be generated on the fly
// from 32-bit floats (so we'll lose precision).
// Internal Tf object of the filter is still null.
// Memory efficient, but this filter maybe lacks precision for filter analysis.
// If we need just to filter data, this is the best solution.

// Note. Filter coeffs are always 32-bit.
// Although we've passed two arrays of 64-bit doubles they
// will be cast to floats anyway.

More illustrations:

// FIR:
// y[n] = x[n] - 0.2x[n - 1] + 0.6x[n - 2]

var fir = new FirFilter(new [] { 1, -0.2, 0.6 });

var kernel = fir.Kernel;       // { 1, -0.2, 0.6 } (floats)
var kernel = fir.Tf.Numerator; // { 1, -0.2, 0.6 } (doubles but cast from floats)

// IIR:
// y[n] = x[n] - 0.2x[n - 1] + 0.6y[n - 1] - 0.8y[n - 2]

var iir = new IirFilter(new[] { 1, -0.2 }, new [] { 1, -0.6, 0.8 });

var tf = iir.Tf;
var b = tf.Numerator;     // { 1, -0.2 }      (doubles but cast from floats)
var a = tf.Denominator;   // { 1, -0.6, 0.8 } (doubles but cast from floats)


// Same IIR:
// y[n] = x[n] - 0.2x[n - 1] + 0.6y[n - 1] - 0.8y[n - 2]

var trf = new TransferFunction(new[] { 1, -0.2 }, new [] { 1, -0.6, 0.8 });
var fiir = new IirFilter(trf);

// fiir.Tf = trf (they are the same objects)

var bf = fiir.tf.Numerator;     // { 1, -0.2 }      (doubles)
var af = fiir.tf.Denominator;   // { 1, -0.6, 0.8 } (doubles)

Let's take a closer look at the TransferFunction class:

tf.Numerator
tf.Denominator
tf.Gain
tf.FrequencyResponse
tf.ImpulseResponse
tf.GroupDelay()
tf.PhaseDelay()
tf.Zeros
tf.Poles
tf.Normalize
tf.ZpToTf
tf.TfToZp
tf.StateSpace
tf.Zi
tf = tf1 + tf2
tf = tf1 * tf2
tf.ToCsv / tf.FromCsv

Impulse response is the signal produced by LTI filter in response to a unit impulse signal {1, 0, 0, 0, ...}. By default the length of calculated impulse and frequency responses is set to 512.

Frequency response describes behavior of LTI filter in frequency domain. It can be computed as the FFT of impulse response.

Roots of the polynomial in TF numerator are called zeros, and roots of the polynomial in TF denominator are called poles.

The TransferFunction object can be constructed both from numerator/denominator coefficients and from zeros/poles:

// y[n] = x[n] - 0.2x[n - 1] + 0.72x[n - 2] - 0.6y[n - 1]

double[] b = { 1, -0.2, 0.72 };
double[] a = { 1,  0.6 };

var tf1 = new TransferFunction(b, a);


// 2 zeros: 0.2 + 0.3i and 0.2 - 0.3i
// 1 pole:  0.6


// old-style - via ComplexDiscreteSignal:

var zeros = new ComplexDiscreteSignal(1, new[] { 0.2, 0.2 }, new[] { 0.3, -0.3 });
var poles = new ComplexDiscreteSignal(1, new[] { 0.6 },      new[] { 0 });

var tf2 = new TransferFunction(zeros, poles);


// new-style - via Complex[]:

var zeros = new Complex[] { new Complex(0.2, 0.3), new Complex(0.2, -0.3) };
var poles = new Complex[] { (Complex)0.6 };

var tf2 = new TransferFunction(zeros, poles);

Zeros and poles are stored as Complex[].

TransferFunction class provides static methods ZpToTf() and TfToZp() for switching between zero/poles and numerator/denominator coefficients. Note, for root finding the Durand-Kerner algorithm is used, and it becomes numerically unstable if the polynomial order exceeds value about 50. It's an iterative algorithm - you can change max number of iterations if necessary (by default, it's 25000).

var b = TransferFunction.ZpToTf(zeros);
var a = TransferFunction.ZpToTf(poles);

zeros = TransferFunction.TfToZp(b);
poles = TransferFunction.TfToZp(a);

var tf = new TransferFunction(b, a);
tf.CalculateZpIterations = 100000;
var z = tf.Zeros;
var p = tf.Poles;

GroupDelay() and PhaseDelay() methods evaluate the corresponding characteristics for a given transfer function. Both methods accept the size (that can be any positive number).

// y[n] = x[n] + 0.5x[n - 1] + 0.2x[n - 2] + 0.8y[n - 1] - 0.3y[n - 2]

var tf = new TransferFunction(new [] {1, 0.5, 0.2}, new [] {1, -0.8, 0.3});

var zeros = tf.Zeros;
var poles = tf.Poles;

var gd = tf.GroupDelay();
var pd = tf.PhaseDelay(256);

var impulseResponse = tf.ImpulseResponse(256);
var magnitudeResponse = tf.FrequencyResponse().Magnitude;
var phaseResponse = tf.FrequencyResponse().Phase;

Transfer functions can be combined sequentially and parallelly:

var tfp = tf1 + tf2;   // parallel combination
var tfs = tf1 * tf2;   // sequential combination

FIR filters may have pretty long kernels which, furthermore, may be designed in some external tools like MATLAB or sciPy. Hence, TransferFunction class provides methods for reading/writing TF numerators and denominators from/to csv files:

TransferFunction tf;

// read from csv:

using (var stream = new FileStream("tf.csv", FileMode.Open))
{
    tf = TransferFunction.FromCsv(stream);
}

// write to csv:

using (var stream = new FileStream("tf.csv", FileMode.Create))
{
    tf.ToCsv(stream);
}

To sum up:

LtiFilter subclasses:

FirFilter
IirFilter
ZiFilter (conceptually, also an IIR filter)

Their implementations are optimized according to their nature and purposes. So use them accordingly. For example, the following construction is completely valid, and all three filters will produce the same results:

var filter1 = new FirFilter(kernel);
var filter2 = new IirFilter(kernel, new[] { 1.0 });
var filter3 = new ZiFilter(kernel, new[] { 1.0 });

However, obviously, the first version should be preferred, since: 1) it's the FIR filter conceptually; 2) it works faster.

ZiFilter should be preferred over IirFilter in cases when number of numerator coeffs is equal or close to the number of denominator coeffs. Conceptually they don't differ, but the performance of the former is much better in these cases. But what is ZiFilter, anyway?

State space

ZiFilter is the special implementation of IIR filter based on state vector (instead of delay lines).

Any LTI filter with order M has a state space representation A, B, C, D, for which the output y of the filter can be expressed as:

z[n + 1] = A * z[n] + B * x[n]
y[n]     = C * z[n] + D * x[n]

More info about state space representation can be found here.

Transfer function can be constructed from state space:

var ss = new StateSpace { A=a, B=b, C=c, D=d };
var tf = new TransferFunction(ss);

And vice versa - state space can be obtained from TF:

var tf = new TransferFunction(num, den);
var ss = tf.StateSpace;
// then get ss.A, ss.B, ss.C, ss.D

Main features of ZiFilter:

allows setting initial state (initial conditions for filter delays)

var n = Math.Max(num.Length, den.Length);
var zi = new float[n];

// fill zi
// ...

var filter = new ZiFilter(num, den);
filter.Init(zi);
var y = filter.ApplyTo(x);


// sciPy analog:

y = lfilter(num, den, x, zi=zi)

provides method for zero-phase filtering:

var filter = new ZiFilter(num, den);
var y = filter.ZeroPhase(x);

// also you can set padLength:
// var y = filter.ZeroPhase(x, padLength)

// by default padLength is 3*(max(len(num), len(den))-1)
// as in MATLAB, as well


// sciPy analog:

// p = 3*(max(len(num), len(den))-1)
// y = filtfilt(num, den, x, padlen=p)

significantly faster than IirFilter in case when TF numerator and denominator length are equal (e.g. Butterworth, Elliptic, Chebyshev, Bessel filters are implemented as ZiFilters as of ver.0.9.5).

FIR filter design

FIR filters can be designed using following well-known techniques:

Window method (sinc)
Frequency sampling
Remez / Parks-McClellan method (equiripple)

Sinc-window method

// design LP, HP, BP, BS filters:

double[] lpKernel = DesignFilter.FirWinLp(123, 0.12, WindowType.Kaiser);
double[] hpKernel = DesignFilter.FirWinHp(97,  0.32);
double[] bpKernel = DesignFilter.FirWinBp(317, 0.32, 0.42, WindowType.Kaiser);
double[] bsKernel = DesignFilter.FirWinBs(239, 0.32, 0.42);

// for filtering:
var lpFilter = new FirFilter(lpKernel);

// for FDA:
var hpFilter = new FirFilter(new TransferFunction(hpKernel));

Fractional delay design:

var delay = 0.25;

// design LP, HP, BP, BS filters:

double[] lpKernel = DesignFilter.FirWinFdLp(123, 0.12, delay, WindowType.Kaiser);
double[] hpKernel = DesignFilter.FirWinFdHp(97,  0.32, delay);
double[] bpKernel = DesignFilter.FirWinFdBp(317, 0.32, 0.42, delay, WindowType.Kaiser);
double[] bsKernel = DesignFilter.FirWinFdBs(239, 0.32, 0.42, delay);

// + AP (all-pass) filter:

double[] apKernel = DesignFilter.FirWinFdAp(101, delay);

Frequency sampling

Function DesignFilter.Fir() designs FIR filter using frequency sampling approach. It's essentially identical to sciPy's firwin2 function and MATLAB's fir2 function. You need to provide array of ordered frequencies (frequency sampling points) in range [0..0.5] and array of gains (desired magnitude response) at these frequencies). The method linearly interpolates magnitude response from given points to fftSize points, then does IFFT and windowing.

By default, the FFT size is auto-computed. If it is set explicitly, then fftSize/2 + 1 must exceed the filter order.

Array of frequencies can be null. In this case the FFT size must be provided and the array of gains must contain fftSize/2 + 1 values. Frequencies will be uniformly sampled on range [0..0.5].

var order = 19;

// Assuming the following frequencies in sciPy / MATLAB:
// freqs = [0, 0.25, 0.5, 0.54, 0.75, 1.0 ]

// in NWaves frequencies are in range [0..0.5], so divide by 2:
double[] freqs = { 0, 0.125, 0.25, 0.27, 0.375, 0.5 };

// Lowpass filter gain
double[] gains1 = { 1, 1, 1, 0.9, 0, 0 };
// Highpass filter gain
double[] gains2 = { 0, 0, 0.9, 1, 1, 1 };

var kernel1 = DesignFilter.Fir(order, freqs, gains1);
var kernel2 = DesignFilter.Fir(order, freqs, gains2, fftSize: 256);

// for filtering:
var lowpass = new FirFilter(kernel1);

// for FDA:
var highpass = new FirFilter(new TransferFunction(kernel2));

Equiripple FIR filter design

var order = 57;
var freqs = new double[] { 0, 0.15, 0.17, 0.5 };
var response = new double[] { 1, 0 };
var weights = new double[] { 0.01, 0.1 };

var remez = new Remez(order, freqs, response, weights);

var kernel = remez.Design();

// We can monitor the following properties:

// - remez.Iterations
// - remez.ExtremalFrequencies
// - remez.InterpolatedResponse
// - remez.Error


// methods for computing weights from decibel specifications:

var passDb = 3;//dB
var attenuateDb = 40;//dB

var passWeight = Remez.DbToPassbandWeight(passDb);
var attenuateWeight = Remez.DbToStopbandWeight(attenuateDb);


// estimate order of the filter:
int order = Remez.EstimateOrder(freqs, weights);

Or the following methods can be called:

var lpKernel = DesignFilter.FirEquirippleLp(125, 0.12, 0.14, 0.05, 0.95);
var hpKernel = DesignFilter.FirEquirippleHp(211, 0.36, 0.378, 0.05, 0.95);
var bpKernel = DesignFilter.FirEquirippleBp(311, 0.12, 0.14, 0.36, 0.378, 0.05, 0.95, 0.08);
var bsKernel = DesignFilter.FirEquirippleBs(253, 0.12, 0.14, 0.36, 0.378, 0.05, 0.95, 0.08);

LP-HP and BP-BS conversions

var lpKernel = DesignFilter.FirWinLp(345, 0.15);

// HP filter can be obtained from LP with the same cutoff frequency:
var hpKernel = DesignFilter.FirLpToHp(lpKernel);
// and vice versa:
var kernel = DesignFilter.FirHpToLp(hpKernel);

// BP -> BS
var bpKernel = DesignFilter.FirWinBp(123, 0.05, 0.15);

// BS -> BP
var brKernel = DesignFilter.FirBpToBs(bpKernel);
var newKernel = DesignFilter.FirBsToBp(bsKernel);

IIR filter design

The following well-known filters are designed by the same procedure:

Butterworth
Chebyshev-I
Chebyshev-II
Elliptic
Bessel

The procedure is:

compute analog poles (and zeros for some filters) of the prototype LP filter
do bilinear transform to obtain z-domain poles from analog poles
normalize frequency response to have max magnitude = 1.0

Each of these filters is contained within its own namespace. There's a class for every band-form (LP/HP/BP/BS). Example code:

using NWaves.Filters;

vat butter = new Butterworth.BandPassFilter(freq1, freq2, order);
var cheb1 = new ChebyshevI.HighPassFilter(freq, order);
var cheb2 = new ChebyshevII.BandStopFilter(freq1, freq2, order);
var bessel = new Bessel.LowPassFilter(freq, order);

// example how to convert linear scale specifications to decibel scale:

var deltaPass = 0.96;
var deltaStop = 0.04;

var ripplePassDb = Utils.Scale.ToDecibel(1 / deltaPass);
var attenuateDb = Utils.Scale.ToDecibel(1 / deltaStop);

var ellip = new Elliptic.LowPassFilter(freq, order, ripplePassDb, attenuateDb);

Note. For Bessel filters usually the design procedure is different. For better phase behaviour the so called matched z-transform is applied instead of bilinear transform. Hence, resulting filter may be not exactly the one you expect.

Finally, let's consider additional IIR design methods that came from sciPy/MATLAB world:

var freq = 300.0/*Hz*/ / signal.SamplingRate;

var notchTf = DesignFilter.IirNotch(freq, q);
var peakTf = DesignFilter.IirPeak(freq, q);
var combNotchTf = DesignFilter.IirCombNotch(freq, q);
var combPeakTf = DesignFilter.IirCombPeak(freq, q);

var filter = new IirFilter(notchTf);
// or
// var filter = new ZiFilter(notchTf);

Combining filters

Just like transfer functions, LTI filters can be combined sequentially and parallelly. But don't forget that TF objects operate with doubles while filters operate with floats. So, for filtering the following code is OK, but the 64-bit precision will be lost:

Sequential:

var cascade = filter * firFilter * notchFilter;
var filtered = cascade.ApplyTo(signal);

// equivalent to:

var filtered = filter.ApplyTo(signal);
filtered = firFilter.ApplyTo(filtered);
filtered = notchFilter.ApplyTo(filtered);


// same but with double precision ("FDA-first", or "TF-first", approach):

var cascadeTf = filter.Tf * firFilter.Tf * notchFilter.Tf;
var cascadeFilter = new IirFilter(cascadeTf);
var filtered = cascadeFilter.ApplyTo(signal);

Parallel:

var parallel = filter1 + filter2;
filtered = parallel.ApplyTo(signal);


// same but with double precision ("FDA-first", or "TF-first", approach):

var parallelTf = filter1.Tf + filter2.Tf;
var parallelFilter = new IirFilter(parallelTf);
var filtered = parallelFilter.ApplyTo(signal);

Second-order sections

There are also methods for working with Second-Order Sections (SOS) in DesignFilter class:

// get array of SOS from TF:
TransferFunction[] sos = DesignFilter.TfToSos(tf);

// now we can create IIR filters from each TF in sos array...

var filters = sos.Select(tf => new IirFilter(tf));

// then we can apply filters sequentially manually:
foreach (var s in samples)
{
    var y = s;
    foreach (var f in filter)
         y = f.Process(y)
}

// but it's better to use FilterChain class

var filter = new FilterChain(sos);       // constructor from IEnumerable<TransferFunction>
// or
//var filter = new FilterChain(filters); // constructor from IEnumerable<IOnlineFilter>

var y = filter.ApplyTo(x);

// or process samples online:
//    ... outSample = filter.Process(sample);

Reverse operation:

tf = DesignFilter.SosToTf(sos);

// this simply does: sos[0] * sos[1] * ... * sos[n - 1]

Polyphase filters

PolyphaseSystem provides methods for polyphase decimation and interpolation:

// decimate by 4, interpolate by 3:

var firKernel = DesignFilter.FirWinLp(111, 0.5 / Math.Max(3, 4));

var decimator = new PolyphaseSystem(firKernel, 4);
var interpolator = new PolyphaseSystem(firKernel, 3, type: 2);

var decimated = decimator.Decimate(signal);
var interpolated = interpolator.Interpolate(decimated);

Polyphase filter kernels can be retrieved from two properties:

        /// <summary>
        /// Polyphase filters with transfer function E(z^k).
        /// 
        /// Example:
        /// h = [1, 2, 3, 4, 3, 2, 1],  k = 3
        /// 
        /// e0 = [1, 0, 0, 4, 0, 0, 1]
        /// e1 = [0, 2, 0, 0, 3, 0, 0]
        /// e2 = [0, 0, 3, 0, 0, 2, 0]
        /// </summary>
        public FirFilter[] Filters { get; private set; }

        /// <summary>
        /// Polyphase filters with transfer function E(z) used for multi-rate processing.
        /// 
        /// h = [1, 2, 3, 4, 3, 2, 1],  k = 3
        /// 
        /// e0 = [1, 4, 1]
        /// e1 = [2, 3, 0]
        /// e2 = [3, 2, 0]
        /// </summary>
        public FirFilter[] MultirateFilters { get; private set; }

Filtering

LTI filtering can be carried out one of the following ways:

call Process(sample) in a loop
directly (straightforward implementation of the difference equation)
overlap-add block convolution
overlap-save block convolution
custom modes

The second parameter in offline filtering method ApplyTo() is of type enum FilteringMethod that defines the following constants:

Auto
DifferenceEquation
OverlapAdd
OverlapSave
Custom

By default mode is set to Auto, so NWaves will decide what method to use depending on the filter type and kernel size, and most of the time there's no need to tweak this parameter. Rules are simple: with Auto mode, all IIR filters and those FIR filters that have kernel size less than 64, always call Process(sample) in a loop. FIR filters with kernel size starting from 64 samples switch to OverlapSave mode (it'll be faster).

Mode Custom is reserved.

Note that in case of FIR offline filtering the resulting signal has length signal.Length + kernel.Length - 1 according to the definition of convolution operation.

// Process(sample) in loop -  for small kernels, OverlapAdd - for big kernels
// filtered.Length = signal.Length + kernel.Length - 1
var filtered = fir.ApplyTo(signal);

// difference equation
filtered = fir.ApplyFilterDirectly(signal, FilteringMethod.DifferenceEquation);

// overlap save / overlap add
filtered = fir.ApplyTo(signal, FilteringMethod.OverlapSave);
filtered = fir.ApplyTo(signal, FilteringMethod.OverlapAdd);


// Process(sample) in loop; filtered.Length = signal.Length
filtered = iir.ApplyTo(signal);

// overlap save / add (not common option for IIR, since impulse reponse will be truncated)
filtered = iir.ApplyTo(signal, FilteringMethod.OverlapSave);
filtered = iir.ApplyTo(signal, FilteringMethod.OverlapAdd);

You can change the minimum kernel length for automatic switching to OverlapSave mode:

var fir = DesignFilter.FirWinLp(31, 0.16);

fir.KernelSizeForBlockConvolution = 31;
var filtered = fir.ApplyTo(signal);
// OverlapSave

Note, the FFT size in block convolution modes will always be autocomputed as the nearest power of two (kernelSize * 4). If you want to set other parameters then use OlsBlockConvolver or OlaBlockConvolver directly and configure them in constructors.

Online filtering

Online processing is supported by all classes that implement the IOnlineFilter interface. Currently, all filters, block convolvers (OlaBlockConvolver, OlsBlockConvolver) and audio effects contain the Process(sample) and Process(inBuffer, outBuffer) methods responsible for online processing.

Just process data sample after sample:

var outputSample = filter.Process(sample);

Or prepare necessary buffers (or just use them if they come from another part of your system)::

float[] output;

...

void NewChunkAvailable(float[] chunk)
{
    filter.Process(chunk, output);
}

Filters maintain delay lines, so you can easily overwrite the same buffer:

void NewChunkAvailable(float[] chunk)
{
    filter.Process(chunk, chunk);
}

Emulate the frame-by-frame online processing in a loop (only for illustration):

var frameSize = 128;

// big input array (it'll be processed frame-by-frame in a loop)
var input = signal.Samples;

// big resulting array that will be filled more and more after processing each frame
var output = new float[input.Length];

for (int i = 0; i + frameSize < input.Length; i += frameSize)
{
    filter.Process(input, output, frameSize, inputPos: i, outputPos: i);
}

If FIR filter order exceeds approx. 60-65, don't call filter.Process(sample) (filtering will be slow, since in case of online filtering FirFilter does not provide options for using a (faster) block convolution). Switch to block convolvers (OlsBlockConvolver or OlaBlockConvolver) instead:

var kernel = DesignFilter.FirWinLp(431, 0.12);
var filter = new FirFilter(kernel);

var blockConvolver = OlsBlockConvolver.FromFilter(filter, 2048);
// usually the second parameter is approx. 4N, where N is filter order

// or without creating filter object:
var blockConvolver = OlsBlockConvolver(kernel, 2048);

// processing loop:
// while new input sample is available
{
    var outputSample = blockConvolver.Process(sample);
}

blockConvolver.Reset();

// or:
// while new input buffer is available
{
    blockConvolver.Process(input, output);
    // blockConvolver.Process(input, input);  // or in-place
}

LTI filters allow changing their coefficients on-the-fly keeping the same delay line (state). It's very useful if the filter needs to be changed online, during its work:

// while processing
{
    // ...

    // if necessary:
    fir.ChangeKernel(newKernel);

    iir.Change(tf); // transfer function

    // or separately numerator / denominator:
    iir.ChangeNumeratorCoeffs(new[] { 1, 0.2f });
    iir.ChangeDenominatorCoeffs(new[] { 1, 0.5f, -0.3f });

}

New arrays must have the same size as previous arrays of filter coefficients, otherwise no execption will be thrown, but the operation will be ignored.

Also, there's a class called FilterChain that can be used for online and offline filtering if several filters should be applied one after another (i.e. connected sequentially):

var chain = new FilterChain(new[] { filter1, filter2 });

//while (new sample is available)
{
    var outputSample = chain.Process(sample);
    // ...
}

// or chunks:
void GotNewChunk(float[] chunk)
{
    chain.Process(chunk, chunk);
}

chain.Reset();

At any time new filter can be added to the chain or removed:

chain.Add(filter3);
chain.RemoveAt(0);

Numerical problems

In FDA and filtering the numerical problems may arise. Some of them (and how to address them in NWaves) are described here.

IIR filters

IIR filters are normalized by default (i.e. all filter coefficients are divided by Denominator[0]). Still, in some cases maybe it will be necessary to do normalization manually. For this purpose the public method Normalize() is added to TransferFunction class.

One pole filters

These are the simplest IIR filters. Their difference equation is simply:

Transfer function:

The coefficients are calculated automatically based on the given cutoff frequency:

var hpFilter = new OnePole.HighPassFilter(0.2f);

var a = hpFilter.A;
var b = hpFilter.B;

We may change filter coefficients without re-creating the filter itself (in addition to base class methods for changing coefficients directly):

//while (new sample is available)
{
    filter.Process(sample);

    // let's say, if user changed frequency slider position
    filter.Change(newFrequency);
}

One particular example of one pole IIR filter is de-emphasis filter, where a1 coefficient is usually set to some value in range [-0.97, -0.9].

BiQuad filters

Their difference equation is:

Transfer function:

BiQuad filters include:

Lowpass filter
Highpass filter
Bandpass filter
Notch filter
Allpass filter
Peak filter
Low shelf filter
High shelf filter

Coefficients of these filters are computed according to Cookbook formulae for audio EQ biquad filter coefficients by Robert Bristow-Johnson.

Examples:

var q = 0.8;
var frequency = 550.0/*Hz*/;
var notchFilter = new BiQuad.NotchFilter(frequency / signal.SamplingRate, q);

var gain = 6;
var peakFilter = new BiQuad.PeakFilter(frequency / signal.SamplingRate, q, gain);

// by default, gain = 1, q = 1:
var lowShelfFilter = new BiQuad.LowShelfFilter(frequency / signal.SamplingRate);

Parameters q and gain define the shape and height of the frequency response.

Like one-pole filters, BiQuad filters have Change() method.

//while (new sample is available)
{
    peakFilter.Process(sample);

    // let's say, if user changed frequency and gain sliders positions
    peakFilter.Change(newFrequency, gain: newGain);
}

And overloaded version that allows changing all coefficients manually:

filter.Change(b0, b1, b2, a0, a1, a2);

Feedback comb filter

Difference equation:

Transfer function:

var combFilter = new CombFeedbackFilter(25, 1, 0.5);
//var combFilter = new CombFeedbackFilter(m: 25, b0: 1, am: 0.5);


// it also can change its coefficients online:

//while (new sample is available)
{
    combFilter.Process(sample);

    // let's say, if user changed coefficients
    combFilter.Change(newB0, newAm);
}

DC removal filter

Transfer function:

.

r is the parameter of filter (usually in range [0.9, 1]):

var dcrem = new DcRemovalFilter(0.99);

RASTA filter

RASTA filter is the filter developed in speech processing.

Transfer function:

The coefficient a = 0.98 by default (as was proposed in original paper). However, sometimes it can be changed (for instance, in PLP algorithm a = 0.94).

var rasta1 = new RastaFilter();
var rasta2 = new RastaFilter(0.94);

FIR filters

Moving average filter

Moving average filter of length N is very simple:

var maFilter = new MovingAverageFilter(7);
var smoothedSignal = maFilter.ApplyTo(signal);

For longer kernel lengths the formula above can be replaced with more efficient equation, since new sample can be recursively calculated from previous samples:

This equation is implemented in class MovingAverageRecursiveFilter inherited from IirFilter:

var maFilter = new MovingAverageRecursiveFilter(11);
var smoothedSignal = maFilter.ApplyTo(signal);

IIR (recursive) implementation is expectedly much faster. See benchmarks.

Feedforward comb filter

Difference equation:

Transfer function:

var combFilter = new CombFeedforwardFilter(25, 1, -0.5);
//var combFilter = new CombFeedforwardFilter(m: 25, b0: 1, bm: -0.5);


// it can change its coefficients online:

//while (new sample is available)
{
    combFilter.Process(sample);

    // let's say, if user changed coefficients
    combFilter.Change(newB0, newBm);
}

Pre-emphasis filter

It's a very simple filter widely used in speech processing:

var filter = new PreEmphasisFilter(0.97);

// note the sign is different! the kernel will be: {1, -0.97}

Usually, pre-emphasis filters are not normalized (i.e. max value of magnitude repsonse is not 1.0). Set the optional parameter normalize in constructor to true for normalization.

Savitzky-Golay filter:

Supported window sizes are [5, 7, 9, .., 31] (derivative kinds: 0, 1, 2):

var savgol1 = new SavitzkyGolayFilter(9);
var savgol2 = new SavitzkyGolayFilter(17, 2);

Adaptive filters

All adaptive filters are inherited from the abstract base class AdaptiveFilter which, in turn, is derived from FirFilter class. Subclasses implement Process(sample, desired) method. Adaptive filters adapt their weights (kernels) based on a particular algorithm:

LmsFilter
NlmsFilter
SignLmsFilter
VariableStepLmsFilter
LmfFilter
NlmfFilter
RlsFilter

Example:

float mu = 0.95f;
AdaptiveFilter filter = new NlmsFilter(5, mu);

// initialize weights if necessary (zeros by default)
// filter.Init(weights);

// suppose, we have input and desired signals

var adapted = Enumerable.Range(0, input.Length)
                        .Select(i => filter.Process(input[i], desired[i]));

var adaptedWeights = filter.Kernel;

Miscellaneous

Median filter

There are 2 implementations of median filtering: MedianFilter and MedianFilter2. Both are identical to scipy.signal.medfilt(). The former class should be preferred in most cases. The latter works slightly faster for small sizes of filter (size <= 5, approx.). See benchmarks.

Wiener filter

Actually, not an adaptive Wiener filter. Implementation is identical to scipy.signal.wiener().

64-bit filters

The following classes are available for working with double precision:

IFilter64
IOnlineFilter64
FirFilter64
IirFilter64
OlsBlockConvolver64
OlaBlockConvolver64
StereoFilter64
FilterChain64

Algorithmically they are completely identical to their 32-bit counterparts. Method ApplyTo() accepts double[] instead of DiscreteSignal.

var butter = new Butterworth.LowPassFilter(0.12, 4);
var filter = new IirFilter64(butter.Tf);

double[] samples/* = GetSamples()*/;

double[] filtered = filter.ApplyTo(samples);

//while (new sample is available)
{
    double outputSample = filter.Process(sample);
    // ...
}

filter.Reset();

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Filters

General notes

Filter design and analysis (FDA)

Transfer function and difference equation

State space

FIR filter design

Sinc-window method

Frequency sampling

Equiripple FIR filter design

LP-HP and BP-BS conversions

IIR filter design

Combining filters

Second-order sections

Polyphase filters

Filtering

Online filtering

Numerical problems

IIR filters

One pole filters

BiQuad filters

Feedback comb filter

DC removal filter

RASTA filter

FIR filters

Moving average filter

Feedforward comb filter

Pre-emphasis filter

Savitzky-Golay filter:

Adaptive filters

Miscellaneous

Median filter

Wiener filter

64-bit filters

Clone this wiki locally