Skip to content

Filters

Jump to bottom
Tim edited this page Feb 25, 2019 · 19 revisions

Contents:

  1. General notes
  2. Filter design and analysis
    1. Transfer function and difference equation
    2. Design lowpass, highpass, bandpass, bandreject filters
    3. Combining filters
  3. Filtering
  4. IIR filters
    1. One pole filters
    2. BiQuad filters
    3. Feedback comb filter
    4. Butterworth filter
    5. DC removal filter
    6. RASTA filter
  5. FIR filters
    1. Moving average filter
    2. Feedforward comb filter
    3. Pre-emphasis filter
  6. Miscellaneous
    1. Median filter
    2. Wiener filter

General notes

In the most general sense, Filter, or 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 the Process() method sample after sample. And this is how most filters are implemented in NWaves. However, there are also certain filters (including some audioeffects) that offer only offline functionality so far.

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.

Main concepts related to LTI filters are:

  • impulse response (ltiFilter.ImpulseResponse(256))
  • frequency response (ltiFilter.FrequencyResponse(256))
  • transfer function (ltiFilter.Tf)

Impulse response is the signal produced by LTI filter in response to a unit impulse signal {1, 0, 0, 0, ...}.

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

By default the length of calculated impulse and frequency responses is set to 512.

Filter design and analysis

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:

img

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:

img

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:

img

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, which, in turn, implements interfaces IFilter and IOnlineFilter. In order to construct them, the corresponding difference equation coefficients must be specified:

// 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 b = iir.B;             // { 1, -0.2 }
var a = iir.A;             // { 1, -0.6, 0.8 }


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

TransferFunction class:

  • tf.Numerator
  • tf.Denominator
  • tf.GroupDelay()
  • tf.PhaseDelay()
  • tf = tf1 + tf2
  • tf = tf1 * tf2

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

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

Note, that in filter design and analysis we operate with double precision. Zeros and poles are given as objects of convenient class ComplexDiscreteSignal.

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.

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

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

GroupDelay() and PhaseDelay() methods evaluate the corresponding characteristics for a given transfer function. Both methods accept the size of FFT as parameter.

All LTI filters have the property Tf. Underlying object is created and filled automatically during filter construction depending on a type of the filter.

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

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

var zeros = filter.Tf.Zeros;
var poles = filter.Tf.Poles;

var gd = filter.Tf.GroupDelay();
var pd = filter.Tf.PhaseDelay();

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

Transfer functions can be combined sequentially and parallelly:

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

Design lowpass filter

There are several ways to design lowpass filter using window method:

var lpFilter1 = DesignFilter.Fir(43, magnitudeResponse);
var kernel1 = lpFilter1.Tf.Numerator;

var lpFilter2 = DesignFilter.Fir(43, magnitudeResponse, phaseResponse, windowTypes.Kaiser);
var kernel2 = lpFilter2.Tf.Numerator;

var lpFilter3 = DesignFilter.FirLp(123, 0.12f, window: windowTypes.Kaiser);
var kernel3 = lpFilter3.Tf.Numerator;

Fir method accepts filter order, array of magnitude response samples and optionally: array of phase response samples and the window type. FirLp method accepts filter order, cutoff frequency (as a fraction of sampling rate) and optionally the window type. In both methods by default the Blackman window is used.

Design highpass, bandpass, bandreject filters

var lpFilter = DesignFilter.FirLp(345, 0.15f);

// HP filter can be obtained from LP with the same cutoff frequency:
var hpFilter = DesignFilter.LpToHp(lpFilter);
// and vice versa:
var lowpass  = DesignFilter.HpToLp(hpFilter);

// design BP filter
var bpFilter = DesignFilter.FirBp(123, 0.05f, 0.15f);

// design BR filter
var brFilter = DesignFilter.FirBr(201, 0.08f, 0.23f, WindowTypes.Kaiser);

In all cases filter kernel size must be odd number, otherwise exception will be thrown.

Combining filters

Just like transfer functions, LTI filters can be combined sequentially and parallelly (they just delegate this work to their Tf objects).

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

Parallel:

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

Filtering

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

  • directly (implement difference equation in a loop)
  • 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: IIR filters always use DifferenceEquation mode. FIR filters with kernel size less than 64 use DifferenceEquation mode, otherwise OverlapAdd.

Mode Custom is reserved.

// difference equation for small kernels, OverlapAdd for big kernels
var filtered = fir.ApplyFilterDirectly(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);


// difference equation always
filtered = iir.ApplyFilterDirectly(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);

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.

Simply 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);
}

Sample after sample:

var outputSample = filter.Process(sample);

Emulate the frame-by-frame online processing in a loop:

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);
}

Also, there's a class called FilterChain that can be used for online filtering if several filters should be applied one after another:

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

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

chain.Reset();

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

chain.Add(filter3);

IIR filters

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

One pole filters

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

img

Transfer function:

img

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 not only create the filter but also call MakeTf() method to see what the coefficients will be (without creating the filter itself):

using NWaves.Filters.OnePole;

var b = new double[2];
var a = new double[2];

LowPassFilter.MakeTf(0.12f, b, a);

// arrays b and a now contain the corresponding coefficients of the filter

Also, we can change filter coefficients without re-creating the filter itself (thus preserving its state). It's very useful if the filter needs to be changed online, during its work:

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

img

Transfer function:

img

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() and MakeTf() methods.

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

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

Feedback comb filter

Difference equation:

img

Transfer function:

img

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);
}

Butterworth filter

Butterworth LP filter coefficients are calculated automatically from order and cutoff frequency:

var butterworth = new ButterworthFilter(0.2, 7);
//var butterworth = new ButterworthFilter(freq: 0.2, order: 7);

DC removal filter

Transfer function:

img.

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:

img

var rasta = new RastaFilter();

FIR filters

FIR filters have a couple of specific methods.

Firstly, FIR filters may have pretty long kernels which, furthermore, may be designed in some external tools. Hence, FirFilter class provides methods for reading/writing kernels from/to csv files:

FirFilter filter;

// read from csv:

using (var stream = new FileStream("kernel.csv", FileMode.Open))
{
    filter = FirFilter.FromCsv(stream);
}

// write to csv:

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

Secondly, FIR filters can be converted to IIR (with TF denominator {1}):

var iir = fir.AsIir();

Moving average filter

Moving average filter of length N is very simple:

img

Size N must be odd number.

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:

img

This equation is implemented in class MovingAverageRecursiveFilter:

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

Feedforward comb filter

Difference equation:

img

Transfer function:

img

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:

img

var filter = new PreEmphasisFilter(0.97);

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

Miscellaneous

Median filter

Identical to scipy.signal.medfilt().

Wiener filter

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

Clone this wiki locally