FFT questions 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??

Below is the coding of FFT. load a1.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:1998/2)).^2/1998; % magnitude of Y squared is the power Phase=angle(Y); % Calculates the phase of the FFT nyquist=1/2; % nyquist frequency freqmins=(1:1998/2)/(1998/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,';')

[Walter Roberson]

Format, please!

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[pete yang]

This is what command window says

Warning: Could not find an exact (case-sensitive) match for 'TimeSeries'. C:\Program Files\MATLAB\R2011a\toolbox\matlab\timeseries\@timeseries\timeseries.m is a case-insensitive match and will be used instead. You can improve the performance of your code by using exact name matches and we therefore recommend that you update your usage accordingly. Alternatively, you can disable this warning using warning('off','MATLAB:dispatcher:InexactCaseMatch'). This warning will become an error in future releases. > In FFTGeneral at 2 ??? Undefined function or method 'fft' for input arguments of type 'timeseries'.

Error in ==> FFTGeneral at 3 Y=fft(TimeSeries); % vector of real and imaginary components

[Walter Roberson]

Please use markup to format what you have written, as it is fairly difficult to read. The description of how to use markup is at the link I posted above.

Answer 1

[Walter Roberson]

The error message you show appears to indicate that you are attempting to take the fft() of an object of class timeseries. fft() is, however, not defined for such objects.

You could extract the data from the timeseries object using getdatasamples(). However, if your data is not sampled at exactly regular intervals, fft() applied to the data values would not produce correct results, as fft() assumes equal time intervals.