%**************************************************************************************************************************************************************
%   本程序為主程序,目的是運用統計分析法計算WM1000發動機故障診斷所布各個測點加速度響應的自協方差函數
%   
%   中國北方發動機研究所試驗測試中心  蘇志霄   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';
%**************************************************************************************************************************************************************
%故障診斷統計分析法：計算各個測點加速度信號數據的自協方差函數
for i=1:m;
    figure(i);
    R=zixiefangcha(acc(i,:));
    length=round(2/3*max(size(R)));
    subplot(2,1,1);
    plot(t(1:length),R(1:length));set(gca,'FontSize',7);
    grid;
    xlabel('time/Sec');
    ylabel('Self-Convariance Coefficient');
    title(['The Self-Convariance Curve and its Frequency Spectrum of ',num2str(i),' point']);
    subplot(2,1,2);
    f=5000*(0:511)/1024;
    A=fft(R,1024);
    pu=sqrt(A.*conj(A))/1024;
    plot(f,pu(1:512));xlabel('Frequcncy/Hz');set(gca,'FontSize',7);
    hold on;
    [a n]=max(pu);
    plot(f(n),pu(n),'r.','MarkerSize',15,'EraseMode','none');
    N=num2str(n/512*2500,'%.1f');
    text(f(n)+30,pu(n),[N,'Hz'],'FontSize',8);
    hold off;
end;