Vibration Analysis by using Fast Fourier Transform

Hi! I have approx 65k data taken by accelerometer within 10 minutes and i need to transform it in frequency domain by FFT and also had Ax,Ay,Az,Gx,Gy,Gz columns. that's confusing can I implement fft to every single column and add them over each other while plotting with "hold on"? I mean plotting the transform for column of Ax then hold on and plotting and implementing same code for column of Ay then hold on and plotting etc. Can you help me about how I should start? Thanks in advance

Answer 1

[dpb]

>> help fft
fft Discrete Fourier transform.
  fft(X) is the discrete Fourier transform (DFT) of vector X.  For
  matrices, the fft operation is applied to each column. ...

You don't need to do anything but put all the components into a single array

A=[Ax,Ay,Az,Gx,Gy,Gz];       % presuming they're all column vectors
F=fft(A);                    % fft each by column

See

doc fft  % for example of how to compute the PSD and scale frequency, magnitude.

[dpb]

• Can you tell me about which parameters are you talking about?

Do you mean that when you call prepareSurfaceData, the X, Y, & IntensityOfOneRow are different than XOut, YOut & ZOut respectively?

• Or, do you mean that every time you call the fit function, you get varying coefficients?

You can also contact Technical Support at MathWorks, if you are looking for a quicker resolution.

[Enez Furkan Cihan]

Thanks for your care, sir when I generate a column vector by creating a time vector with 'linspace' as independent my acce. data, the code works well and getting figure. But, when use my acce data column(Ax), it doesn't work. Isn't it [Ax] already time-dependent? I miss something. Thx again!

[dpb]

Don't understand the question precisely, sorry...

The acceleration will have been collected at whatever sampling frequency it was and stored in (presumably a column) vector or array with each column holding a particular trace (x, y, z of a triax accelerometer, for example) and will thus be N rows by M columns.

The time associated with those points will be 0:dt:(N-1)*dt or linspace(0,T,N) to match the length of the acceleration vector(s) in time.

If the data are sampled at a fixed dt (as assumed by FFT analysis), then there's no difference excepting a scale factor of plotting the time trace versus the time vector or the ordinal position in the array--it's simply a pure multiplication by dt of the vector 0:N-1.

When you take the transform, however, that transforms the time domain to frequency and the frequency is based on the sampling rate and length of the sample as df=1/T where T is the overall time of the sampled record. That's what the linspace() line in the example from the FFT documentation does is to compute that associated frequency.

The Matlab FFT returns the two-sided transform from -Fmax:+Fmax with the DC component in the middle; the N/2+1 subscripting simply returns the positive frequency half and the scaling shown accounts for the length multiplier in the magnitude. If you work through the example given there for a specific magnitude and frequency you'll see you get the peak displayed at the proper frequency and its height of the correct magnitude within the discretization effects inherent in finite sampling.

What seems to be the actual issue -- post any error message in context or attach the figure with any concern the above doesn't answer.

[Enez Furkan Cihan]

I'm so thankful for this article, I have thought that I can explain the problem w/ images.

Fs=1000;
t = linspace(0,1,1000); % N=1000, T=1
k = cos(2*pi*250*t)+0.5*randn(size(t))
kdft=fft(k);
kdft=kdft(1:length(k)/2+1); 
freq = linspace(0,Fs/2,length(k)/2+1);
plot(freq,abs(kdft));
xlabel('Hz');ylabel('Magnitude');
figure;
plot(freq,20*log10(abs(kdft)));
xlabel('Hz');ylabel('dB');
title('Plot in dB');

here I generate a random column vector and getting

and

but; when write this;

Fs=97.885;
Axdft=fft(Ax);
Axdft=Axdft(1:length(Ax)/2+1,:); 
freq = linspace(0,Fs/2,length(Ax)/2+1);
plot(freq,abs(Axdft));
xlabel('Hz');ylabel('Magnitude');
figure;
plot(freq,20*log10(abs(Axdft)));
xlabel('Hz');ylabel('dB');
title('Plot in dB');

getting nothing for the figures.
By the way, Ax has 64998x1 value from acce's data and overall time is 664 seconds. Thank you once again.

[dpb]

That figure doesn't match the code; note that your peak frequency is at 25 Hz whereas the signal was generated at 250 Hz. Your plot Fmax --> 50 Hz instead of Fsamp/2--> 500 Hz; ergo you've got a factor of 10 in the length calculation somewhere in the frequency vector you used vis a vis that it should be.

Also, to get the proper magnitude of the PSD peak, you need to scale the FFT output by 1/L where L is the length of the signal. Notice the figure in the demo is a sum of two frequencies with different magnitudes (0.7 and 1.0) and the output PSD peak magnitudes are near those values. Your multiplier is 1 but the magnitude is 500 in round numbers, not 1.

