1
$\begingroup$

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.

enter image description here

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:enter image description here

$\endgroup$
  • 1
    $\begingroup$ An extremely simple solution is to apply a simple thresholding operation on the signal, but maybe this is not suitable here (because of some signal properties that we dont know). Additionally, can you plot the fourier transform of the signal? $\endgroup$ – Maximilian Matthé Nov 23 '16 at 13:03
  • $\begingroup$ Do you know anything about the square waves, eg a possible range of frequencies? $\endgroup$ – Marcus Müller Nov 23 '16 at 13:06
  • $\begingroup$ Maximilian: I've added the fft of the original signal. As for thresholding: what would I threshold against? And you are right, this particular signal is a square wave because it's the simplest example I've found, but there are more complexe signals in there. $\endgroup$ – user247243 Nov 24 '16 at 0:18
  • $\begingroup$ Marcus: As mentioned in my comment to Maximillian: I think the square waves are the unknown part (i.e. I've picked a simple example from the data, but some of the other signals are more complex in shape and they are stochastic as well) However, I think there might be hope with the interference noise itself: there's chances that it's always the same frequency* and amplitude (with an added noise component to it). * or at least from a pool of known frequencies. $\endgroup$ – user247243 Nov 24 '16 at 0:20
  • $\begingroup$ I've edited my answer to also contain a python code that works similarly. $\endgroup$ – Maximilian Matthé Nov 25 '16 at 11:00
0
$\begingroup$

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) > 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');

program output

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 < w0 < 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 > 1.0 or w0 < 0.0:
        raise ValueError("w0 should be such that 0 < w0 < 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: python 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.

$\endgroup$
  • $\begingroup$ Thank you. Those output signals look really promising. I will try and reproduce this in python. $\endgroup$ – user247243 Nov 24 '16 at 12:43
  • $\begingroup$ Once you have it veryfied in Python, you might consider closing this question by marking the correct answer. $\endgroup$ – Maximilian Matthé Nov 24 '16 at 18:40
  • $\begingroup$ Amazing. Thank you for putting in this much time! (There is no way I would have figured that out on my own) $\endgroup$ – user247243 Nov 25 '16 at 12:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.