MATLAB: Calculate the frequency response of acceleration data

hi to all,

I have a vector in Matlab which contains samples acquired at 10Hz from an accelerometer mounted on a frame's vehicle.
This vector contains the acceleration values along the Z-axis.
The acquisition time is 103 seconds.

I would like to plot the frequency values in Hz since I'm trying to study the terrain frequency response.

I tried to use the fft function available in Matlab, but I do not know how to proceed.
Can you help me, please?

"The solutions and answers provided on Experts Exchange have been extremely helpful to me over the last few years. I wear a lot of hats - Developer, Database Administrator, Help Desk, etc., so I know a lot of things but not a lot about one thing. Experts Exchange gives me answers from people who do know a lot about one thing, in a easy to use platform." -Todd S.

Did you witness this test? If so, how would you describe the vehicle motion?
Have you plotted the raw acceleration data? How well does this relate to the observed motion?
Do you understand how to interpret the output of the an FFT?

Rocco GalatiR&D EngineerAuthor Commented:

Thank you for your help, d-glitch.

I plotted the raw data and this is the result:

The vehicle was moving in the mud, so I expect to see "soft and smooth" values compared to the asphalt case.
The peak at the end of the graph is caused by a terrain change (from mud to grass).
I studied the FFT when I was at the university even if I do not remember it very well (i studied it in 2009).

I would like to find the Fundamental frequency in order to check if I'm able to compare the frequency in different tests.
If this is true, it would be possible to use the fundamental frequency as a criteria to detect the terrain type.

Azure has a changed a lot since it was originally introduce by adding new services and features. Do you know everything you need to about Azure? This course will teach you about the Azure App Service, monitoring and application insights, DevOps, and Team Services.

You were accelerating in the mud for 7.5 seconds. The signal is noisy but almost always positive.

You transition from mud to asphalt between 7.5 and 8.5 seconds. The signal peaks, since you suddenly have more power than you need.

You drive on asphalt between 8.5 and 10.2 seconds. Minimal and smooth acceleration.

I think you can probably interpret the raw data better than the FFT.
What sort of "fundamental frequency" were you expecting to see?

What sort of time resolution do you need for the terrain classification?
In this case it looks like the terrain transition takes 1 second. How fast were you moving?

10 seconds is clearly the wrong size chunk of data to analyze.
The sample size should probably cover no more that the time required to travel one wheel base distance.

Rocco GalatiR&D EngineerAuthor Commented:

I'm working at 10Hz, so I expect to have 10 samples per second and this should be enough for terrain classification.
The raw acceleration data comes from the Z-axis, so the graph of the raw data shows the interaction between the terrain and the robot.

I expect to see high values when the vehicle travels on gravel and lower values when it moves on sandy terrains.
The vehicle moves about 0.8 - 1 m/s and it uses tracks and not wheels so in this case I expect to see lower values along the Z axis since the tracks offer more stability than the wheels.

I would like to check if the fundamental frequency changes significantly depending on the terrain.

Rocco GalatiR&D EngineerAuthor Commented:

Take a look at these data coming from a test over the asphalt.
The first picture is related to the raw data while the second comes from the FFT.

I think there is no question that this will work, but you have some work to do.
The closest thing you have to a fundamental frequency is track segments per second.
The right sample size for the FFT is still going to be the time required to travel one vehicle length.

If you run FFT's for the the first 128 sample you will find a MUD signature.
If you run FFT's for the the last 128 sample you will find a ASPH signature.
You may have to do tests for gravel, sand, snow, ...

And you may be able to find signatures for the transitions: half the track in mud, and half on gravel.

Rocco GalatiR&D EngineerAuthor Commented:

I already did all the test on different terrains/surfaces and now I would like to find a way to compare them by using some kind of signature.

Honestly, I already successfully considered the currents, the torque, the pitch and roll angle, the slip effects and with these parameters I'm already able to detect the terrain, but I would like to add also the information about the frequency.

Please, can you explain me what you mean with:

