Combining Bandpass filters on 1D signal to obtain better SNR?
If you know that the interference is only in a given frequency range, you might use a notch filter to filter out this interference. If it is in several frequency ranges, you can either filter with different filters in parallel and use the output that works best for you, or filter the signal serially through all filters (but this might distort the signal shape a lot, depending on the Q-factors and number of filters). Here I show an example with a single filter and a known interferer frequency.

clear all;

Fs = 1000; % 100Hz sampling frequency
t = -5:1/Fs:5;

f_sig = 0.5;

f_intf = 15;
A_intf = 0.05;

% Generate two exemplary signals and pollute them with interference
signal1 = double(sin(2*pi*f_sig*t) &gt; 0);
signal2 = cumsum(signal1 - 0.5) / (Fs/f_sig);
noise1 = A_intf * sin(2*pi*(f_intf)*t);
noise2 = A_intf * sin(2*pi*(f_intf)*t);

received1 = signal1 + noise1;
received2 = signal2 + noise2;

% create a notch filter to filter out this intreference.
w0 = f_intf/(Fs/2);
bw = w0/1; % 1 is the Q-factor here, you can play with it to see the effect
[num, den] = iirnotch(w0, bw);

filtered1 = filter(num, den, received1);
filtered2 = filter(num, den, received2);

subplot(2,2,1);
plot(t, received1)
title('received 1');
grid;

subplot(2,2,2);
plot(t, received2);
title('received2');

subplot(2,2,3);
plot(filtered1);
title('filtered 1');

subplot(2,2,4);
plot(filtered2);
title('filtered2');




Edit: I have ported your code to python, using the correct filter. Here's the filter design function that was taken from scipy v0.19 (not yet officially released), you can find the source here: Scipy Source

def _design_notch_peak_filter(w0, Q, ftype):
    """
    Design notch or peak digital filter.
    Parameters
    ----------
    w0 : float
        Normalized frequency to remove from a signal. It is a
        scalar that must satisfy  ``0 &lt; w0 &lt; 1``, with ``w0 = 1``
        corresponding to half of the sampling frequency.
    Q : float
        Quality factor. Dimensionless parameter that characterizes
        notch filter -3 dB bandwidth ``bw`` relative to its center
        frequency, ``Q = w0/bw``.
    ftype : str
        The type of IIR filter to design:
            - notch filter : ``notch``
            - peak filter  : ``peak``
    Returns
    -------
    b, a : ndarray, ndarray
        Numerator (``b``) and denominator (``a``) polynomials
        of the IIR filter.
    """

    # Guarantee that the inputs are floats
    w0 = float(w0)
    Q = float(Q)

    # Checks if w0 is within the range
    if w0 &gt; 1.0 or w0 &lt; 0.0:
        raise ValueError("w0 should be such that 0 &lt; w0 &lt; 1")

    # Get bandwidth
    bw = w0/Q

    # Normalize inputs
    bw = bw*np.pi
    w0 = w0*np.pi

    # Compute -3dB atenuation
    gb = 1/np.sqrt(2)

    if ftype == "notch":
        # Compute beta: formula 11.3.4 (p.575) from reference [1]
        beta = (np.sqrt(1.0-gb**2.0)/gb)*np.tan(bw/2.0)
    elif ftype == "peak":
        # Compute beta: formula 11.3.19 (p.579) from reference [1]
        beta = (gb/np.sqrt(1.0-gb**2.0))*np.tan(bw/2.0)
    else:
        raise ValueError("Unknown ftype.")

    # Compute gain: formula 11.3.6 (p.575) from reference [1]
    gain = 1.0/(1.0+beta)

    # Compute numerator b and denominator a
    # formulas 11.3.7 (p.575) and 11.3.21 (p.579)
    # from reference [1]
    if ftype == "notch":
        b = gain*np.array([1.0, -2.0*np.cos(w0), 1.0])
    else:
        b = (1.0-gain)*np.array([1.0, 0.0, -1.0])
    a = np.array([1.0, -2.0*gain*np.cos(w0), (2.0*gain-1.0)])

    return b, a


Here's your corrected code using this function:

length = 1200
frequency = .00259
interference = 0.0537
interference_amplitude = 0.2

def iirfilter_filter(data):
        nyq = 1/2.
        bottom = interference/nyq
        top = interference*10/nyq
        #b, a = iirfilter(10, [bottom,top], rs=60, btype="bandstop" , analog=False , ftype="cheby2" )
        b, a = _design_notch_peak_filter(interference*2, 1, 'notch')
        return lfilter(b,a, data)


time_series = np.linspace(1, length, length )

sine_wave = np.sin(2 * np.pi * frequency * time_series )
square_wave = (sine_wave / abs(sine_wave))

sine_wave = sine_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) - 1
square_wave = square_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) + 1

noisy_sine_wave = sine_wave + ( numpy.random.randn(length) * interference_amplitude)
noisy_square_wave = square_wave + ( numpy.random.randn(length) * interference_amplitude)



a = plt.subplot(2,2,1)
b = plt.subplot(2,2,2)
c = plt.subplot(2,2,3)
d = plt.subplot(2,2,4)


a.plot( sine_wave )
a.plot( square_wave )
b.plot( iirfilter_filter( sine_wave ) )
b.plot( iirfilter_filter( square_wave ) )
c.plot( noisy_sine_wave )
c.plot( noisy_square_wave )
d.plot( iirfilter_filter( noisy_sine_wave ) )
d.plot( iirfilter_filter( noisy_square_wave ) )

plt.show()


And the output:


Note that you cannot filter the white noise from randn (corrected from rand), as it contains all frequencies and the notch filter only removes a single frequency.