Combining Bandpass filters on 1D signal to obtain better SNR? I have a signal which I'm trying to extract information from (specifically, features to recognize using DTW). Parts of it have very clear and unique features which are unfortunately only sometimes superimposed with a higher frequency low power interference. Below is an example of a square wave feature.



In this picture,


Black line represents the original signal.
Blue line is a lowpass Butterworth filter (I have added 3 to the value so that the graph looks more legible). 
Red line is the signal minus the result of a highpass filter (again offset by +5). And the 
green line is a simple box smoothing filter (specifically [1,3,5,7,9,11,9,7,5,3,1]) on the original signal.


The blue line in the second graph is the high pass filter's result.

My questions are: is there any way to disentangle what appears to be a very crisp square wave from a very periodic hum? This hum isn't mains hum, by the way, and it isn't always the same strength or frequency. When it comes on, it'll remain constant all along, but the next time it comes on, it might be a different frequency.

I'm a bit disappointed that I can't get a result that's better than my box filter without severely distorting the resulting signal into something that barely resembles a square signal. Already the low pass filter result (blue graph) is starting to move well past square.

I've gone so far as to thinking of picking the dominant high end frequency, and creating an artificial sine wave to try and "manually remove"  the hum. Is this completely misguided? And what hope do I have of getting a clear signal if I won't have the luxury of fine tuning my filters at design time?

I appreciate any general recommendations. 



Edit 2: 

Here's a functional equivalent to my problem in python. Any comments as to why I'm not obtaining a clean output signal? I've played around with the cutoff values but I'm not very familiar with signal processing so I'm not sure if I'm approaching it right.

length = 1200
frequency = .00259
interference = 0.0237
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" )
        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.rand(length) * interference_amplitude)
noisy_square_wave = square_wave + ( numpy.random.rand(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()


Output:
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.