"Simplest Lowpass Filter" in Python

I'm currently working my way through Introduction to Digital Filters with Audio Applications by Julius O. Smith III. The first section deals with what the book calls the simplest lowpass filter. After the initial mathematical introduction (both via tedious trigonometry and a complex approach) there's an actual implementation. The book uses Matlab, but I much prefer Python for this type of job even though I should probably use GNU Octave for this job but I really prefer Python. These are the notes I took when converting that initial script.

The filter is described by:

$$ y(n) = x(n)+x(n-1), n = 1,2,...,N $$

This is the implementation in Matlab:

 0% simplpm1.m - matlab main program implementing
 1%              the simplest lowpass filter:
 2%
 3%                    y(n) = x(n)+x(n-1)}
 4
 5N=10;       % length of test input signal
 6x = 1:N;    % test input signal (integer ramp)
 7B = [1,1];  % transfer function numerator
 8A = 1;      % transfer function denominator
 9
10y = filter(B,A,x);
11
12for i=1:N
13  disp(sprintf('x(%d)=%f\ty(%d)=%f',i,x(i),i,y(i)));
14end

There's not much to it, we need to figure out what the scipy alternative is to the Matlab filter function and make sure the indices are OK.

This is the translation I came up with:

 0from scipy.signal import lfilter
 1import numpy as np
 2
 3N=10
 4x = np.arange(1.0,N+1)
 5B = [1,1]
 6A = 1
 7
 8y = lfilter(B,A,x);
 9
10for i in range(N):
11    print(f"x({i+1})={x[i]} \t y({i+1})={y[i]}")

As far as I can tell the Matlab function filter correspons the most to the Python function lfilter. There's also the need to convert the code from 1-indexed to 0-indexed.

This is the same result. scipy.signal.lfilter takes three parameters: b, a and x. The documentation mentions that scipy.signal.sosfilt is actually the preferred function. That's something to figure out at a later point though.

A second part of the exercise is a block oriented approach of the filter, where the filter state is saved and added as an additional parameter on the next call of the function.

 0% simplpm2.m - block-oriented version of simplpm1.m
 1
 2N=10;       % length of test input signal
 3NB=N/2;     % block length
 4x = 1:N;    % test input signal
 5B = [1,1];  % feedforward coefficients
 6A = 1;      % feedback coefficients (no-feedback case)
 7
 8[y1, Sf] = filter(B,A,x(1:NB));      % process block 1
 9      y2 = filter(B,A,x(NB+1:N),Sf); % process block 2
10
11for i=1:NB   % print input and output for block 1
12  disp(sprintf('x(%d)=%f\ty(%d)=%f',i,x(i),i,y1(i)));
13end
14
15for i=NB+1:N % print input and output for block 2
16  disp(sprintf('x(%d)=%f\ty(%d)=%f',i,x(i),i,y2(i-NB)));
17end

The python implementation is a bit more finicky though:

The main difference is that you have to explicitly set the initial condition for the filter for the scipy.signal.lfilter function to actually return the state zf necessary for the next block. Setting it to [0.0] basically means that the filter is at initial rest. It would probably be slightly more correct to use scipy.signal.lfiltic for this, but it wouldn't make any difference. After that it's just a matter of passing zf into the computation of the second block.

I found it easier to concatenate both blocks instead of having two loops. Because these are numpy.ndarray you can't do the typical python trick of using the + operator. Using numpy.concatenate does the trick though.

 0import numpy as np
 1
 2N=10
 3NB=int(N/2)
 4x = np.arange(1.0,N+1)
 5
 6B = [1,1]
 7A = 1
 8
 9y1, zf = lfilter(B,A,x[0:NB],zi=[0.0])
10y2, _ = lfilter(B,A,x[NB:N],zi=zf)
11
12y = np.concatenate((y1,y2))
13
14for i in range(N):
15    print(f"x({i+1})={x[i]} \t y({i+1})={y[i]}")

The script gives the following output:

This is exactly the same as above. Scipy and Numpy really are amazing tools and in my experience really good at replacing something like Matlab, but there are some pitfalls you have to watch out for when trying to do the conversion.