vibration analysis FFT with accelerometers

Hello,

I am doing vibration analysis with accelerometers by converting time domain to frequency domain. And fortunately, my coleage gave me a matlab fft file. However, I am a super amateur with Matlab. If possible, would you like to review and give me some tips??

Below is my questions. 1 - do I need to change N?? N is the number of time signals from accelerometers?? In my experiment,2k/second is sampling rate and did the experiment for 3 seconds. so N should be 6k?? 2 - Do I also need to modify Y??

3 - Do I need to freqhours to seconds??

%load Random_Numbers.txt; % loads two months of synthetic wind

N=length(TimeSeries) % number of data points, minutes in two years

Y=fft(TimeSeries); % vector of real and imaginary components

Y(1)=[]; % remove first component because it is the sum

power=abs(Y(1:N/2)).^2/N; % magnitude of Y squared is the power

Phase=angle(Y); % Calculates the phase of the FFT

nyquist=1/2; % nyquist frequency

freqmins=(1:N/2)/(N/2)*nyquist; %Not actually frequency in mins

SamplingFrequency=2000; % Sampling frequency in Hz

freqhours=freqmins*3600*SamplingFrequency; % frquency in cycles per hour

% periodhours=1./freqhours % peiod in hours

logfreqhours=LOG10(freqhours); % log base 10 of period in hours

scalepower=power.*freqmins'; % Scale the power to give power density

%plot(logfreqhours,scalepower), axis([-5 2 0 300]), grid on; %plot of spectrum

%ylabel('power density');

%xlabel('cycles per hour');

%title('Periodogram');

% save march06 freqhours logfreqhours power scalepower

%slogfreqhours=(logfreqhours-mean(logfreqhours))./std(logfreqhours);

%spower=(scalepower-mean(scalepower))./std(scalepower);

%p=polyfit(slogfreqhours,spower,8);

%spower8=polyval(p,slogfreqhours);

%power8=spower8.*std(scalepower)+mean(scalepower);

%figure,plot(logfreqhours,power8);

%figure,plot(logfreqhours,scalepower);

% %

%The above didn't work, probably because there are too many points

%and they are too close together at the high frequency end.

%Therefore I will group the points in intervals of 0.05 on a log10 scale.

I=1;

shortgraph=[];

ShortGraphAngle=[]; increment=0.01; %Logarithmic increement. E.g. 0.05 is a 12% increase from bin to bin

try;

    countlogfreq=2.7+increment/2;  %Log10 of 20 year's worth of hours is 5.24

    for count=1:100000;   %Far more than enough for the range of log 
frequencies

        countI=I;       %First subscript in the averaging bin

        totpower=0; %Value of spectrum before any data points are found

        TotAngle=0;

        % Calculate default value of spectrum when there are no data points in 
the bin:

        if I>1

            InterpFactor=(countlogfreq+increment/2-logfreqhours(I-1))/
(logfreqhours(I)-logfreqhours(I-1));

            totpower=InterpFactor*scalepower(I)+(1-InterpFactor)*scalepower(I-
1);

            TotAngle=InterpFactor*Phase(I)+(1-InterpFactor)*Phase(I-1);

        end

        used=0; %Flag to get averaging right whether or not there is a data 
point in the bin

        if logfreqhours(I)<countlogfreq+increment;

            used=1;

            totpower=0; %Initialise the averaging

            TotAngle=0;

            while logfreqhours(I)<countlogfreq+increment;   

                %Loop is only used when the current data point is in the 
current bin

                totpower=totpower+scalepower(I);

                TotAngle=TotAngle+Phase(I);

                I=I+1;

            end;

        end;

        shortgraph(count,1)=countlogfreq+increment/2;

        shortgraph(count,2)=totpower/(I+1-used-countI);

        ShortGraphAngle(count,1)=countlogfreq+increment/2;

        ShortGraphAngle(count,2)=TotAngle/(I+1-used-countI);

        countlogfreq=countlogfreq+increment;

    end;

catch;

    if I>countI+100;        %Put the last bits of data into the last bin but 
only if there is enough

        shortgraph(count,1)=countlogfreq+increment/2;   %data to make a 
meaningfull average

        shortgraph(count,2)=totpower/(I+1-used-countI);

        ShortGraphAngle(count,1)=countlogfreq+increment/2;

        ShortGraphAngle(count,2)=TotAngle/(I+1-used-countI);

    end

% lasterr

% shortgraph

% ShortGraphAngle

    %

% figure,plot(shortgraph(:,1),shortgraph(:,2),'r+-');

% ylabel('power density');

% xlabel('cycles per hour');

% title('Periodogram');

%

% figure,plot(ShortGraphAngle(:,1),ShortGraphAngle(:,2),'r+-');

% ylabel('Phase Angle');

% xlabel('cycles per hour');

% title('Phase Information');

end

'Mean of TimeSeries'

mean(TimeSeries)

'Standard deviation of TimeSeries'

std(TimeSeries)

'Variance of TimeSeries'

std(TimeSeries)^2

%Now integrate the spectrum to check the standard deviation

Area=0;

for J=1:size(shortgraph,1);

    %StartFreq=10^(shortgraph(J,1)-increment/2)

    %EndFreq=10^(shortgraph(J,1)+increment/2)

    %Interval=EndFreq-StartFreq;

    %Var=Var+shortgraph(J,2)*Interval/10^shortgraph(J,1)

    Area=Area+shortgraph(J,2)*increment*log(10);

end

'Calculated variance'

Var=Area*2

%stop

% NB The area needs multiplying by 2 to give the full value of the variance.

% This is because the Discrete Fourier Transform includes all frequencies from the fundamental

% frequency to the sampling frequency, N. However, the above spectrum only includes frequencies up to

% the Nyquist frequency, N/2. The Discrete Fourier transform is symmetrical about the Nyquist

% frequency so the second half of the spectrum was ignored. However, when calculating the variance

% in the random data, both halves must be included. The area of the above graph must therefore be

% doubled. See Bendat and Piersol pages 10 to 12.

%

% Output the spectrum for use in another matlab program or in excell.

% E.g. convolution of the p.d.f. of state of charge of a store

%dlmwrite('TimeSeries_FFT.dat',shortgraph,';')

[Image Analyst]

Is this any different than this: http://www.mathworks.com/matlabcentral/answers/46745-fft-questions-with-accelerometers, other than some of the lines are now formatted correctly (something you could have done by editing your prior post)?