If you run FFT's for the the first 128 sample you will find a MUD signature.
If you run FFT's for the the last 128 sample you will find a ASPH signature.

Why do you suggest to consider the first 128 samples?

The vehicle takes about 2.5 seconds to cover the track length (the track contact surface is about 100 cm)

Your last two plots are examples of way to much data.
That little peak around 750 probably corresponds to track_segments/second. Everything else is noise.

Do you know why you don't see much in the way of negative signals? You may have an offset issue.

The secret to interpreting an FFT is to look at the results of perfect data. For example: Assume 16 samples at 10Hz.
The fastest signal you can see is:
[ +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 ] => 5.000 Hz
The slowest signal you can see is:
[ +1 +1 +1 +1 +1 +1 +1 +1 -1 -1 -1 -1 -1 -1 -1 -1 ] =>0.625 Hz

Rocco GalatiR&D EngineerAuthor Commented:

in effect, I do not know if it is a coincidence or not, but the same peak is present again at 1500, exactly after 750.

What do you mean when you say "track_segments/second"?
Are you referring to the numbers of track links/pad which touch the terrain in a second?

It is not a coincidence. The FFT is symmetric, with the lower half mirroring the upper half.

Rocco GalatiR&D EngineerAuthor Commented:

How the FFT window is connected to the time and the number of samples?
For example, why the FFT goes from 0 to 1000 and then it mirrors from 1000 to 2000?

Looking at your 1024 sample Mud/Asphalt case: You can hope to find a signature here because the FFT combines both cases and smears everything out.

Your 2200 sample case also seems to have a terrain change at the beginning. How many samples were there in this case?
I am not sure what Matlab does with odd sample sizes (not 1024 or 2048) so I don't know exactly what frequency is associate with the 750 spike.

If the sample size was 2048 => 204.8 s then the lowest freq would be 1/204.8 s = 4.883 mHz
That would put the peak at 750 x 4.883 mHz => 3.66 Hz

Are you running Matlab on the robot? I'm sure running smaller samples will save you processing time and give you less confusing results.

Lock again my 16-sample max and min freq case. There are 16 samples, but there are only 8 possible frequency values in the result.

What the 16-sample FFT does is take 16 amplitude values and return 16 complex numbers, corresponding to 8 frequency and 8 phase values.
This matrix is always symmetric, since you can reverse the time signal and not change the FFT.

When you take the absolute value, you basically discard the phase information.

Rocco GalatiR&D EngineerAuthor Commented:

The total time for the second case is 219 seconds and the elements in the raw acceleration vector are 2182.

I don't know why it seems that there is a terrain change since the only difference is that in some area the asphalt was wet.
I'm not running Matlab on the robot since I analyze the data offline on a different computer.

So, what I should do is: running short tests and check for the signature by finding the lowest and the highest frequency for each terrain and then use these two "parameters" to detect the terrain features?

In the 219 second case, some thing clearly happens around 20 seconds. I don't know what it is.

If the robot is 1m long and travels 2.5 m/s, then you should probably look at clean data files of 64, 128, and 256 samples for each of your terains. By clean, I mean the vehicle should spend all the time on the same terrain.
If you can look at the raw data and see any artifacts, then the data is not clean.

Once you have the FFT's of these data files you should be able to look at them, and pick out a few artifacts.
They won't be at the max or min frequencies, they will be in the middle.
And they will be at the same real frequency, no matter what the sample size is.
These few artifacts are the signature you are looking for.

If you rerun clean 64, 128, and 256 samples from your 219 second file (skip the first 20s), you should see the same peak around 3.66 Hz in every case.

Rocco GalatiR&D EngineerAuthor Commented:

Thank you a lot for all your support! I really appreciate it!

If the robot is 1m long and travels 2.5 m/s, then you should probably look at clean data files of 64, 128, and 256 samples for each of your terains. By clean, I mean the vehicle should spend all the time on the same terrain.

I'm sorry if this is a stupid question, but how did you calculate the 64, 128 and 256 values?
I would like to better understand your considerations since they are very interesting!

