28
$\begingroup$

I am fairly new to DSP, and have done some research on possible filters for smoothing accelerometer data in python. An example of the type of data Ill be experiencing can be seen in the following image:

Accelerometer data

Essentially, I am looking for advice as to smooth this data to eventually convert it into velocity and displacement. I understand that accelerometers from mobile phones are extremely noisy.

I dont think I can use a Kalman filter at the moment because I cant get hold of the device to reference the noise produced by the data (I read that its essential to place the device flat and find the amount of noise from those readings?)

FFT has produced some interesting results. One of my attempts was to FFT the acceleration signal, then render low frequencies to have a absolute FFT value of 0. Then I used omega arithmetic and inverse FFT to gain a plot for velocity. The results were as follows:

Filtered signal

Is this a good way to go about things? I am trying to remove the overall noisy nature of the signal but obvious peaks such as at around 80 seconds need to be identified.

I have also tired using a low pass filter on the original accelerometer data, which has done a great job of smoothing it, but I'm not really sure where to go from here. Any guidance on where to go from here would be really helpful!

EDIT: A little bit of code:

for i in range(len(fz)): 
    testing = (abs(Sz[i]))/Nz

    if fz[i] < 0.05:
        Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)
Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

So essentially, ive performed a FFT on my accelerometer data, giving Sz, filtered high frequencies out using a simple brick wall filter (I know its not ideal). Then ive use omega arithmetic on the FFT of the data. Also thanks very much to datageist for adding my images into my post :)

$\endgroup$
  • $\begingroup$ Welcome to DSP! Is the red curve in your second picture a "smoothed" version of the original (green) data? $\endgroup$ – Phonon Aug 2 '12 at 19:41
  • $\begingroup$ The red curve is (hopefully!) a velocity curve generated from fft followed by filtering, followed by omega arithmetic (dividing by 2*pifj), following by inv. fft $\endgroup$ – Michael M Aug 2 '12 at 20:24
  • 1
    $\begingroup$ Perhaps if you include a more precise mathematical expression or pseudocode for what you did would clear things up a bit. $\endgroup$ – Phonon Aug 2 '12 at 20:53
  • $\begingroup$ Added some now, thats the general feel of the code.. $\endgroup$ – Michael M Aug 4 '12 at 10:55
  • 1
    $\begingroup$ My question would be: what do you expect to see in the data? You won't know whether you have a good approach unless you know something about the underlying signal that you expect to see after filtering. In addition, the code that you showed is confusing. Although you don't show the initialization of the fz array, it appears that you're applying a highpass filter instead. $\endgroup$ – Jason R Aug 6 '12 at 14:13
13
+200
$\begingroup$

As pointed out by @JohnRobertson in this post, Total Variaton (TV) denoising is another good alternative if your signal is piece-wise constant. This may be the case for the accelerometer data, if your signal keeps varying between different plateaux.

Below is a Matlab code that performs TV denoising in such a signal. The code is based on this paper. The parameters $\mu$ and $\rho$ have to be adjusted according to the noise level and signal characteristics.

If $y$ is the noisy signal and $x$ is the signal to be estimated, the function to be minimized is $\mu\|{x-y}\|^2+\|{Dx}\|_1$, where $D$ is the finite differences operator.

function denoise()

f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);
fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);
rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
    f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
    v = D(f)+(1/rho)*y;
    u = max(abs(v)-1/rho,0).*sign(v);
    y = y - rho*(u-D(f));
end

function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Results:

enter image description here enter image description here enter image description here

$\endgroup$
  • $\begingroup$ Really like this answer, gonna go ahead and try it. Sorry it took me so long to reply! $\endgroup$ – Michael M Aug 28 '12 at 16:51
  • $\begingroup$ Excellent answer. Thanks for the details. I am looking for the C version of this code. Anyone here ported this matlab code to C they would like to share? Thanks. $\endgroup$ – pixbroker Jan 12 '16 at 18:24
  • $\begingroup$ What does piece-wise constant mean? $\endgroup$ – tilaprimera Jun 1 '18 at 20:25
