FFT on data acquired from an accelerometer?

107 次查看（过去 30 天）
André Rocha
André Rocha 2019-5-6
回答： Ivo Bastos 2023-3-21
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).
ensaio 2.PNG
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!!!

请先登录，再回答此问题。

采纳的回答

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

更多回答（6 个）

KSSV
KSSV 2019-5-6
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');
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.
André Rocha
André Rocha 2019-5-6
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!!!
André Rocha
André Rocha 2019-5-6
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!
AADIL
AADIL 2023-2-8
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.
Ivo Bastos
Ivo Bastos 2023-3-21
@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]);

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

产品


版本

R2015b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by