Removal of noise from vibration recordings

Question

Problem
I am analyzing accelerometer recordings that are contaminated by 'spikes' caused by issues in the digital circuitry. The spikes are problematic, because I wish to determine the peak-amplitude from these signals using an automated procedure from a large number (hundreds, perhaps thousands) of these recordings. Hence, manually determining the amplitudes is possible, but not preferred.

Background
The acceleration recordings are oscillatory in nature and reach frequencies of ~200 Hz. I am planning on using a low-pass filter in Matlab to remove the contaminating spikes in the signal. Sample rate is around 730 Hz. An FFT-based filter may be complicated by the fact that the signal gradually increases and decreases in frequency over time.

Here is a picture of an example recording including the pesky spikes (the noise) in the first 100 ms:

Question
I am using Matlab for signal processing and I am looking for advice on which digital filter for use as a high-pass filter? For example, Butterworth and Chebyshev filters are quite commonly applied if I am not mistaken, but why? Will wavelets do me any good?

Answer 1

If you want to detect the peak value of the time-domain signal, it is important not to change phase of each frequency signals. Otherwise, the peak value will be affected by it and will slightly change.

I think you already have recorded signal and you can process it after the recording. If so, it is better to use zero-phased filters, like this function in matlab.

http://uk.mathworks.com/help/signal/ref/filtfilt.html

In this case, the needed function is to reduce the spike, which contains so much frequencies. In the help page I have shown before, there is an example of filtering this kind of spikes from the relatively low-frequency signal. I think it's what you need.

Edit

I recommend Butterworth low pass filter for its flat gain characteristics. The gain ripple of the Chebtshev would affect badly. The phase change will be corrected by using the filter backward, which the filtfit does. As you imagine, the effect of the phase change would be very little in this case, but it cannot be guaranteed. Just for your information, Bessel would be a candidate when you do the filtering online.

If you still doubt about the effect of the phase shift (not the time-shift) by the filter, you can imagine the time-streched pulse and impulse signal in the image and code below. These signals are almost the same but only the phase are different. Besides the amount of the effect, the phase shift does affect to the amplitude of the signal.

I recomend you to use zero-phase filter if you want to measure the actual amplitude precisely.

clear all;
close all;

% generate signal with the same phase
N = 2^15; % number of samples for this signal
for k = 1:N/2;
    S(k+1) = exp(0);
end
S(N/2+2: N) = conj(S(N/2: -1 :2));

signal1 = real(ifft(S));

% generate signal with different phase from above
% this signal is called time-streched pulse.
N = 2^15;
for k = 1:N/2;
    S(k+1) = exp(-1j*2*pi*N/2*(k/N)^2);
end
S(N/2+2: N) = conj(S(N/2: -1 :2));

signal2 = real(ifft(S));

% and let's see them.
subplot(211)
plot(signal1)
xlim([-1e4, length(signal1)])
xlabel('time')
ylabel('all the same phase')
subplot(212)
plot(signal2)
xlim([-1e4, length(signal1)])
xlabel('time')
ylabel('only the phases are different form above')

Answer 2

As I said in the comment, a high pass filter is not what you want.

Below is some scilab code that generates the plot. The plot shows four signals:

• The original signal.
• The signal high-pass filtered.
• The signal median filter.
• The signal low-pass filtered.

The median filter is a non-linear filter sometimes used in image processing to remove spikes. However, because it is non-linear it may have undesirable effects on your signal.

The low-pass filter still gets impacted by the spike a bit, but it may be preferable to the median filter.

---

CODE:

    x = rand(1,1000,"norm");
    omega = 2*%pi*0.0789234;
    y = filter(1,[1 -2*cos(omega) 1],x);

    mx = max(y);

    z = y;
    z(100) = 3*mx;
    clf
    subplot(411)
    plot(z)
    title("ORIGINAL SIGNAL")

    hp=eqfir(33,[0 .2;.25 .35;.4 .5],[0 1 1],[1 1 1]);

    zhp = filter(hp,1,z);
    subplot(412);
    plot(zhp);
    title("HIGH PASS FILTERED")

    M = 3;
    for k=1:1000,
        lw = max(1,k-M);
        hgh = min(1000,k+M);
        zmdn(k) = median(z(lw:hgh))
    end
    subplot(413);
    plot(zmdn);
    title("MEDIAN FILTERED")

    zlp = filter(0.01,[1 -0.99],z);
    subplot(414);
    plot(zlp);
    title("LOW PASS FILTERED")