3
$\begingroup$

I'm pretty well versed in statistics, but not really digital signal filtering. I have a data scenario where I expected to be able to pretty easily filter out some noise (human pulse) that's at a known frequency band, but I'm having a lot of trouble using the standard tools in the scipy.signal library and think I must be misunderstanding how to design digital filters. I have a notebook here that walks through my explorations thus far, but the gist is that the standard scipy filters seem to cause large distortions at the start and end of my signal, with the precise behaviour dependent on the phase of the noise signal I'm trying to subtract. Just in case the above binder link goes down, I'll include some of the key points below as well:

First generating some synthetic data that's similar to my real data:

#generate time vector
samples_per_sec = 10.0
total_time = 40.0
time = np.linspace(0, total_time, int(total_time*samples_per_sec))

#generate the pulse signal
pulse_hz = 1.0
pulse_phase = np.radians(0)
pulse = np.sin(time*(2*np.pi)*pulse_hz - pulse_phase)

#generate the BOLD signal (just something that goes up then down)
dist = stats.beta(2, 2)
bold = dist.pdf((time-10)/20) / 10.0 # division by 10 to make bold a small signal

#combine
pulse_plus_bold = pulse+bold
plt.plot(time, pulse_plus_bold);

combined signals

Try a 1st order butterworth:

#1st order butterworth filter in ba mode
ba1 = signal.butter(
    output = 'ba'
    , N = 1 #needs to be low if using output='ba', else use output='sos' and sosfiltfilt
    , Wn = [0.5,1.5]
    , btype = 'bandstop'
    , fs = samples_per_sec
)
filtered_ba1_nopad = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_plus_bold
    , padtype = None
)
plt.plot(time, filtered_ba1_nopad, 'b');
plt.plot(time, bold, 'r--');
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

enter image description here

First-order butterworth with even padding:

filtered_ba1_pad_even = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_plus_bold
    , method = 'pad'
    , padtype = 'even'
)
plt.plot(time, filtered_ba1_pad_even, 'b');
plt.plot(time, bold, 'r--');
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

enter image description here

First-order butterworth with odd padding:

filtered_ba1_pad_odd = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_plus_bold
    , method = 'pad'
    , padtype = 'odd'
)
plt.plot(time, filtered_ba1_pad_odd, 'b');
plt.plot(time, bold, 'r--');
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

enter image description here

This latter looks really good! But after playing around I discovered that whether odd or even (or either) padding works better seems to be contingent on the phase of the signal being filtered-out. As an example, while the above obtained excellent filtering with odd-padding, here's the same scenario but with a phase-shift added to the pulse signal that yields edge artifacts in both odd and even:

phase = np.radians(45)
pulse_shifted = np.sin(time*(2*np.pi)*pulse_hz - phase)
pulse_shifted_plus_bold = pulse_shifted+bold

filtered_shifted_ba1_pad_odd = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_shifted_plus_bold
    , method = 'pad'
    , padtype = 'odd'
)
filtered_shifted_ba1_pad_even = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_shifted_plus_bold
    , method = 'pad'
    , padtype = 'even'
)

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(5, 3))
axes[0].plot(time, filtered_shifted_ba1_pad_odd, 'b')
axes[0].plot(time, bold, 'r--')
axes[1].plot(time, filtered_shifted_ba1_pad_even, 'b')
axes[1].plot(time, bold, 'r--')
fig.tight_layout()
plt.title('Odd (left) and Even (right)')
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

enter image description here

I've also tried the 'gust' padding method as well as higher order filters (using sos of course), and observe the same phase-dependent edge artifacts in everything I've tried. Any tips?

$\endgroup$
  • $\begingroup$ Is your real signal going to be continuous, or will you be using short chunks like you are here? $\endgroup$ – TimWescott Aug 6 at 16:04
  • $\begingroup$ @TimWescott, thanks for your interest. Short chunks as I have here. $\endgroup$ – Mike Lawrence Aug 6 at 16:46
3
$\begingroup$

Your basic problem is that filtfilt (and most other linear filtering routines) take filters that are designed for infinitely long time expanses, and apply them to a chunk of data as if the data were extended infinitely in both directions with zeros.

So you have a legitimate bandpass filter, and it's "seeing" a legitimate jump in the signal at the starting point of your signal.

There's three basic approaches you can take; the first two are ad-hoc and easy, the third one is difficult if you're starting from first principles. It's certainly been solved out there someplace, but a brief search here on "filter finite-length data" didn't find me joy.

Approach 1: window the input data

Take your input data, and multiply it by something that'll make it taper off at the ends. I.e. a ramp from 0 to 1 over 10 samples at each end, or $\frac{1}{2}\left (1 - \cos \frac{\pi n}{N} \right)$ for N samples at each end (suitably reversed on the left end). You'll have some artifacts (a rising sine wave isn't the same as a steady one, after all), but they'll be attenuated. Here's python code implementing a cosine edge attenuation with the ability to customize what central % of the signal is kept as 1:

def attenuate_edges(signal,time,edge_attenuation_percent):
  start = int(np.floor(len(time)*edge_attenuation_percent))
  end = int(len(time)-start)
  ramp = (1-np.cos(np.pi*(np.arange(start)/start)))/2
  edge_attenuator = np.ones(len(time))
  edge_attenuator[0:start] = ramp
  edge_attenuator[end:len(time)] = np.flip(ramp)
  return(signal*edge_attenuator)

Approach 2: Trim the output data

Do what you're doing now, and lop off the nastiness at the ends. This is probably the easiest, and if you can just collect a bit more data, doesn't lose you anything.

Approach 3: Do a proper estimate of the interfering signal, and subtract it out

This will be fun-fun if you love math and have the time. Basically you'll use the fact that the value of your interfering signal at time $n$ correlates in a specific way with the values of your interfering signal at time $k$ for all values of $n$ and $k$ in your data set. You'll probably end up with something that looks a lot like a Wiener or a Kalman filter, that takes the end-effects into account. Your estimate will be worse at the ends, but this will show up as a bit of noise on the ends -- not as those honkin' big pulses.

If I couldn't figure out the search terms for this, it would take me a day to do and another day to verify, and supposedly an expert. OTOH, Gauss or Laplace probably invented it in the 19th century, and may even have thought it important enough to write down, somewhere. So I'm sure the method exists.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Further searching for estimating signals on finite-duration samples didn't find anything relevant, either. The blessing of signal processing is that you can almost always assume signals of infinite duration and linear filters. The curse of signal processing is that when those assumptions don't hold, you either have to dig out some really obscure prior research, or you have to invent your own method -- knowing that Gauss or Laplace beat you to it, under a different name. $\endgroup$ – TimWescott Aug 6 at 17:46
  • $\begingroup$ Thanks! I added code to your answer for folks that come here later. Seems to work really well. I'm surprised they don't talk about both the edge distortion and solutions thereto like this in the scipy docs. $\endgroup$ – Mike Lawrence Aug 7 at 15:24

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.