From your 2nd picture it looks like you need a linear-phase FIR, as opposed to an IIR. The main difference is that FIRs can have perfectly linear phase, thus a perfectly flat group delay, which means all the frequencies are delayed by the same amount, while IIRs have non-linear phase by design, thus a nonlinear group delay, resulting in smearing in time domain.
For a [0.2, 1.6] Baud passband you could use a 0.4 Baud transition width (\$\omega_{tw}\$), so that the first band will be [0, 0.4] Baud (thus \$\omega_c^1=0.2\$), and the second [1.4, 1.8] Baud (\$\omega_c^2=1.4\$).
You linked to Python's least squares FIR, and that's a good start. If the coefficients need to be calculated in real-time, the algorithm might pose problems, because it involves a bit of math, but the closest alternative is the Kaiser window (which is a closed form converging towards a LS). I am not very familiar with Python, but I'll show a possible filter that you could use.
The order of the filter is directly related to the transition width, and for an LS there is no closed form, but the LS and the Kaiser window come very close in terms of response, and Kaiser has a known formula for determining the order. Considering an attenuation of 60 dB and a sampling of 8 Baud:
$$\lceil\frac{A_s-7.95}{2.285\omega_{tw}}\rceil=\lceil\frac{A_s-7.95}{2.285\pi\frac{0.4\,\text{Baud}}{4\,\text{Baud}}}\rceil=73$$
If the sampling frequency would have been 4, the order would have been halved; the opposite for higher sampling rates. This is how a 37th order Kaiser would look like in frequency domain:
And some random sequence, vaguely similar to a bi-phase differential encoding:
The output seems to wobble a bit because of the transient response. If the input were steady for more periods, it would stabilize. Anyway you look at it, it gives a correct band-limited signal, that looks much like what you have in your 2nd picture, except less smooth. Note that I have no idea what signal is used there, if there is oversampling or not, or whether that's the analog output for the sampled version of the signal. This is meant to show that it's a linear phase filter that's needed in order to correctly preserve the shape of the signal. Otherwise, considering a 4x oversampling, the difference would look like what you see in the lower graph, in the picture above.
Finally, as I'm not very familiar with Python, in Octave, a Kaiser window FIR can be designed like this:
w = [0 0.4 1.4 1.8];
[N, w, b, ftype] = kaiserord(w, [0 1 0], [1 1 1]*0.001, 8)
h = fir1(N-1, w, kaiser(N, b));
freqz(h) # check the response
You could also concoct a signal and use filt(), or conv() to show the response. I already did it above, but here's what Octave shows for a similar pseudo coded input:
# a quick'n'dirty setup
r = randi(2, 1, 20) - 1;
r = [r; r](:).';
s = [ones(1, 20); zeros(1,20)](:).';
p = xor(r, s);
q = [p; p; p; p](:).';
stem([zeros(1, ceil(N/2)), q]);
hold on;
stem(conv(h, q));
hold off;
In a similar fashion you can use firls(), or remez() if you want equiripple, or whatever other filter, but, unless you have a different plan, you'll most probably need linear phase. You've seen what happens with a Butterworth filter, which is an IIR. As a last example, here's how the same input is filtered by a 5th order inverse Chebyshev (to try to match the frequency response of the FIR), and how the 4x interpolated shows up:
