DSPRelated.com
Forums

Deriving velocity from accelerometer data

Started by KousikBarathwaj 5 days ago12 replieslatest reply 4 days ago155 views

Hello, I am currently using an accelerometer for measuring the vibration signatures in rotating machines. But the requirement has changed now, to have the velocity to be compared against the ISO 10816-3 standard. I am collecting 2048 samples at a sampling rate of 10 kHz. So I am trying to derive velocity from the discrete time acceleration data as below. 

a

Where:
u = initial velocity (I am considering initial velocity to be 0 for the first sample)
v = final velocity
a = acceleration
t = time  (0.0001s as my sampling rate is 10kHz).

I expect the velocity output to have an oscillating pattern, similar to the acceleration. But to my surprise, for the given acceleration i got velocity as a ramp of decreasing slope as below. 

image (1)_7661.png image_78984.png

After a bit of googling, i found that the DC offset can cause this behavior. So as a next step i tried to remove the DC offset by high pass filtering in the frequency domain by following the steps as follows. 

1. FFT of the acceleration signal.

2. Apply high pass filter by multiplying with rectangular high pass window (0 from DC-10Hz, 1 thereafter. Basically i made the first two bins 0, given my bin resolution of 5HZ)

3. IFFT of the complex signal and taking the magnitude (Absolute value complex number)

4. Velocity computation as explained above. 

This approach has not solved my problem and now my output looks as follows,

image (2)_3501.pngimage (3)_62778.png

Observation and Queries:

1. The effect of calculating the absolute magnitude has shifted my negative acceleration to the positive side.  How do i reconstruct my acceleration signal properly with oscillating pattern, having the negative terms intact? 

2. Velocity is still a ramp with increasing slope. What am i missing here? 

Any help on this would be highly appreciated. 

Note: 

1. I am just getting started in signal processing and excuse if it seems novice.

2. I am only interested in the RMS magnitude of the velocity signal currently and not the phase of the velocity signal. 

3. I do have come across FIR and IRR for time domain filtering but the approach mentioned above seems more understandable to me. 

Best Regards,

Kousik K





[ - ]
Reply by djmaguireJuly 9, 2024
You don't remove DC with a lowpass filter.  You remove DC with a highpass filter. 
[ - ]
Reply by KousikBarathwajJuly 9, 2024

My bad. I actually meant High pass filter. 

[ - ]
Reply by djmaguireJuly 9, 2024
Got it to the mistype, Kousik.  Admittedly, I didn't read much beyond that in case it wasn't a mistype.

I know that you mentioned being more comfortable with the frequency domain filtering, but time domain filtering can be a bit more intuitive with respect to filter characteristics such as settling time.  That can impact subsequent operations like integration.  I'm not suggesting that the leads that folks are providing are not spot on.  It is just food for thought.


[ - ]
Reply by LKoppJuly 9, 2024

Hello Kousik,

It seems to me that you are experiencing what is called as "Brownian" motion. Integration white noise (zero mean) should provide an estimate with a variance proportional to time whatever you do..  but is there is a bias (drift or whatever) the problem is worse because you are integrating a deterministic phenomenon..

Also this recalls me of what is called the "Allan Variance" that allow to analyse the variance of flicker noise (Clock fluctuations)..You should find easily links on that topic..

Rgards

Laurent

[ - ]
Reply by KousikBarathwajJuly 9, 2024

Thanks for your guidance. Will look into it. 


[ - ]
Reply by hjnitzponJuly 9, 2024

You may do your integration in the time or in the frequency domain. If you do it in the time domain, a simple method is the trapezoidal integration (the example below uses multiple axes), here implemented as an IIR filter:

# ----- integration ----- #

def acc2vel(self):

coef = 1 / (2 * self.sample_rate)

b = [coef, coef]

a = [1, -1]

x = np.asmatrix(self.df)

x = x[:, 1:4]

x = signal.detrend(x, axis=0, type='constant')

x = signal.lfilter(b, a, x, axis=0)

x = signal.detrend(x, axis=0, type='constant')

self.df.iloc[:, 1:4] = x

It is important to eliminate DC offsets and eventually low frequency drifts. Either using a simple detrending or you may use a high-pass filter.


Alternatively, you may do the integration in the frequency domain:

1) FFT

2) divide complex bins by: FFT./(2*pi*f) (f = frequency vector of your FFT)

3) IFFT


[ - ]
Reply by ombzJuly 9, 2024

Hi. Even you mentioned in another reply that you mistyped low-pass, I believe that you actually did low-pass filter your acceleration data. It seems that on your first velocity plot that some small level oscillation overlaps the strong drift. On your second velocity plot, the signal seems more smooth and those oscillations seem gone (which somewhat points me to believe that you integrate a smoothed bias value, and not oscillations+noise).

I don't know what kind of FFT implementation you use, but be aware that the center of a usual FFT-output corresponds to f_s/2 and not frequency 0. Or in other words: I think you mistakenly take the wrong FFT-bins into account. Btw: in your FFT domain you should see a small peak (or bump) corresponding to your oscillations: you want to inverse-FFT the bins making up that peak (or bump) by zeroing out the other bins. But then: be aware of introducing oscillations in your filtered acceleration signal by using a rectangular window in the Fourier domain... But I think you get the idea.

[ - ]
Reply by PlatybelJuly 9, 2024

Hello Kousik:

Your original result is probably correct. Think of velocity and acceleration as related through integration and differentiation.  Both these operations affect the spectrum of the time series.  Differentiation acts as a high pass filter, highlighting the high frequencies.  

Your expression of converting the acceleration to velocity is, in effect, integration. You are, in effect, applying a low pass filter.  In your velocity output the low frequency is dominant (high magnitude); the high frequencies are reduced in magnitude, but not totally eliminated.  You can see the high frequencies in the output by enlarging the vertical scale of the plot.  The fluctuations that you see in the acceleration data are still there in the velocity output, but not so visible, as the plot is now dominated by high magnitude low frequencies.

Lalu Mansinha


[ - ]
Reply by njswedeJuly 9, 2024

1. Like others have pointed out, it looks like you lowpass filtered the data, not highpass.

2. A simple DC blocker can be realized by calculating the average and subtracting it from the signal. Since you seem to have the 2048 samples all at once, you don't have to bother with streaming implementations of averaging etc. 

[ - ]
Reply by jbrowerJuly 9, 2024

Kousik-

First very well presented question and attention to detail. You may be a DSP novice but you will go far.

Second, I concur with DJ Maguire, to debug this you might try some time domain things. If dc offset is really the problem, experiment with subtracting sum/2048 from your input and see what happens. That's a simple step and eliminates the evident confusion about whether you high-passed or low-passed your data. If that works or is an improvement then you can think about next steps.

-Jeff

[ - ]
Reply by Mannai_MuraliJuly 10, 2024

I think after IFFT you have taken absolute value for acceleration and has fed the absolute value to the velocity calculation.Actually you should take real part and give to velocity calculator.If you feed the absolute value all negative acceleration will become positive.From your graph it seems the acceleration has + ve and -ve going ramp when integrated it become a parobola!

Also make sure in FFT you set 0 not only +ve frequency but also -ve frequencies.ie +k and N-k.

[ - ]
Reply by Kenneth43July 10, 2024

Hi Kousik

If you have an accelerometer on the nonrotating part of a machine with a rotating part, the way to find out the rotational speed is to find a repeated pattern in the accelerometer signal. The rotational frequency is the same as the repetition frequency.