Filters

Contents:

General notes
Filter design and analysis
Filtering
IIR filters
FIR filters
Adaptive filters
Miscellaneous
1. Median filter
2. Wiener filter

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 the Process() method sample after sample. And this is how most filters are implemented in NWaves.

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 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 of this design decision is (in short): FirFilter / IirFilter classes are used only for filtering. Everything related to filter design & analysis as of ver.0.9.2 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 = tf1 + tf2
tf = tf1 * tf2
tf.ToCsv / tf.FromCsv
tf.ZpToTf
tf.TfToZp

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

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

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. 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 of FFT as parameter (512, by default).

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

FIR filter design

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

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

Due to NWaves limitations, in all cases filter kernel size must be odd number, otherwise exception will be thrown.

Frequency sampling

You need to specify uniformly sampled frequency response (magnitude and optionally phase):

double[] kernel1 = DesignFilter.Fir(43, magnitudeResponse);
double[] kernel2 = DesignFilter.Fir(43, magnitudeResponse, phaseResponse, windowTypes.Kaiser);

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

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

Sinc-window method

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

double[] lpKernel = DesignFilter.FirWinLp(123, 0.12, windowTypes.Kaiser);
double[] hpKernel = DesignFilter.FirWinHp(97,  0.32);
double[] bpKernel = DesignFilter.FirWinBp(317, 0.32, 0.42, windowTypes.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));

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.

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

Parallel:

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

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

TransferFunction[] sos = DesignFilter.TfToSos(tf);

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

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:

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.

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:

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

Some filters allow changing their coefficients on-the-fly keeping the same delay line (state). In next version I'm planning to add this capability to all LTI filters.

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

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

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

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:

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:

This equation is implemented in class MovingAverageRecursiveFilter inherited from IirFilter:

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

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}

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

LmsFilter
NlmsFilter
SignLmsFilter
VariableStepLmsFilter
LmfFilter
NlmfFilter
RlsFilter

Example:

AdaptiveFilter filter = new NlmsFilter(5, mu);

// 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.Weights;

Miscellaneous

Median filter

Identical to scipy.signal.medfilt().

Wiener filter

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

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

FIR filter design

Frequency sampling

Sinc-window method

Equiripple FIR filter design

LP-HP and BP-BS conversions

IIR filter design

Combining filters

Polyphase filters

Filtering

Online filtering

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

Clone this wiki locally