%**************************************************************************************************************************************************************
%    該程序所要達到的目的是對發動機振動加速度信號進行二進制小波分解，繪制各層信號的Fourier頻譜圖和各子空間能量的棒圖,
%    通過各層細節信號來識別判斷發動機的故障特征
%
%    中國北方發動機研究所試驗測試中心  蘇志霄   2006年10月
%**************************************************************************************************************************************************************
m=13;                                                             % m取值為測點數目
n=6;                                                              % n為小波分解的尺度
fsamp=5000;                                                       % fsame為采樣頻率
%讀入時間加速度數據,各測點在每個測試時間內時間數據列是一致的，因此統一標識
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';
%**************************************************************************************************************************************************************
%對分解結構[c,l]中的各層高頻部分進行重構,并繪制各層信號及其Fourier頻譜圖，為了清晰可見，僅繪制前2秒的時域信號
f=fsamp*(0:1023)/2048;
length=round(2*fsamp);
for i=1:m;
    R=xiaobofenxi(acc(i,:),n);
    figure(i);
    subplot(n+1,2,1);
    plot(t(1:length),R(1,1:length));
    set(gca,'XTickLabel',{''});
    set(gca,'YTickLabel',{''});
    Ylabel(['a',num2str(n)],'FontSize',7);
    set(gca,'FontSize',7);
    title(['Wavelet signal in ',num2str(n),'-scales,' num2str(i) '-point'],'Fontsize',8);
    set(gca,'FontSize',7);
    A=fft(R(1,:),2048);
    pu=sqrt(A.*conj(A))/2048;
    subplot(n+1,2,2);
    plot(f,pu(1:1024));
    Ylabel(['FFT(a',num2str(n),')'],'FontSize',7);
    set(gca,'XTickLabel',{''});
    set(gca,'YTickLabel',{''});
    set(gca,'FontSize',5);
    hold on;
    [a n1]=max(pu);
    plot(f(n1),pu(n1),'r.','MarkerSize',7,'EraseMode','none');
    title(['Freq. spectrum of wavelet in ' num2str(n) '-scales,' num2str(i) '-point'],'Fontsize',8);
    N=num2str(n1/1024*2500,'%.1f');
    text(f(n1)+40,pu(n1),[N 'Hz'],'FontSize',7,'Color',[1 0 0]);
    hold off;
    for j=1:n;
        subplot(n+1,2,2*j+1);
        plot(t(1:length),R(j+1,1:length));
        set(gca,'YTickLabel',{''});
        if (j~=n);
            set(gca,'XTickLabel',{''});
        end;
        if (j==n);
            xlabel('time/Sec','FontSize',7);
        end;
        Ylabel(['d',num2str(n+1-j)],'FontSize',7);
        set(gca,'FontSize',5);
        A=fft(R(j+1,:),2048);
        pu=sqrt(A.*conj(A))/2048;
        subplot(n+1,2,2*j+2);
        plot(f,pu(1:1024));
        Ylabel(['FFT(d',num2str(n+1-j),')'],'FontSize',7);
        set(gca,'YTickLabel',{''});
        if (j~=n);
            set(gca,'XTickLabel',{''});
        end;
        if (j==n);
            xlabel('Frequency/Hz','FontSize',7);
        end;
        set(gca,'FontSize',7);
        hold on;
        [a n1]=max(pu);
        plot(f(n1),pu(n1),'r.','MarkerSize',7,'EraseMode','none');
        N=num2str(n1/1024*2500,'%.1f');
        text(f(n1)+40,pu(n1),[N,'Hz'],'FontSize',7,'Color',[1 0 0]);
        hold off;
    end;
end;
%**************************************************************************************************************************************************************
%繪制小波分解各子空間能量的棒圖
for i=m+1:m+round(m/6)+1;
    figure(i);
    for j=1:6;
        if ((6*(i-m-1)+j)<=m);
            R=xiaobofenxi(acc(6*(i-m-1)+j,:),n);
            for k=1:n+1;
                W(k)=sum(R(k,:).^2);
            end;
            totleenergy=sum(W);
            subplot(3,2,j);
            bar(1:n+1,W,0.3);
            colormap hsv;
            set(gca,'XTickLabel',{'a6','d6','d5','d4','d3','d2','d1'}); %重要：此表達式為6層小波分解的格式，若是改變分解尺度，需要調整橫坐標的刻度表達，因本人有些懶得編制通用程序了，所以這樣應付一下
            set(gca,'YTickLabel',{''});
            set(gca,'FontSize',7);
            hold on;
            for ii=1:n+1;
                percent=num2str(100*W(ii)/totleenergy,'%.1f');
                text(ii-0.4,W(ii)+max(W/10),[percent,'%'],'FontSize',7,'color',[0 0 0]); %-0.4,+max(W/10)為標記位置調整，應隨具體情況確定相應調整值
            end;
            hold off;
            ylabel([num2str(6*(i-m-1)+j),' point'],'FontSize',8);
        end;
    end;
end;