FFT on data acquired from an accelerometer?

Good evening!

Please, I have a collection of values acquired from an accelerometer (MPU-6050), sampled over time at a sampling period of 50 ms, evenly spaced.

I intend to do a FFT to see the frequency spectrum of that data vectors, but the results don't seem to be correct (image below).

I did not write the code in here, because of the vectors' lenght (I didn't use xls or cvs).

Here is the code: https://pastebin.com/mUtxjyKA

Thanks a lot in advance!!!

Answer 1

[Sulaymon Eshkabilov]

Hi,

First of all, why are you plotting FFT of your recorded accelerations' magnitude from x, y, z axes againt time signal? It has to be MFFT vs. frequency values.

There are some errs in your scripts: the frequency values are not generated correctly: f = fs/2*linspace(0, 1, N/2+1);

If you need the time space then it has to be w.r.t. fs = 1/.050 = 20 Hz and thus, dt = 1/fs. t = 0:dt:(1-LL)*dt, where LL is the length of your recorded signal for X_vect1 or Y_vect1 or Z_vect1.

Another unclear point is why you are performing fft of acceleration magnitudes recorded from X, Y, Z axes. It makes more sense you would need to treat acceleration collected from each axis in separate.

Here is the code for FFT - one-sided amplitude spectrum:

fs = 1/0.050; % sampling frequency
LL = length(X_vect1);
N=1024;  % Block size
f = fs/2*linspace(0, 1, N/2+1);  % Frequency values
subplot(311)
X = fft(X_vect1, N)/LL;          % Accel. along X axis
plot(f, 2*abs(X(1:N/2+1)));  grid on; 
title('One-sided Amplitude Spectrum') , shg
subplot(312)
Y = fft(Y_vect1, N)/LL;          % Accel. along Y axis    
plot(f, 2*abs(Y(1:N/2+1)));  grid on; 
subplot(313)
Z = fft(Z_vect1, N)/LL;         % Accel. along Z axis 
plot(f, 2*abs(Z(1:N/2+1)));  grid on; 
figure(2)
plot(f, 2*abs(X(1:N/2+1)), 'r', f, 2*abs(Y(1:N/2+1)), 'g' , f, 2*abs(Z(1:N/2+1)), 'b');  grid on; 
xlabel('f, frequency, [Hz]')
legend('Acc X ','Acc Y','Acc Z')
title('One-sided Amplitude Spectrum') , shg

Answer 2

[KSSV]

L = length(time_vect) ;         % signal length 
n = 2^nextpow2(L);
M_vect1_FFT = fft(M_vect1-mean(M_vect1),n);

f = f*(0:(n/2))/n;
plot(f, abs(2*M_vect1_FFT(1:n/2+1)));
grid on;
title('3 axis - frequency');
xlim([0 100])
xlabel('f (Hz)');
ylim([0 1000])
ylabel('Power');

Answer 3

[Sulaymon Eshkabilov]

A second parameter N is the block size that is needed for fft to get better resolution of your spectrum analysis. From your collected data, it seems to me you ahd better collected the data - in a bit high frequency - let's say 100 Hz. The division by LL is a signal normalization related issue.

abs and *2 are needed just becasue fft of a signal has two components, viz real and imaginary parts, and to get only its magnitude spectrum, they are needed. (1:N/2+1) takes only one half of the spectrum and the other half will be its reflection and thus, we atking only half to get one-sided spectrum.

Accelerations in X, Y, Z axes need to be treated in separate; otherwise, it does not make much sense. If you would like to compute acceleration impact level or something like that then you can choose the max of X acc, max of Y acc, max of Z acc. The selected max. acc. shall be squared, summed and taken out of sqrt as you did with your all signals. Good luck.

Answer 4

[André Rocha]

Hi Sulaymon Eshkabilov!

First of all, thank you very much!

Now, about your questions:

"First of all, why are you plotting FFT of your recorded accelerations' magnitude from x, y, z axes againt time signal? It has to be MFFT vs. frequency values."

Sorry, I don't think if I completely understand what you mean by "MFFT vs frequency values". The idea of this excercise is all about the code you gently helped me out with. I need to analyze the frequency of the vibration acquired from the accelerometer.

"There are some errs in your scripts: the frequency values are not generated correctly: f = fs/2*linspace(0, 1, N/2+1); "

Sure, thank you! I'm currently learning MATLAB and I make mistakes like that kinda often...

"If you need the time space then it has to be w.r.t. fs = 1/.050 = 20 Hz and thus, dt = 1/fs. t = 0:dt:(1-LL)*dt, where LL is the length of your recorded signal for X_vect1 or Y_vect1 or Z_vect1."