6
$\begingroup$

The problem is that your noise has a flat spectrum. If you assume white Gaussian noise (which turns out to be a good assumption) its power spectrum density is constant. Roughly speaking, it means that your noise contains all frequencies. That's why any frequency approach, e.g. DFT or low-pass filters, is not a good one. What would be your cut-off frequencies since your noise is all over the spectrum?

One answer to this question is the Wiener filter, which requires knowledge of the statistics of your noise and your desired signal. Basically, the noisy signal (signal + noise) is attenuated over the frequencies where the noise is expected to be grater than your signal, and it is amplified where your signal is expected be grater than your noise.

However, I would suggest more modern approaches that use non-linear processing, for example wavelet denoising. These methods provide excellent results. Basically, the noisy signal is first decomposed into wavelets and then small coefficients are zeroed. This approach works (and DFT doesn't) because of the multi-resolution nature of wavelets. That is, the signal is processed separately in frequency bands defined by the wavelet transform.

In MATLAB, type 'wavemenu' and then 'SWT denoising 1-D'. Then 'File', 'Example Analysis', 'Noisy signals', 'with Haar at level 5, Noisy blocks'. This example uses Haar wavelet, which should work fine for your problem.

I'm not good at Python, but I believe you can find some NumPy packages which perform Haar wavelet denoising.

$\endgroup$
  • 1
    $\begingroup$ I would disagree with your first statement. You're assuming that the signal of interest covers the full bandwidth of the input sequence, which is unlikely. It is still possible to obtain improved signal-to-noise ratio using linear filtering in this case, eliminating the out-of-band noise. If the signal is highly oversampled, then you may obtain a large improvement with such a simple approach. $\endgroup$ – Jason R Aug 7 '12 at 1:37
  • $\begingroup$ It's true, and this is achieved by the Wiener filter, when you know the statistics of your signal and your noise. $\endgroup$ – Daniel R. Pipa Aug 7 '12 at 10:57
  • $\begingroup$ Although the theory behind wavelet denoising is complicated, the implementation is as simple as the approach you described. It involves only filter banks and thresholding. $\endgroup$ – Daniel R. Pipa Aug 7 '12 at 12:08
  • 1
    $\begingroup$ Im researching into this now, will post my progress above, thanks to both of you and Phonon for all of your help so far! $\endgroup$ – Michael M Aug 7 '12 at 14:11
  • $\begingroup$ @DanielPipa I don't have access to the matlab packages in question. Can you provide a paper or other reference that describes the method that corresponds to your matlab code. $\endgroup$ – John Robertson Aug 8 '12 at 22:44
0
$\begingroup$

As per suggestion of Daniel Pipa, I took a look at wavelet denoising and found this excellent article by Francisco Blanco-Silva.

Here I have modified his Python code for image processing to work with 2D (accelerometer) rather than 3D (image) data.

Note, the threshold is "made up" for the soft-thresholding in Francisco's example. Consider this and modify for your application.

def wavelet_denoise(data, wavelet, noise_sigma):
    '''Filter accelerometer data using wavelet denoising

    Modification of F. Blanco-Silva's code at: https://goo.gl/gOQwy5
    '''
    import numpy
    import scipy
    import pywt

    wavelet = pywt.Wavelet(wavelet)
    levels  = min(15, (numpy.floor(numpy.log2(data.shape[0]))).astype(int))

    # Francisco's code used wavedec2 for image data
    wavelet_coeffs = pywt.wavedec(data, wavelet, level=levels)
    threshold = noise_sigma*numpy.sqrt(2*numpy.log2(data.size))

    new_wavelet_coeffs = map(lambda x: pywt.threshold(x, threshold, mode='soft'),
                             wavelet_coeffs)

    return pywt.waverec(list(new_wavelet_coeffs), wavelet)

Where:

  • wavelet - string name of wavelet form to be used (see pywt.wavelist(), e.g. 'haar')
  • noise_sigma - standard deviation of noise from data
  • data - array of values to filter (e.g. x, y, or z axis data)
$\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.