%**************************************************************************************************************************************************************
%   本程序為主程序,目的是對WM1000發動機故障診斷所布各個測點的加速度響應進行FFT，并繪制時頻域響應圖和功率譜密度函數圖
%   
%   中國北方發動機研究所試驗測試中心  蘇志霄   2006年10月
%**************************************************************************************************************************************************************
m=13;                                                             %m取值為測點數目
%讀入時間加速度數據,各測點在每個測試時間內時間數據列是一致的，因此統一標識
t=acceleration_X;
len=max(size(t));
acc=zeros(m,len);
acc(1,:)=acceleration_Y';acc(2,:)=acceleration2_Y';acc(3,:)=acceleration3_Y';acc(4,:)=acceleration4_Y';acc(5,:)=acceleration5_Y';acc(6,:)=acceleration6_Y';
acc(7,:)=NONE_Y';acc(8,:)=NONE2_Y';acc(9,:)=NONE3_Y';acc(10,:)=NONE4_Y';acc(11,:)=NONE5_Y';acc(12,:)=NONE6_Y';acc(13,:)=NONE7_Y';
%**************************************************************************************************************************************************************
%將各測點加速度信號進行Fourier變換
f=5000*(0:1023)/2048;
A=zeros(m,2048);
pu=zeros(m,2048);
for i=1:m;
    A(i,:)=fft(acc(i,:),2048);
    pu(i,:)=sqrt(A(i,:).*conj(A(i,:)))/2048;
end;
%**************************************************************************************************************************************************************
%繪制各個測點的時域頻域信號，并標出時域最大加速度值和頻域最大加速度對應的頻率值
for i=1:m;
    figure(i);
    subplot(3,1,1);
    plot(t,acc(i,:));set(gca,'FontSize',7);grid;xlabel('time/Sec');ylabel('acceleration/g');
    title(['The time,frequcncy response and power spectral density of  ',num2str(i), ' point, maxacc=',num2str(max(abs(acc(i,:))),'%.1f'),'g']);
    hold on;
    [a n]=max(abs(acc(i,:)));
    plot(t(n),acc(i,n),'r.','MarkerSize',15,'EraseMode','none');
    hold off;
    subplot(3,1,2);
    plot(f,pu(i,1:1024));set(gca,'FontSize',7);xlabel('Frequcncy/Hz');ylabel('acceleration/g');
    hold on;
    [a n]=max(pu(i,:));
    plot(f(n),pu(i,n),'r.','MarkerSize',15,'EraseMode','none');
    N=num2str(n/1024*2500,'%.1f');
    text(f(n)+30,pu(i,n),[N,'Hz'],'FontSize',8);
    hold off;
    subplot(3,1,3);
    periodogram(acc(i,:),[],'onesided',512,5000);set(gca,'FontSize',7);xlabel('Frequcncy/Hz');ylabel('PSD/dB');title([]);
    
end;