Also note in the example the FFT of length=NFFT is scaled by L, the time series length. This is owing to the fact it used nextpow2() to expand to the next power of two but those additional values are all zeros so there's only power in the L original values.

I again commend that example to you to follow VERY carefully and keep reading over and over until you fully understand what they're doing and why.

As for the other data and its plot, since you didn't show it can't really judge but I'd venture it's owing to scaling issues similar to the above -- quite possibly you've got an axis limits (could be either x- or y-axis) that simply don't encompass where the frequency/output PSD values are.

First, however, ensure that your computed f-->0:Fs/2, the number of elements matches N/2 and the max(abs(Y)) ~ the input energy magnitude at the dominant frequency.

ADDENDUM

Actually, the frequency axis looks like it would be correct for the "real" data. If your Fsamp is just under 100 Hz, then a Fmax < 50 Hz is correct and that looks to be where the plots end. Maybe you mixed up between having compute that sampling frequency vector and the other that would go with the 1000 Hz sampling rate and then the first 50 Hz on a 500 Hz axis range would be all scrunched up on the left and a scaling of 1/65000 points compared to 1000 would make the signal amplitude 1/65th of the other on that basis alone so if you plotted onto another existing axis the "real" data are just so small in magnitude in comparison they just aren't of a scale to show up.

[Enez Furkan Cihan]

I have carefully read and sure that 0:fs/2:L/2+1 matches N/2, but I didn't get max(abs(Y)) where you mean actually max(abs(Ax)) for my code and it returns 'NaN'(not a number) so far so good here? As you said, there's a fact that axis limits(could be either x or y axis) simply dont encompass. I've tried to fix it with axis([0 Fs/2 0 1000]) where appointed y=1000 but whatever I appoint, it doesn't change anything. BTW I haven't tried to plot onto another column yet, just trying on column of Ax (64998x1) . Isn't 'Axis' command enough to fix axis limits? Here I guess it's not enough. Right?

Thank you once again.

[dpb]

Well, NaN isn't good; that means something blew up somewhere...w/o the data I don't know there's much more I can tell you other than I'd replicate the example code precisely with no more than a variable change. I note for some reason you didn't go ahead and compute the PSD but only in the plot statement; why? Makes it harder to debug when you don't have the variables to look at.

PLOT() without any limits set will auto-range to the ranges of the variables in the x-, y- locations. But, if they're all NaN, NaN's aren't plotted which would explain there not being anything visible.

What does

sum(isnan(abs(Axdft)))

return? It should be zero, but if it is same as length(Axdft) then that means the whole array is NaN and something really went south. I've got to run, see if

fft([NaN rand(1,127)])

returns all NaN, maybe. If the one NaN propogates through, maybe that's the problem in the original data series?

ADDENDUM

Back for a sec...

>> F=fft([nan; rand(127,1)]);
>> all(isnan(F))
ans =
   1
>>

So, indeed, if there were a single NaN in your input time series, the whole FFT will be NaN and nothing would show up...but, if you don't rescale the axes, after a plot the axes limits would be 0-1 for both x- and y- axes as there wasn't any data to auto-scale from. But, if you set xlim/ylim manually, they would have those limits.

[dpb]

Attach an m-file with the data for the particular sensor you're having so much trouble with...I'll try to find time to take a look at it and see what gives...

[Enez Furkan Cihan]

Thank you so much for everything! I have succeeded and got figure as dB and I have already wanted it. In excel file, among approx 65k data, there were some NaNs and slipped lines. Actually, never thought I would examine 65k data but I did. That's why it didn't work. As a result, I got that whatever it's long or short data, always control it. Special thanks to you for your incredible attention.

Answer 2

[Bjorn Gustavsson]

Yes, go ahead and fft individual columns, for plotting it seems most likely you would plot the absolute of the fft, or possibly the log of the absolute:

ftAx = fft(Ax);
plot(log10(abs(ftAx)),'b')
hold on
ftAy = fft(Ay);
plot(log10(Ax),'r')

You should also make sure that you have a regular sampling intervals.

Finally you might be interested in the fftshift function and mathworks documentation about what frequency-scale you get.

HTH

Answer 3

[Robert Hughes]

I have seen data like yours. Ax Ay Az are accelerations in 3 orthogonal axes. Gx Gy Gz are gyro outputs. They must be handled differently to be useful. There may be a "transformation" needed to get data to represent actual physical motion. The data is most likely is in "counts" and would need to be converted into useful engineering units. You can do an fft of each column and that's what you have--6 ffts in counts.

[dpb]

While true perhaps he's operating on A/D counts, that really only makes any difference in scaling to engineering/absolute units; the FFT would be same to a scale factor so that isn't the root problem.