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):
But applying it backward (precompensated_signal
) to precompensate yields an oscillating signal:
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
):
and behaves a little better in the backward direction (precompensated_signal_2
):
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
):
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
- Explaining what the modified filter corresponds to
- Finding the filter from the answer from scipy
- 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).
b = [0.0305, 0.0289]
anda = [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$