Delay filters in Python using for loop and $\tt lfilter$

Question

I have these two filters as given in recursive form: \begin{align} y_1[i] &= y[i] - y[i-d] + 1.2y_1[i-d] - 0.8y_1[i-2d]\\ y_2[i] &= y[i] - 2y[i-d] + y[i-2d] \end{align}

where $d$ is an integer for setting the delay. I wanted to have a look and see how such filters change input signals depending on the value of $d$, but I faced the problem with implementing the actual filters in Python: I tried iterations and scipy.signal.lfilter, and each method was giving me different results. Let me show you the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lfilter
# The test signal is just a sine of the length l, constant amplitude of 1
# and varying frequency w
l = 1024
x = np.array(range(l))
w = x/l
y = np.sin(w*x)
d = 4

# Let's caclulate the signals filtered with IIR filter
## First using forward iteration
y11 = y.copy()
for i in range(2*d, l):
    y11[i] = y[i] - y[i-d] + 1.2*y11[i-d] - 0.8*y11[i-2*d]
## Second using backward iteration
y12 = y.copy()
for i in reversed(range(2*d, l)):
    y12[i] = y[i] - y[i-d] + 1.2*y12[i-d] - 0.8*y12[i-2*d]
## Third using lfilter
y13 = lfilter([1, 0, 0, 0, 1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)

# Let's caclulate the signals filtered with FIR filter
## First using forward iteration
y21 = y.copy()
for i in range(2*d, l):
    y21[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Second using backward iteration
y22 = y.copy()
for i in reversed(range(2*d, l)):
    y22[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Third using lfilter
y23 = lfilter([1, 0, 0, 0, 2, 0, 0, 0, -1], [1], y)

# Let's start the figure and plot the signal in the top row, then the filtered
# signals below it
fig = plt.figure()
ax1 = plt.subplot(421)
plt.title('Test signal')
plt.plot(x, y)
ax2 = plt.subplot(422)
plt.title('Test signal')
plt.plot(x, y)

# The left column is for the IIR-filtered signals
plt.subplot(423, sharex=ax1)
plt.title('Forward iter.')
plt.plot(x, y11)
plt.subplot(425, sharex=ax1)
plt.title('Backward iter.')
plt.plot(x, y12)
plt.subplot(427, sharex=ax1)
plt.title('lfilter')
plt.plot(x, y13)

# The right column is for the FIR-filtered signals
plt.subplot(424, sharex=ax2)
plt.title('Forward iter.')
plt.plot(x, y21)
plt.subplot(426, sharex=ax2)
plt.title('Backward iter.')
plt.plot(x, y22)
plt.subplot(428, sharex=ax2)
plt.title('lfilter')
plt.plot(x, y23)

plt.tight_layout()
plt.show()

And here's the resulting figure:

Judging by the look of the plots in rows 2 to 4 one can see that each calculation produced different results except the 2nd and 3rd plot in the right column where backward and forward iteration worked out similarly. Normally I would rely on the results from Scipy's lfilter, but there I'm not sure if I specified the coefficients correctly. Somebody please take look and point out the mistakes.

Answer 1

After a trial-and-error session I figured what was wrong in my code:

• The coefficients I fed to lfilter had incorrect signs, because only the signs of the coeffs from the a array have to be reversed and those in the array b stays the same.

• The backward iteration did not work properly for the IIR filter, because when running backwards the delayed values were not yet filtered and in that case equal to the input signal.

Here is the fixed code:

l = 1024
x = np.array(range(l))
w = x/l
y = np.sin(w*x)
d = 4

# Let's caclulate the signals filtered with IIR filter
## First using forward iteration
y11 = y.copy()
for i in range(2*d, l):
    y11[i] = y[i] - y[i-d] + 1.2*y11[i-d] - 0.8*y11[i-2*d]
## Second using backward iteration
y12 = y.copy()
for i in reversed(range(2*d, l)):
    y12[i] = y[i] - y[i-d] + 1.2*y12[i-d] - 0.8*y12[i-2*d]
## Third using lfilter
# was:
#y13 = lfilter([1, 0, 0, 0, 1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)
y13 = lfilter([1, 0, 0, 0, -1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)

# Let's caclulate the signals filtered with FIR filter
## First using forward iteration
y21 = y.copy()
for i in range(2*d, l):
    y21[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Second using backward iteration
y22 = y.copy()
for i in reversed(range(2*d, l)):
    y22[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Third using lfilter
# was:
#y23 = lfilter([1, 0, 0, 0, 2, 0, 0, 0, -1], [1], y)
y23 = lfilter([1, 0, 0, 0, -2, 0, 0, 0, 1], [1], y)

# Let's start the figure and plot the signal in the top row, then the filtered
# signals below it
fig = plt.figure()
ax1 = plt.subplot(421)
plt.title('Test signal')
plt.plot(x, y)
ax2 = plt.subplot(422)
plt.title('Test signal')
plt.plot(x, y)

# The left column is for the IIR-filtered signals
plt.subplot(423, sharex=ax1)
plt.title('Forward iter.')
plt.plot(x, y11)
plt.subplot(425, sharex=ax1)
plt.title('Backward iter.')
plt.plot(x, y12)
plt.subplot(427, sharex=ax1)
plt.title('lfilter')
plt.plot(x, y13)

# The right column is for the FIR-filtered signals
plt.subplot(424, sharex=ax2)
plt.title('Forward iter.')
plt.plot(x, y21)
plt.subplot(426, sharex=ax2)
plt.title('Backward iter.')
plt.plot(x, y22)
plt.subplot(428, sharex=ax2)
plt.title('lfilter')
plt.plot(x, y23)

plt.tight_layout()
plt.show()

Now on the plots it is clearly visible that iterations produce the same results as scipy's lfilter except for the backward iteration attempt at the IIR filter as it was explained above.