However, there is another aspect I'm not able to understand very well: when the robot moves, the number of track pads which are in contact with the terrain is always the same and also the track's fingerprint is always the same size.

These are the tracks I'm using for my robot (in total, I use two tracks, one for each side and the robot is a skid-steer)
The second picture shows the links for each tracks and if, for example, the robot moves straight forward, the number of links on the terrain is always the same.
So, why do you suggest to consider the "track_segments/second"? shouldn't this value be always the same?

Your robot moves 2.5 m/s ==> 5.6 m/hr which is pretty fast.
You have 1 meter of track on the ground and probably 2.5 meters of track total.

Is the track a straight belt with a coupling link or a continuous loop?
And how many links in total? And how many teeth on the drive and guide wheels? I will assume 30 and 12 for now.

What frequencies do you expect to see?
If you have one bad link you will see a bump every 2.5 s or 0.4 Hz
If you have one bad tooth on the drive wheel you will see a bump every 1.0 s or 1.0 Hz

The links hit the ground at a rate of 30 every 2.5 s which is 12 Hz.
To see this component clearly you would have to sample faster, at least 24 /s and probably 40 /s.

You should also look at the RPM of you motor and gears.
These will contribute even higher frequencies which you can't resolve either. But there will be noise.
It would be good to look at the output of the accelerometer and do an FFT at 10 kHz.
If you have a lot of noise at freq greater than 5 Hz you need to filter the signal before running the FFT. You can do this in Matlab as well.

Suppose you have a bad tooth on a drive wheel. You will see a spike on the FFT at 1 Hz on the FFT when you are running at 2.5 m/s. But this is speed dependent. The spike location will move with speed. In fact the whole signature will move with speed.

Do you expect to see more vibration at harder terrain from links hitting the ground? If so this would be my prescription:
Set your sample size to 256. This will give 20 Hz max frequency and 20/128 => 0.16 Hz resolution.
Once you find a signature, you can change the filter size to speed up processing or increase resolution.

Rocco GalatiR&D EngineerAuthor Commented:

Thank you again for your time.

These are the real tracks I'm using.
The vehicle can move with a velocity up to 1 m/s, its weight is about 300 Kg.
It uses two tracks with 36 links and the main sprocket has 15 teeth, the weighs for each track is 34 Kg.
The tracks use a continuous chain and they are made in 100% natural rubber.
The overall length the each rubber track is 254 cm and each lug is 3 centimeters height.
The track has a width of 33 cm and the contact ground patch has a length of 100 cm.
The pitch is about 7 cm and so there are 100/7 = 14 links on the ground.

The ratio of the gear box is 1:60 and the maximum velocity for the sprocket is 330RPM (the maximum RPMs for the motors are 2200).

The tracks are pretty new so they do not have any bad link or bad tooth, but each one moves in a very specific way for various reasons (different tension on the belts, different velocity, different rubber response to stress,..).

I expect to have high vibrations on harder terrains: on the sand, for example, I expect less vibrations than on the asphalt.

I acquire the accelerations at 160 Hz, so I have 160 samples for each second.

Here you can find the velocity vectors (vs, vd), the vibration vector and the raw Z acceleration vector sampled at 160 Hz.

Unfortunately, I still do not understand how you are able to calculate the values and the frequencies.

For example, with my data, if I have a bad link, I should see a peak every 1.5 seconds or at 0.66Hz.
Then, the links hit the ground in a group of 13 each 1.5 seconds or at 19.5Hz.
Is it correct?

EDIT:
This is what I get when I plot the fft of the raw acceleration data without any filtering and at 160 Hz (asphalt case - the second one)

>> If you run FFT's for the the first 128 sample you will find a MUD signature.
>> If you run FFT's for the the last 128 sample you will find a ASPH signature.

