qingxiexiaobofenxi.m

自己编写的比较全面的故障诊断matlab函数程序

%**************************************************************************************************************************************************************
%    该程序所要达到的目的是对发动机振动加速度信号进行二进制小波分解，绘制各层信号的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;