1
$\begingroup$

Background: I was trying to remove the DC component of a load cell signal using a high-pass filter. I don't care about that anymore. I'm curious why a digital filter and an analog filter behave differently under the same conditions.

Here (link 1, link 2) is the data I'm working with, and below is my Python code:

import numpy as np
import math as m
from scipy import signal
from scipy import fftpack
from myLib import myFFT
import matplotlib.pyplot as plt

# load data

filename = "raw_data.txt"
fSampling = 1E3
fNyquist = 0.5*fSampling
data = np.loadtxt(filename,skiprows=0)
i1 = 1000
i2 = 5000
t = data[:,0]
x = data[:,3]


# remove DC component

fCutoff = 0.01/fNyquist
bd, ad = signal.butter(4, fCutoff, 'high', analog=False)
xd = signal.filtfilt(bd,ad,x)
ba, aa = signal.butter(4, fCutoff, 'high', analog=True)
xa = signal.filtfilt(ba,aa,x)


# plot results

f, axarr = plt.subplots(3, sharex=True)
axarr[0].plot(t, x)
axarr[0].set_title("raw")
axarr[0].set_ylabel("x [-]")
axarr[1].plot(t, xd)
axarr[1].set_title("digital")
axarr[1].set_ylabel("x [-]")
axarr[2].plot(t, xa)
axarr[2].set_title("analog")
axarr[2].set_xlabel("t [s]")
axarr[2].set_ylabel("x [-]")

plt.show()

Problem The digital filter doubles the values, while the analog filter does not. all plots

Why does the digital filter behave differently? Is something weird about my data, or do digital and analog filters behave differently? I'm new to signals processing (this all originated from me leaving the analog flag with the default "False"), but I thought analog and digital filters were generally different means to the same end.

Edit: I did some research online, and I've learned that I should use a digital filter to filter regularly sampled data. I don't know why yet (digital filters perform better at a known range?), but this is just making even more curious why my digital filter increased the gain of my signal. I analyzed the frequency response of the digital filter, and theoretically I should get a gain of 0: freq response of digital HPF

I want to chalk this up to numerical errors since the digital filter's numerator and denominator are so close, but so are the analog filter's numerator and denominator: 

b digital = [ 0.99991791 -3.99967164  5.99950746 -3.99967164  0.99991791]
a digital = [ 1.         -3.99983581  5.99950745 -3.99950746  0.99983583]
b analog = [ 1.  0.  0.  0.  0.]
a analog = [  1.00000000e+00   5.22625186e-05   1.36568542e-09   2.09050074e-14
   1.60000000e-19]
|improve this question
$\endgroup$
  • $\begingroup$ Can you please upload the sampled data file in matlab readable format? $\endgroup$ – Fat32 Jul 12 '17 at 19:06
  • $\begingroup$ I already have. The "here" in "here is the "raw raw" data..." is a hyperlink to pastebin where you can download a tab-delimited copy of the data. $\endgroup$ – techSultan Jul 12 '17 at 19:10
  • $\begingroup$ sorry I havent seen it, but I have problems of downloading from there, can you upload it into another web space provider ? $\endgroup$ – Fat32 Jul 12 '17 at 19:32
  • $\begingroup$ Here you go: zerobin.net/… $\endgroup$ – techSultan Jul 12 '17 at 19:35
  • $\begingroup$ ok I have downloaded the data and selected the 4th column. But its y-value is not about 63 as in your plot but about 0.428, so you have scaled you data ? Also a DC notch may not be the best option here... because the initial result is not satisfactory. $\endgroup$ – Fat32 Jul 12 '17 at 19:54
1
$\begingroup$

I was trying to remove the DC component of a load cell signal

The simplest way would be x = x - mean(x), but this won't line up perfectly if you're processing multiple consecutive chunks.

Problem The digital filter doubles the values,

Your 4th-order filter is too aggressive, and one pole has moved outside the unit circle due to numerical error:

unit circle

Zoomed-in:

zoomed in to DC

So the behavior is not going to be sensible. Try using cascaded second-order sections instead:

sos = signal.butter(4, fCutoff, 'high', analog=False, output='sos')
xd = signal.sosfiltfilt(sos, x)