As you search online about frequency scaling , what may be confusing is that sometimes you'll see +0.5 or +pi. These are just two-sided values of the same thing; namely, the highest frequency content - call it Fmax (Hz). Fmax = Fs/2 where Fs is the sampling rate (samples/second). So, .5 or pi are just scaled values that represent Fmax (Hz).

Wildly busy yesterday, but I am back. phoffric seems to be much better at Matlab than I am. Useful info.

Your last plot is excellent. Probably all the information you need.
It shows a strong fundamental frequency and the odd harmonics. Classic.

35k samples at 160 samples per second [sps] ==> 218.8 seconds
The fastest freq you can see is (160Hz)/2 ==> 80 Hz
The slowest freq you can see is 1/(218.8s) ==> 4.571 mHz <== This is also the resolution of the x-axis.

Your plot shows the fundamental (1st harmonic) and the 3rd, 5th, and 7th. And then the mirror image.

The 7th harmonic is at approx 16,000 ==> 16k * 4.571 mHz ==> 73.1 Hz
And the fundamental is at 10.45 Hz

The sort of information you might consider for a signature is the amplitude of the 1st harmonic and the ratio of the 1st and 3rd. And let's say you are satisfied with a resolution of 0.1 Hz.
The minimum sample size required for this resolution is T = 1/(0.1 Hz) = 10 s ==> 1600 samples.
You should see a similar plot, with the same pattern if you run this smaller sample.

There are several opportunities for optimization:

You might discover that you can get by with 0.2 or even 0.5 Hz resolution. That would let you use 5 or 2 second samples for faster decisions.

If you don't need all the harmonics, you can filter and decimate the data to reduce your FFT computation requirements.

If 0.2 Hz resolution and two harmonics are enough:
Sample time ==> 5s ==> 800 samples Filtering and decimation let you run a 400 sample FFT

If all you care about is the amplitude of the fundamental and 0.5 Hz resolution is enough:
Sample time ==> 2s ==> 320 samples Filter at 20Hz and decimate by 5 and for a 64 sample FFT.
Even a PIC or Arduino can handle this pretty easily.

It is crucial to filter before you decimate the data. You can wind up with severe aliasing problems if you forget.

Rocco GalatiR&D EngineerAuthor Commented:

Hello guys, I'm sorry for my delay, but I spent several hours during this afternoon to try to understand all your suggestions.

First of all, thank you for all your assistance on this topic.

I set the frequency on the x-axis (as phoffric suggested) and I changed my Matlab code in this way:

and I run it over different dataset and these are the results.
i changed the time windows according to the condition of each test, later, I will define a standard time window which will be always the same for all the tests.

Case 1: SAND

I get the first harmonic with a maximum value of 3102, the third at 1077, the fifth at 711 and the last at 588.

Case 2: GRAVEL

In this case, the data are: #1: 2087, #3: 683, #5:481, #69.85

Case 3: ASPHALT

In this case, the data are: #1: 6707, #3: 2394, #5:1665, #1433

Case 4: HEAVY MUD

This time, the data are: #1: 279, #3: 103, #5:70, #59

These values change significantly in each test and so I think I can use them as some kind of fingerprint for each terrain.
May be, I can consider only the first harmonic and then the ratio of the #3 and 5# harmonic as d-glitch suggested few posts ago.

Do you think these considerations can be correct?

Is there a way to accurately define the frequency of each harmonic?
At the moment, I can only approximately estimate their frequency, but how can I understand their correct value?
Should I consider the peak of each harmonic as its "position" on the frequency axis?

The amplitude values are still m/s^2 ?

P.S. Since I also have the pitch, yaw and roll rate and X,Y,Z accelerations I was thinking to include all of them in these considerations.

>> Should I consider the peak of each harmonic as its "position" on the frequency axis?
Depends how precise you need to be. Looking at the main frequency, there are two frequencies near the true peak. That is, the energy is divided almost equally in two adjacent frequency bins. So the true peak is likely at a frequency in between those two adjacent bins. You can append (i.e., pad) your time domain with lots of zeros (i.e., the number of samples appears to increase with the end being all zeros). Then your fft should have finer grained frequency interpolated resolution.

