2
$\begingroup$

Question

Is there an invertible low-pass filter built into scipy.signal (or other python package)? If so, what is it? If not, why not (is there something particularly difficult about inverting a low-pass filter)?

Elaboration

I need an invertible (digital, first-order, for concreteness) low-pass filter, such as a butterworth filter. I need to apply the inverse of it to precompensate an analog signal with $(x,y) = (\mathrm{time}, \mathrm{voltage})$. The signal is converted from digital to analog after filtering.

Work done so far

Naively, (as far as I understand) I can do the following:

import numpy as np 
import scipy
# define signal here
num_order = 1
cutoff_frequency = 1e4 # Hz 
time_constant = 1/(2*np.pi*cutoff_frequency) 
lpf_b, lpf_a = scipy.signal.butter(num_order, btype='low', Wn=1/time_constant/(2*np.pi*sample_rate/2))
precompensated_signal = scipy.signal.lfilter(lpf_a, lpf_b, signal)

where an example signal is defined, as (e.g.)

sample_rate = 1e6
t_period = 2*1e-3
num_samples = int(round(t_period*sample_rate))
amplitude = 0.23
signal = amplitude*np.ones(num_samples)

This approach works when inverting a high-pass filter.

Applying the filter in the forward-direction (filtered_signal)

filtered_signal = scipy.signal.lfilter(lpf_b, lpf_a, signal)

produces the expected signal (check by running the code or see image): scipy.signal.butter in forward direction

But applying it backward (precompensated_signal) to precompensate yields an oscillating signal:

scipy.signal.butter in backward direction

Looking at the coefficients, I find

lpf_a = [ 1.0, -0.93906251]
lpf_b = [0.03046875, 0.03046875]

Seemingly having two identical coefficients in b as the first and only coefficients is what renders the filter non-invertible.

There is a related question and answer from 2016. Based on the answer, taking the filter coefficients given by scipy.signal.butter and modifying them as follows

lpf_b, lpf_a = scipy.signal.butter(num_order, btype='low', Wn=1/time_constant/(2*np.pi*sample_rate/2))
lpf_b_2 = [1 + lpf_a[1]] # note that lpf_b_2[0] = lpf_b[0] + lpf_b[1] 
filtered_signal_2 = scipy.signal.lfilter(lpf_b_2, lpf_a, signal)
precompensated_signal_2 = scipy.signal.lfilter(lpf_a, lpf_b_2, signal)

yields a filter which behaves identically to the original one in the forward direction (filtered_signal_2):

precompensated signal with modified filter (only one b coefficient)

and behaves a little better in the backward direction (precompensated_signal_2):

enter image description here

though the resulting signal is still questionable. The only coefficient that differs from amplitude is the first one (and this seems to be true even when changing the cutoff frequency), and worryingly then the voltage at the first datapoint gets very large. This feature seems to be true for the filter in the above linked answer.

Edit

ZR Han suggested (thanks!) to shift the pole of the filter. As far as I understand, it can be done as follows:

z, p, k = scipy.signal.tf2zpk(lpf_b, lpf_a)
# the zero z = [-1.0] as he suggested
z = [-0.95]
lpf_b, lpf_a = scipy.signal.zpk2tf(z, p, k)

This results in

lpf_b = [0.03046875, 0.02894531]
lpf_a = [ 1.0, -0.93906251]

But unfortunately this did not remove all of the oscillation, as can be seen by plotting the precompensated signal (scipy.signal.lfilter(lpf_a, lpf_b, signal) with the modified coeffs):

precompensated signal with zero moved towards the unit circle

Perhaps this is caused by a poor choice in z/my lack of formal understanding of what is "a small shift". Since the precompensated signal is converted to analog, ill behaviour is not tolerated in the solution unfortunately.

Directions to go from here

  1. Explaining what the modified filter corresponds to
  2. Finding the filter from the answer from scipy
  3. Finding another filter which is better suited to be used with signal that is converted to analog (smooth and finite)

Context: I have no formal background in digital signal processing, but I can learn what I need (faster if I'm pointed to a source).

$\endgroup$
7
  • $\begingroup$ you say you need a digital filter but it will be applied to an analog signal – I'm confused. Are you planning to build an analog filter from a digital prototype to apply it to an analog signal, or planning to digitize an analog signal to apply a digital filter to it? You can't directly apply a digital filter to an analog signal. $\endgroup$ Jul 4 at 19:40
  • $\begingroup$ Regarding the directions to go from there – it might really help to tell us what you need this for! We might be able to help with that :) $\endgroup$ Jul 4 at 19:45
  • 1
    $\begingroup$ Do you mean $H_1(z) \cdot H_2(z) = 1$ by inverse of a digital filter? Exchanging the denominator and the numerator gives you a stable inversion as long as the zeros of the original transfer function are all inside the unit circle. $\endgroup$
    – ZR Han
    Jul 5 at 1:16
  • 1
    $\begingroup$ Your filter has a zero at $z=-1$ and a pole at $z=0.939$ so it can't be invertible. You can move the zeros slightly inside the unit circle such as $z=-0.95$ and then the filter coefficients would be b = [0.0305, 0.0289] and a = [1.0000, -0.9391], which is also a low-pass filter but its attenuation at $\pi$ is a little bit lower compared with the original filter. $\endgroup$
    – ZR Han
    Jul 5 at 1:36
  • 1
    $\begingroup$ As Hilmar suggests, a fundamental knowledge of DSP would help you to understand the problem. The filter you designed has a zero at $z=-1$ to obtain a sharp stopband at Nyquist frequency. A digital filter is stable iff all its poles are inside the unit circle. As you swap the numerator and denominator, which means that poles and zeros exchange their position, results in a pole on the unit circle and consequently an unstable filter. $\endgroup$
    – ZR Han
    Jul 6 at 4:56
