Convert Accelerometer data to displacement values

[tjakeman]

Hello,

I'm sure this topic has been covered before, but I'm still struggling with it.

The Story:

I have signal data from an accelerometer mounted to a car - one lap around a track. After filtering this data (high pass 13Hz) to remove ridgid body motion, I used it to drive a shaker rig with a component attached. An accelerometer on the conpoent was used to measure the response acceleration.

In order to convert to displacement (which I gather isn't an ideal process) I FFT the response, divide by the samping frequency squared then inverse FFT...?

My displacements are tiny, not really sure why. This is my m-file:

Fs = 512; % Sampling frequency
T = 1/Fs; % Sample time
L = 53760; % Length of signal(no. of data lines)
t = (0:L-1)*T; % Time vector

a = load('TrackReplay Z-axis.txt');

x = a,2);
y = a,3);
z = a,4);

x = x-mean(x*9.81);
x = detrend(x, 'linear');
y = y-mean(y*9.81);
y = detrend(y, 'linear');
z = z-mean(z*9.81);
z = detrend(z, 'linear');

X = fft(x,L)/L; % FFT of x
Y = fft(y,L)/L;
Z = fft(z,L)/L;

Xd = X/(2*pi*Fs)^2;
Yd = Y/(2*pi*Fs)^2;
Zd = Z/(2*pi*Fs)^2;

dx = ifft(Xd);
dy = ifft(Yd);
dz = ifft(Zd);

figure(1)
plot(t,dy)

Any help greatly appreciated

Tom

 

[IRstuff]

Huh? You have to integrate acceleration twice to get displacement; THAT's where the difficulties with the conversion lie.

TTFN

FAQ731-376

 

[tjakeman]

Hello,thanks for your response.

That is initially what I did, but there was a great deal of drift etc. Saw someone on the internet suggesting this approach. I think the logic is:

convert to frequency domain, then:

(1) x=Asin(wt)
(2) v=wAcos(wt)
(3) a=-w^2Asin(wt)

Substitute (1) in (3): a /-w^2 = x

I had misunderstood last night what I used my sampling frequency for w. I think I should have divided by a frequency vector from 0-512 Hz. This gives a sensible plot but the values are all too small.

Have a completely misunderstood this subject?

 

[GregLocock]

No, you are now on the right track.

Bear in mind that your results from double integrating an acceleration spectrum in the frequency domain should match the spectrum of the signal resulting from double integrating the original signal in the time domain, which is easy peasy to program.

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[FeX32]

That is only valid for steady state sinusoidal signals


Fe

 

[tjakeman]

Isn’t that what the FFT’d signal is though, a collection of steady state sine waves?

Well, I thought it seemed too good to be true.

 

[DavidShih]

Values are too small? Did you check the engineering unit conversion? For example, if you convert g to mm by double integration, you should multiply you data by 9800.

 

[IRstuff]

How small are they? After high-pass filtering at 13 Hz, there shouldn't be a whole lot of amplitude left since the car's suspension is low pass filter to start with.

TTFN

FAQ731-376

 

[tjakeman]

Thanks everyone for your help.

By small I mean microns. Its quite a flappy component too, you can see it moving +- 15mm on the shaker rig. First mode is 77 Hz so I should easily have kept the peak displacements.

David, good suggestion but I think my units are consistant.

I was hoping someone might have done this before and tell me straight away that i'm either totally wrong, or point out what I'm doing wrong..

 

[IRstuff]

If you see it moving that much, then the object is pretty much at its first mode, and you should see a sinusoidal time series and a frequency spike at around your first mode; do you see those? If not, then there's something wrong with your data. You should be seeing something on the order of 360-g peak accelerations.

TTFN

FAQ731-376

 

[tjakeman]

Yep, and yes - massive spike at mode one in the frequency domain. Peak g is around 350-400. I believe my data is good, just struggling to convert to displacements! What is your suggested method for this IRstuff?

 

[FeX32]

If you are sure that the data is approx. steady state then your a /-w^2 = x is the way to go.
Just make sure that the accel. is calibrated properly, otherwise you will get wacky numbers


Fe

 

[MikeyP]

FeX32, steady state or not makes no difference the integration algorithm doesn't know that the data has come from a vibrating system.

tjakeman - 3 things to consider

Run a perfect sinusoid through it and see if it gives you the expected result from hand calcs.

When you perform the ifft the resulting numbers should be real (exept for a small round-off error of 1e-15 or so), If there is any imaginary component then you have done something wrong

Check the scaling of the amplitudes following FFT. You are dividing by L. I *think* you should be dividing by L/2 (although it is always best to check you scaling by running a sinusoid through the code).

Here is a link that may be useful

http://www.mathworks.co.uk/support/tech-notes/1700/1702.html

M

--
Dr Michael F Platten

 

[tjakeman]

Hi Mike,

Thanks for your input. Even with a pure sine wave I don’t get the expected results. The imaginary part after ifft is of the same magnitude as the real part. I followed the FFT example you linked me to and did the following:

Fs = 512; % Sampling frequency
T = 1/Fs; % Sample time
L = 4096; % Length of signal(no. of data lines)
t = (0:L-1)*T; % Time vector

% x =Asin(wt)
% v =wAcos(wt)

f = 10;
w = 2*pi*f;
A = 20;
y2=-w^2*A*sin(w*t);

y = detrend(y2, 'constant');

nfft = 2^(nextpow2(length)); % Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
ffty = fft(y,nfft); % Take fft, padding with zeros so that length(fftx) is equal to nfft

NumUniquePts = ceil((nfft+1)/2); % Calculate the number of unique points
ffty = ffty(1:NumUniquePts); % FFT is symmetric, throw away second half

my = abs(ffty)/length; % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
my = my.^2; % Take the square of the magnitude of fft of x.

% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
if rem(nfft, 2) % odd nfft excludes Nyquist point
my(2:end) = my(2:end)*2;
else
my(2:end -1) = my(2:end -1)*2;
end

f = (0:NumUniquePts-1)*Fs/nfft; % This is an evenly spaced frequency vector with NumUniquePts points.

% Generate the plot, title and labels.
figure(1)
subplot(211)
plot(t,y)
subplot(212)
plot(f,my);
title('Power Spectrum of response in y');

w2 = w.^2;
Yd = my./-w2;

dy = ifft(Yd);

t2 = [0:8/(L/2):8];
figure(2)
plot(t2,dy)



I’m not sure if its appropriate to only take the absolute values after FFT in this case, but it doesn’t seem to affect the results either way…

 

[FeX32]

"steady state or not makes no difference the integration algorithm doesn't know that the data has come from a vibrating system"

Of course, but is x=Asin(wt) not a steady state response only?
imo, power spectrum vs. frequency is a different story.

So are we saying that: dy = ifft(Yd) is valid for any signal?


Fe

 

[FeX32]

One last thought,
Check you disp. values by high pass filtering and numerically integrating twice.....usually works nicely, ss or not ss.


Fe

 

[IRstuff]

When you ifft, you should ostensibly multiply by L to undo the scaling for the time domain, otherwise the fft-ifft sequence is asymmetrical.

TTFN

FAQ731-376

 

[GregLocock]

Your 9 april script has at least two errors.

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[GregLocock]

Looping back, the good news is that a simple script agrees that a 350 g acceleration at 77 Hz gives approximately 15mm of deflection, whether I integrate in the time domain or the frequency domain, and that both techniques will work for arbitrary and ugly step function type acceleration signals as well.


Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?