Rocco GalatiR&D EngineerAuthor Commented:

You can append (i.e., pad) your time domain with lots of zeros (i.e., the number of samples appears to increase with the end being all zeros). Then your fft should have finer grained frequency interpolated resolution.

You can interpolate the FFT by zero padding. Zero padding enables you to obtain more accurate amplitude estimates of resolvable signal components. On the other hand, zero padding does not improve the spectral (frequency) resolution of the FFT. The resolution is determined by the number of samples and the sample rate.

Here is some of your revised code showing how to zero pad your data.
Data file created using posted data in https:#a42430623
I also modified the raw data so that your plots in https:#a42430458 would look better.

raw_accelearation = csvread('data_42430623.txt');% raw_accel = dlmread('data_42430623.txt');N1 = length(raw_accelearation)% remove the data before the motionraw_accel2 = raw_accelearation(210:N1);N = length(raw_accel2)% remove the DC component in the frequency spectrum % to improve the scaling in the plotsraw_accel = raw_accel2 - mean(raw_accel2);figureplot(raw_accel);xlabel('Sample Number');ylabel('Raw Acceleration');Fs=160;Y = fftshift( abs( fft(raw_accel) ) );fshift = (-N/2:N/2-1)*(Fs/N); % zero-centered frequency rangefigureplot(fshift,Y)grid on;xlabel( ['Hz (' num2str(N) ' bins)'] );ylabel('Freq Energy');Y = fftshift( abs( fft(raw_accel, 32768) ) );N3 = length(Y)fshift_pad = (-N3/2:N3/2-1)*(Fs/N3); % zero-centered frequency rangefigureplot(fshift_pad,Y);grid on;ylabel('Freq Energy');xlabel( ['Hz (' num2str(N3) ' bins)'] );% Here is how to zero pad when using only 1 arg to fft:zeroPad = zeros(32768-N, 1);accel_pad = [ raw_accel; zeroPad ];N2 = length(accel_pad)Y_pad = fftshift( abs( fft(accel_pad) ) );fshift_pad = (-N2/2:N2/2-1)*(Fs/N2); % zero-centered frequency rangefigureplot(fshift_pad,Y_pad);grid on;ylabel('Freq Energy');xlabel( ['Hz (' num2str(N2) ' bins)'] );

Here are zoomed in sections of the fft using different number of bins. As the number of bins increase, you see a smoother interpolation of the discrete-time Fourier transform (DTFT) because more points are sampled from the DTFT.

Experts Exchange Solution brought to you by

Your issues matter to us.

Facing a tech roadblock? Get the help and guidance you need from experienced professionals who care. Ask your question anytime, anywhere, with no hassle.

Why do you subtract the mean values from the same vector?

The vector posted in https:#a42430623 was a downsampled version of the original vector acquired at 160Hz.
I uploaded the original vector - here - and it contains 35103 elements.

Unfortunately, I am unable to load your .mat files. I subtracted the mean from the signal because the fft that you plotted in post ID: 42430458 had a very large zero frequency energy that swamped out other data that you are interested in viewing.

Suppose that you had a very tiny signal on a DC voltage line. The mean of that signal is just going to be the DC value. By subtracting out the mean of the signal, all you are left with is the essential aspects of the signal. When taking the fft, the signal will no longer be swamped out by the DC energy level.
it may be that the mean and variance of your signal are also parameters that you can use to classify the terrain.

You are in fact applying the window called a rectangular window . But it's spectrum is a sync function which is noted in the above link .

It's more than this solution.Get answers and train to solve all your tech problems - anytime, anywhere.Try it for freeEdge Out The Competitionfor your dream job with proven skills and certifications.Get started todayStand Outas the employee with proven skills.Start learning today for freeMove Your Career Forwardwith certification training in the latest technologies.Start your trial today

X = [ Your data ... ]Y = fft( X)This the complex FFT of your dataplot( abs( Y))You take the absolute value to get amplitudeThis is the plot: