Calculating Intensity/Strength of Vibration with 3DOF I want to calculate the intensity/strength of vibration at a given location. I have measured the acceleration at this location, using an accelerometer. So my measures look for example like:
t1 - x:-3.81 y:-10.13 z:5.82
t2 - x:-5.81 y:-9.13 z:3.82
t3 - ...

These are measures along the x,y,z axis in m/s^2. My assumption is that the actual strength of the vibration at t1 can be calculated by calculating the distance in 3d space from (0,0,0). I calculate the distance using the Pythagorean theorem for 3 dimensional space. To get the average intensity I just take a lot (e.g. 1.000 .000) measurements and calculate the average of the distances in 3d space of these measures.
Is this correct/valid? Are there better ways to do this?
Would it be better to use the Jerk?
If I use the Jerk, how would I calculate it? Using the distance in 3D seems to be wrong, since the distance is always positive and is not associated with any direction. So I am not sure if it would be correct.
Edit
As I was asked in the comments about my usecase:
I want to measure the vibration intensity of a vehicle at different positions of the vehicle. In the end, I want to find the sweet spot of the vehicle, where the least vibration occurs. For example I would want to mount a vibration sensitive device there or add a driver seat to that position to reduce the vibration for the driver.
Edit2
This is a plot of approx 1400 samples. I added the 3 axis and the distance in 3d space. The distance is green, z is yellow, y is orange, x is blue.

Edit3
Ok so I just created an FFT analysis of the above data. Attached are 2 different plots. One for the 3d distance (left) and one for only the x-axis (right). Essentially the results tell me that the strongest acceleration occurs in the low frequencies, right? I should definitely get finer grained data.

BTW: I used the following function in R plot.frequency.spectrum(fft(accMeasures$x0), xlimits=c(0,1500)) the function itself can be found here.
Edit4
This is the normalized FFT Plot of the sample magnitude (left) and the  x-axis (right) of my actual samples here:

So I tested my accelerometer by hanging it to a rubber band and let it bounce. I sampled with ~100Hz. You can see two spike at 20Hz and 100Hz using the FFT analysis. When I "subtract from every value the calculated average" (I will use the word normalize for that - I hope that is correct). When I compare the non-normalized magnitude (3d distance) to the normalized magnitude, the frequency spikes change from 20 and 100Hz to 40 and 80Hz. This seems to be weird, since the spikes for every axis on its' own are at 20 and 100Hz.
This plot shows the x-axis (the axis the main movement of the rubber band occurred). On the left you can see the non-normalized FFT and on the right you can see the normalized FFT. This looks like what I would expect from normalizing the values.

This plot shows the magnitude, left is non-normalized, right is normalized. The change of frequency spikes is weird IMHO.

 A: To formalize the comments (now in chat here):
Associated jerk is probably what you want to calculate, as it is the measure of how violently something is shaken.$^1$ Jerk is the derivative of the acceleration with respect to time. To properly calculate this, you would use the formula $\left| a \right| = \sqrt{a_x^2+a_y^2 +a_z^2}$ (from the Pythagorean Theorem) to calculate the distance in 3-D. Then, using the formula $j = \sqrt{ja^2 + jb^2 + jc^2}$ you would calculate the overall jerk.$^2$
It is important to note that you probably do not need to measure short peaks as these would be absorbed by the driver's seat. Further minimizing what you need to measure, only the 0 to 80 Hz range is of importance for what you are doing. You can use mechanical damping or a low-pass filter to get just the range you want. This site might help you in analyzing your signals.$^3$
Finally, the reason for the strange change of the frequency spikes is explained as follows:

in taking the magnitude of the acceleration, you effectively take the absolute value of the acceleration. The absolute value of a sine wave gives a constant + a sine wave with the double frequency. This is why you see a peak at 0Hz and a doubling of the frequency. The alias signal shifts accordingly.$^2$

Hope this helps!
Sources:
Thanks to @lemon, @Crimson, @Previous, and @rmhleo for putting your time into the comments. This answer is the formalization of your hard work. I cited your comments by putting footnotes for those sentences that are derived from a user's particular comment. The one quote in this answer is cited in the same way. Even if certain comments are not directly cited, they helped refine the end result.
$^1$@lemon
$^2$@Crimson
$^3$@Previous
A: "Comfortable ride" is a tricky thing to quantify. Jerk is not the right metric to use. The reason it works for roller coaster design is the fact that in a roller coaster, you brace yourself against the rather large low-frequency acceleration. If you make a sharp turn to the left, you will want to lean left. If you then suddenly make a large turn to the right, you will want to turn right. The jerk of the acceleration measures how quickly you need to shift your weight to compensate.
When you are dealing with vibrations in the Hz range and above, as you are doing in this case, I don't think that jerk is the right metric - as you don't have enough time to adjust. In essence, you are driving the body at a sufficiently high frequency that the only thing that matters to your comfort is the force felt - which is proportional to the acceleration, not its derivative.
The degree to which vibrations are transmitted to the passenger will be a sensitive function of the characteristics of the chair - direction, amplitude and frequency will all affect the degree of damping, and thus the passenger's sense of "comfort".
But let's look at the physics for a moment. You were sampling a bouncing accelerometer at 100 Hz, and saw a component at 20 Hz and 80 Hz, as well as some low-frequency components.
Now when you take a Fourier Transform of a finite-width sample, you need to use a window on the data to prevent spill-over of signal into neighboring bins of the FFT. In essence, sampling a waveform for a short time means that there is considerable uncertainty in the frequency, and this shows up as a broadening of the peaks. The longer the sample, the better the frequency resolution; but even with short samples, you can improve the accuracy by applying an appropriate window (Hamming, Hanning, ...) - see the linked article above.
The next thing to keep in mind is aliasing: when you sample at a frequency of 100 Hz, frequencies above 50 Hz are meaningless; what is particularly strange here is that your frequency axis goes beyond 100 Hz, which suggests to me that something is not right; either you did not perform the FTT correctly, your sampling frequency isn't really 100 Hz, or there is a problem with the way you are plotting things.
A properly calculated FFT with 100 Hz sampling frequency should not show frequencies above 50 Hz.
To demonstrate a simple way to do this (and illustrate the problem with your plot) I wrote a short Matlab script that generates the following output.

Here is the content of the script. I recommend that you study it and ask any questions that may arise.
%FFT demo
Fs = 100;   % 100 Hz
Ts = 1./Fs; % sampling interval
L = 1024;   % number of samples

t = Ts * (0:L-1); % sampling time points

s = 2*rand(size(f)) + 12.34*sin(2*pi*20.*t); % noise with a 20 Hz signal

f1 = fft(s);

% to plot the frequency content properly we need to do a few things:
% scale the DC and other bins
% and get the frequency bins right
P1 = 2*abs(f1(1:L/2+1)); % half the power is in the bins above the middle
P1(1) = 0.5*P1(1);       % except for the DC signal

fBins = Fs *(0:L/2)/L;   % bins corresponding to each frequency

figure; 
% demonstrate the "basic" plot:
subplot(3,1,1)
bar(fBins, P1); title 'unapodized frequency plot'
xlim([0 50]);  % the bar plot by default uses wider axes - not helpful

% demonstrate the "wrong" way of taking the absolute value:
% this is where frequency doubling appears
f2 = fft(abs(s)); % taking the FFT of the absolute value of the signal
P2 = 2*abs(f2(1:L/2+1));
P2(1) = 0.5*P2(1);
subplot(3,1,2)
bar(fBins, P2); title 'unapodized - FFT of absolute value'
xlim([0 50])

% showing the use of a Hamming window to reduce spectral spill    
w = 0.64 - 0.54*cos(2*pi*(0:L-1)/(L-1)); % standard equation
f3 = fft(s.*w);                          % apply window to signal
P3 = 2*abs(f3(1:L/2+1));
P3(1) = 0.5*P3(1);

subplot(3,1,3)
bar(fBins, P3); title 'apodized frequency spectrum'
xlim([0 50])

figure; 
subplot(3,1,1)
bar(fBins, P1); title 'unapodized frequency plot'
set(gca, 'fontsize', 12);
xlim([0 50])

subplot(3,1,2)
bar(fBins, P2); title 'unapodized - FFT of absolute value'
set(gca, 'fontsize', 12);
xlim([0 50])

subplot(3,1,3)
bar(fBins, P3); title 'apodized frequency spectrum'
set(gca, 'fontsize', 12);
xlim([0 50])

As you can see, you have to apply the abs function in the FFT space - if you take the abs of the signal, you get offsets and frequency doubling. Incidentally the original signal doesn't disappear completely if you have an offset in the initial signal (as I do in my example, since the rand function gives a value between 0 and 1 and will therefore result in an average offset).
Final thought: if you add the time resolved XYZ components together using the sum of squares followed by square root, you again do a frequency doubling: $\sin^2\omega t = \frac12(1+2\sin 2\omega t)$. If you do the FFT on the individual components first, then take the sum of squares of the relative amplitudes, you avoid this problem.