4
$\begingroup$

Is there an invertible low-pass filter

No

is there something particularly difficult about inverting a low-pass filter?

Yes. Digital low pass filters (in the most common sense) have a zero at Nyquist which means that the inverse has infinite gain at Nyquist and is unstable.

Seemingly having two identical coefficients in b as the first and only coefficients is what renders the filter non-invertible.

Yes. That creates the zero at Nyquist.

a filter which behaves identically to the original one in the forward direction

It does not. Changing the "b" coefficient to a single one will move the pole away from the unit circle on the real axis which changes the high frequency behavior. The signal only looks the same since it has very little high frequency content

Depending on you application you are probably better off with complimentary shelf filters. A low-pass filter destroys information which is unrecoverable. In practice if any filter has a lot of attenuation the inverse will have a lot of gain which is problematic in terms of noise performance.

Context: I have no formal background in digital signal processing, but I can learn what I need (faster if I'm pointed to a source).

Digital Signal Processing is fairly math heavy and it really helps to have some idea about the theoretical foundations. A good starter could be:

https://ccrma.stanford.edu/~jos/filters/

$\endgroup$
3
  • $\begingroup$ Thank you for this, this information is very helpful. I will see if there is an algorithmic python solution that works well enough, before letting go of the question - the practical implementation is what I'm looking for (with sufficient connection to theory that I can explain it formally if necessary) - hence the part "what is it" in the question. The solution does not need to be exact/ideal, it only needs to add gain at higher frequencies without distorting the signal up to 1-2 O magnitude further. $\endgroup$ Jul 5 at 11:41
  • $\begingroup$ Your point about noise is worth considering, though for now I will assume that the circuit with the low pass filter cannot be modified. Requirement of higher output amplitude at high frequencies will mean that the SNR at higher frequencies is larger, assuming white noise coming out of the instrument? Or do you mean something else? $\endgroup$ Jul 5 at 11:47
  • 1
    $\begingroup$ Shelving filters here: gist.github.com/endolith/5455375#file-biquad_cookbook-py-L455 $\endgroup$
    – endolith
    Jul 5 at 12:56
0
$\begingroup$

Following up the suggestion to use so-called shelving filters by Hilmar (thanks!), with the biquad_cookbook module that endolith linked (thanks!). The package (written in 2013) runs with (conda) python 3.8.10 out of the box.

Mini example:

import numpy as np
import scipy
import biquad_cookbook # need the file in the working directory

precompensated_signals = []
precompensated_and_filtered_signals = []
for g in [3, 6, 9, 12]:
    shelf_b, shelf_a = biquad_cookbook.shelf(Wn=1/time_constant/(2*np.pi*sample_rate/2), 
        dBgain=g, S=1, btype='high', ftype='half', analog=False, output='ba')
    precompensated_signals.append(scipy.signal.lfilter(shelf_b, shelf_a, signal))
    precompensated_and_filtered_signals.append(scipy.signal.lfilter(lpf_b, lpf_a, 
        scipy.signal.lfilter(shelf_b, shelf_a, signal)))

where lpf_a and lpf_b are defined as in the question, to yield the smooth, finite, albeit slightly oscillating precompensated signals (precompensated_signals; note that we applied the shelf filter in the forward direction - the filter is directly the precompensating filter):

enter image description here

Taking these through the original low-pass filter (precompensated_and_filtered_signals modelled as the first order butterworth filter from the question) yields the following signals:

enter image description here

So, this is a partial solution. For example, at 12 dB gain, the signal reaches the desired amplitude in roughly tenth of the time compared with the non-precompensated signal. The tradeoff is that the amplitude will overshoot due to the ripple in the precompensated signal (created by the filter).

$\endgroup$
2
  • $\begingroup$ One could still ask whether the ripple is necessary in this kind of filter? One could naively expect that an (carefully-defined) exponential decay corresponds to a very similar filter, but without the ripple. Is the problem that it would require more coefficients to define, or something else? $\endgroup$ Jul 6 at 16:12
  • $\begingroup$ I found my answer to my above question: exponentially decaying waveform will also overshoot when going through a low-pass filter. $\endgroup$ Jul 6 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.