OK, but one question, that's exactly what I had written, but in a more formal, and better, way, right?

"Another unclear point is why you are performing fft of acceleration magnitudes recorded from X, Y, Z axes. It makes more sense you would need to treat acceleration collected from each axis in separate."

Yeah, the way you have done is perfect. But if I need just one vector that represents the three axes (a resulting vector), is it wrong to do calculate the magnitude like I did?

Thanks again!!!

Answer 5

[André Rocha]

And just one more thing. Like I said, I'm currently learning the basics of MATLAB. If it's not too much to ask, could you please explain with details the two lines below?

Why the second parameter of fft function, and why do you divide by the lenght of the signal?

X = fft(X_vect1, N)/LL; % Accel. along X axis

When you are about to plot the fft, why did you multiplied its magnitude by 2? Because it's single-sided? And what about the vector inside the X variable (1:N/2+1) ?? That will clear all of my questions by now.

plot(f, 2*abs(X(1:N/2+1)));

Thanks!

Answer 6

[AADIL]

hi,

this attached excel file is acceleration values and time.

i dont know how to use matlab.

can u help me do fft plot.

time domain, frequency domain and fft needed.

Answer 7

[Ivo Bastos]

@André Rocha Hi, were you able to get this running properly? I am having exactly the same problem, however, the solutions provided in this discusson did not work for me. I am using a sampling frequency of 500Hz and I am collecting data from the MEMS accelerometer using an esp32 board, via I2C. In the following code I compute two FFTs, one using the imported data from the accelerometers and the other using a generated signal that I tried to make as similar as possible to the imported data. The FFT for the generated signal seems to be computing quite nicely (needs some processing perhaps) but for the imported data it won't even compute, even thought I am using the same parameters for both FFTs. Here is the code and also the plots. I hope someone might have the answer I look for.

% Importing data from MEMS
time = readmatrix("MEMSData/data_345v2.xlsx", "Range", "H1:H2500"); % Imports Time Measurement Values 
acc = readmatrix("MEMSData/data_345v2.xlsx","Range", "P1:P2500");    % Imports Acceleration in Z axis

% DAQ parameters

Fs = 500;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = length(acc);
t = (0:L-1)*T;        % Time vector

% Creating/Importing the signal

output = sin(2*pi*10*t)+10 + sin(2*pi*48*t) + cos(2*pi*125*t);     %Signal

figure; 
subplot(4,1,1);
plot(t, output, 'r'); 
subplot(4,1,2);
plot(t, acc); xlabel("Time"); ylabel("Acceleration (m/s²)"); title("Time Domain Acceleration");
ylim([11 12.5]);

% Setting Parameters for FFT

dataLength = length(acc);
dataLengthOut = length(output);
nextPowerOfTwo = 2 ^ nextpow2(dataLength);
nextPowerOfTwoOut = 2 ^ nextpow2(dataLengthOut);
w = hanning(nextPowerOfTwo);  % Hanning window setup

plotScaleFactor = 2;
% plot is symmetric about n/2
plotRange = nextPowerOfTwo / 2; 
plotRangeOut = nextPowerOfTwoOut / 2;
plotRange = floor(plotRange / plotScaleFactor);
plotRangeOut = floor(plotRangeOut / plotScaleFactor);

yDFT_output= fft(output, nextPowerOfTwoOut);
yDFT = fft(acc, nextPowerOfTwo);
% yDFTwindow = yDFT .* w;        % Hanning window

h = yDFT(1:plotRange);
h_output = yDFT_output(1:plotRangeOut);
% hwindow = yDFTwindow(1:plotRange);  % Hanning window

abs_h = abs(h);
abs_houtput = abs(h_output);
% abs_hwindow = abs(hwindow);     % Hanning window

% Frequency range
freqRange = (0:nextPowerOfTwo-1) * (Fs / nextPowerOfTwo);
freqRangeOut = (0:nextPowerOfTwoOut-1) * (Fs / nextPowerOfTwoOut);
% Only plot up to n/2 (as other half is the mirror image)
gfreq = freqRange(1:plotRange);
gfreqOut = freqRange(1:plotRangeOut);

subplot(4,1,3);
plot(gfreq, abs_h, 'k'); xlabel('Frequency'); ylabel('Amplitude'); title('FFT');
subplot(4,1,4);
plot(gfreqOut, abs_houtput, 'r'); xlabel('Frequency'); ylabel('Amplitude'); title('FFT Output Signal');
ylim([0 1500]);