while the analog filter does not.

You can't do that. filtfilt() is a digital filter. You can't feed it analog filter coefficients and expect it to do anything reasonable.

The only purpose of analog=True in butter() is to then transform them into digital filters, or to calculate component values for analog circuits. So:

ba, aa = signal.butter(4, fCutoff, 'high', analog=True)
xa = signal.filtfilt(ba, aa, x)

is nonsense

|improve this answer
$\endgroup$
  • 1
    $\begingroup$ This is exactly what I was looking for. I tried subtracting the mean, but I have to apply it piece-wise because of hysteresis in the load cell. I never knew about filtfilt() being for digital signals only. Ditto for "analog" filters in butter(). So is designing an analog filter on a computer a moot point? $\endgroup$ – techSultan Jul 13 '17 at 19:06
  • 1
    $\begingroup$ @techSultan Designing an analog filter is useful, but you can't apply it to a digital signal. You know that filtfilt() applies the filter twice, forward and backward in time? I don't think it will line up correctly on the boundaries of chunks, either. If processing a signal in chunks you'd typically use lfilter (or sosfilt in this case) and keep the zi/zf initial conditions from one chunk to the next. $\endgroup$ – endolith Jul 13 '17 at 19:11
  • $\begingroup$ Say I subtract the mean values in my signal, such as x= x - simple moving average (x) or piecewise x = x - mean(x). Why can't I apply A digital LPF via filtfilt() to remove the remaining high frequency noise? How would the forward/backward pass affect the boundaries? $\endgroup$ – techSultan Jul 13 '17 at 19:18
  • 1
    $\begingroup$ @techSultan IIR digital filtering always has a "tail" that extends infinitely (hence the name), outside of the chunk. If you use initial conditions from previous chunk, you are carrying the tail from the previous chunk into the current one, and so on, so it's correct. For this to work with filtfilt, you'd need to do this an infinite number of times forward and backward, which is of course impossible. Depending on the tail, you could truncate it after some time, though, where it won't be noticeable. $\endgroup$ – endolith Jul 13 '17 at 20:07
  • 1
    $\begingroup$ matlab has a function called detrend Remove a linear trend from a vector, usually for FFT processing. Y = detrend(X) removes the best straight-line fit linear trend from the data in vector X and returns the residual in vector Y. If X is a matrix, detrend removes the trend from each column of the matrix. Y = detrend(X,'constant') removes just the mean value from the vector X, Y = detrend(X,'linear',BP) removes a continuous, piecewise linear trend. Breakpoint indices for the linear trend are contained in the vector BP. $\endgroup$ – Stanley Pawlukiewicz Jul 13 '17 at 21:29
1
$\begingroup$

The following coefficients $a$ and $b$ implement a second order digital DC notch filter : 

b = 0.9900   -1.9800    0.9900

a = 1.0000   -1.9800    0.9801

Whose frequency response is (y-axis in dB):

enter image description here

Given the raw-input signal $x[n]$ as a vector of samples in $x$ in matlab, you would apply the following line to filter it:

y = filter(b,a,x);

Note that the filter is not linear phase, hence you will have some phase distorions especially about DC frequencies... if you need a strictly linear phase DC blocking filter please indicate.

Given the raw data with the following plot: enter image description here

The resulting filtered signal plot is the following: enter image description here

It can be seen that initially there is a transient starting from a value of $0.428$ and takes some samples to stabilize to remove the DC. This transient is a direct consequence of long impulse response associated with the IIR filter.

Shorter impulse responses would provide faster response, but more steady state distortions.

Assuming that a shorter filter would do better, here is another set of DC notch coefficients with shorter impulse response:

b = 0.9100   -1.8200    0.9100
a = 1.0000   -1.8200    0.8281

In which case the output becomes: enter image description here

Clearly an improvement in dynamic response, but keep in mind that the steadystate response now gives more distortions.

You can even further get a quicker response with the following set:

b = 0.7000   -1.4000    0.7000
a = 1.0000   -1.4000    0.4900

whose result is: enter image description here

Much faster response, but there could be untolerable distortions on the samples, depending on the application requirements.

|improve this answer
$\endgroup$

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.