How to get mobility from velocity and force

Dear all,

I have a data of velocity (get from accelerometer) and force data (get from hammer). Now, I need to get mobility from them. I already know that mobility=velocity/force. But I don't know how code it. I just begin learning about signal processing for a week.

Thank you all in advance.

Answer 1

[Joe]

I have a background in structural dynamics, not soil mechanics, but this test looks similar enough. You have measurements of acceleration (NOT velocity, unless your accelerometer is doing something it shouldn't) and force in the time domain. You want to get a mobility measurement in the frequency domain. The basic steps you have to follow are:

1) Convert your time domain measurements to the frequency domain

2) Determine the accelerance of your system (accelerance = acceleration/force)

3) Convert from accelerance to mobility

In more detail:

1) Use FFT to go from time domain to frequency domain. MATLAB's default FFT function is a little annoying to use if your background is not signal processing; you have to handle scaling yourself, and you are only going to want the first half of the FFT measurements. The other half are redundant if your signal is real-valued. Do something like this:

 function [X, freqs] = easy_fft(x, fs)
  L = length(x);
  NFFT = 2^nextpow2(length(x));
  X = fft(x, NFFT);
  X = 2/L*X(1:NFFT/2 + 1); % 2/L is the required scaling
  freqs = fs/2*linspace(0,1,NFFT/2+1); % Frequency vector
 end

Then do

 fs = 1/y(2, 1); % Convert time step to sample rate
 [A, freq] = easy_fft(y(:, 2), fs); % Acceleration
 [F, ~] = easy_fft(y(:, 3), fs); % Force

Now you have complex frequency domain vectors A and F. Don't take the real part of anything, they must be complex.

2) Now get the accelerance. We know that A = HF and want H; in your case the simple element-wise division

 H = A./F;

will work. (This ONLY works for single-input, single-output systems such as what you have here). H is now your accelerance.

3) The conversion to mobility is just a scaling factor. Divide the accelerance by i*omega to get the mobility:

 omega = 2*pi*freq; % Convert to rad/s
 mobility = H./(1i*omega); % Acceleration = i*omega*Velocity for a periodic steady-state response

To replicate the plot you showed, take the absolute value of mobility and put it on a log scale vs. frequency.

 semilogy(freq, abs(mobility));

Apologies if there are any bugs in what I wrote, I didn't test it on anything. Hopefully you get the idea.

[vu ngothanh]

Dear Mr. Joe,

Thank you very much for your time. I did as you said. But there is some thing wrong.

In the result, we have

H is a matrix 65537x1 complex double

omega is a matrix 1x65537 complex double.

Therefore, it can not calculate mobility=H./(1i*omega)

So, I transpose omega into a matrix 65537x1. then matlab can calculate H. In your opinion, is it right or wrong. This is my result.

[Joe]

The row vector/column vector issue is due to the way I created the frequency vector in the function I gave you. You're free to transpose it without any issue, as you did.

That frequency response looks pretty standard to me. You can see that after about 5 kHz the readings are garbage, which is not surprising. If you zoom in on the region from 0 to 500 Hz, it looks like there is a response peak with pretty high damping, like you're expecting. You can get the damped natural frequency and damping ratio off of that peak if you so desire. There are higher modes in the response as well, but they don't start until somewhere after 1000 Hz.

The only thing that concerns me is the low-frequency spike in the results. To me this indicates a rigid-body motion of your accelerometer, i.e. it moved at some point, such as when you smacked the soil. I don't think it is a signal-processing issue.

To summarize: I think you have processed the signal correctly and have useful response data. The low-frequency spike is not ideal, but is probably a test artifact. It's up to you whether that spike is acceptable or not.

[vu ngothanh]

Dear Mr Joe.

Your information is very very useful for me. Because I don't have a good basic knowledge about this area, I've just learned about this for a month, So, can I ask you some more question:

1. Is it necessary to make the result more smoothy? If yes, How can we make it more smooth. As you say, the region over 5000 is garbage(so we don't care about it). However, the region from 1000 to 5000 is not smooth.

2. You said about natural damped frequency and damping ration. My main purpose of this test is find them out. Do you have any idea to get them. (or any document about it, I need to learn more about it).

Thank your very much. Thank you for your time.

Wish you have a nice day.

[Joe]

1) I don't further smoothing is necessary for the quality of your results if you are only concerned about the damped natural frequency and damping of the first peak. If you'd like to try and get a smoother response curve, the best way to do that would be to take multiple test measurements and average them together in the frequency domain.

2) If your main purpose is to get the damped natural frequency and damping ratio of the soil region, you have two options.

A) Time domain: Count distance between peaks of the decaying response, the inverse of the period is the damped natural frequency. Then, use the "log decrement" method to find the damping ratio. The Wikipedia page on log decrement does a good job of explaining how to use it.

B) Frequency domain: The location of your primary response peak is the damped natural frequency. To obtain the damping ratio, use the "half power bandwidth" method, explained here. Basically, you go down by a factor of sqrt(2) from the peak, look at the width of the peak there, and use that to find the damping ratio.

If you are going to be doing this type of work more often in the future, there are dedicated texts for soil dynamics - "Principles of Soil Dynamics" seems to be well reviewed.

[vu ngothanh]

Thank you very much.

I will learn more about it through your advices.

Thanks you for your time to help me. I greatly appreciate it.

Wish you have a nice day.

Answer 2

[Image Analyst]

Just add a dot before the slash.

mobility=velocity ./ force;

I'm assuming you know how to get your data into a variable in MATLAB. If you don't know how to do that, then look into csvread(), dlmread(), importdata(), xlsread(), textscan(), fread(), fgetl(), or functions like that.

[vu ngothanh]

Dear Mr

This is my code. As you said. but I don't know about it's result, but it is quite strange with my thinking.

y=load('test_010.lvm');
a=y(:,1);
b=y(:,2);              %accelerometer
c=y(:,3);              %hammer
mobi=b./c;
plot(mobi);

This is my input data:

This is mobility result:

Could you please to help me to fix it.

Thanks you very much.

[Image Analyst]

The code looks correct, though you really should use more descriptive variable names, like I showed you. Why do you think it's not right?

[vu ngothanh]

in my project, it should be look like this:

[Image Analyst]

I can't help you. I don't know that test. Maybe you need to take the absolute value of the velocity, or maybe the what you think is the velocity is really the acceleration - I have no idea. All I know is that the syntax for dividing those two arrays, whatever they may represent, is correct. You need to dig further to find out if the arrays really represent what you think they do.

[vu ngothanh]

Yes. Thanks you very much for your time for me.

Thanks again.

Wish you have a nice day.