time_index1.m

刀具在线监测故障信号的时域频域时频域特征提取matlab程序

%pre-Processing the data: 1,convert to real current 2,substract the zero

%待解决问题1，判断文件内容是否为空，是否需要添加 表头 2，如果查找一个文件夹下的txt文件，并作循环%
clear;
%clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%输入参数
filename='normal_7_1.txt';
startline=570000;
nfft=90000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%固定参数
[c1,c2,c3] = textread(['D:\dataprocess\data\datatxt\',filename], ' %f  %f  %f',  nfft, 'headerlines', startline);
spindlemotor=fopen('D:\dataprocess\data\spindlemotor.txt','at');
servocontrol=fopen('D:\dataprocess\data\servocontrol.txt','at');
servomotor=fopen('D:\dataprocess\data\servomotor.txt','at');

c=[c1,c2,c3];
spindlespeed=500;                                        %主电机转速r/min
hallratio=500;                                           %霍尔传感器变比
ballscrew=5;                                             %5mm/rev
servopulse=7200;                                         %7200pulse/rev
yname={'spindle motor','servo control','servo motor'};   %图像y轴内容显示
fsample=33000;                                           %采样频率
keyfrequency=200;                                        %高低频率的界点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算参数
if findstr(filename,'_7_')>0
    diameter=7;                                              %刀具直径mm
elseif  findstr(filename,'_8_')>0
    diameter=8;                                              %刀具直径mm
elseif  findstr(filename,'_9_')>0
    diameter=9;                                              %刀具直径mm
elseif  findstr(filename,'_10_')>0
    diameter=10.2;                                              %刀具直径mm
elseif  findstr(filename,'_11_')>0
    diameter=11.5;                                              %刀具直径mm
elseif  findstr(filename,'_12_')>0
    diameter=12.1;                                              %刀具直径mm
elseif  findstr(filename,'_13_')>0
    diameter=7;                                              %刀具直径mm
end  

if findstr(filename,'_1.')>0
    feedspeed=500;                                           %进给脉冲周期500us
elseif  findstr(filename,'_2.')>0
    feedspeed=700;                                           %进给脉冲周期700us
elseif  findstr(filename,'_3.')>0
    feedspeed=800;                                           %进给脉冲周期800us
end  

if findstr(filename,'sharp')>0
    state=1;
    resitance=680;                                           %外界电阻欧姆
elseif  findstr(filename,'normal')>0
    state=2;
    resitance=680;                                           %外界电阻欧姆
elseif  findstr(filename,'worn')>0
    state=3;
    resitance=180;                                           %外界电阻欧姆
end                                                          %锋利  1   正常  2   钝  3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
feedspeed=(ballscrew/(servopulse*feedspeed*0.000001))/(spindlespeed/60);        %进给速度mm/rev
cuttingspeed=spindlespeed*pi*diameter/1000;              %切削线速度m/min
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
window1=hamming(32768);                                         %汉明窗
noverlap=16384;
range='onesided';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%指标计算

for i=1:3
    c(:,i)=c(:,i)*(hallratio/resitance);
    c(:,i)=c(:,i)-mean(c(:,i));
    figure(i);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(3,1,1); plot(c(:,i),'.'); ylabel('amplitude /A');title([yname{i},' original signal']);
    %time-domain feature index
    rms=(mean((c(:,i).^2))).^0.5;  % 均方根% 
    peaktime=max(abs(c(:,i))); %最大值
    wavefactor=((mean((c(:,i).^2))).^0.5)/mean(abs(c(:,i))); %波形因数
    kurt=kurtosis(c(:,i)); %峭度
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %frequency-domain feature index
    [Pxx1,f]=pwelch(c(:,i),window1,noverlap,nfft,fsample,range); 
    subplot(3,1,2);plot(f,Pxx1,'.'); ylabel('power welch');title([yname{i},' power spectrum']);
    peakspec=max(Pxx1);    %功率最大值
    energylow=sum(Pxx1(1:keyfrequency));  energyhigh=sum(Pxx1(1:keyfrequency));   %高低频段能量
    barylow=mean(Pxx1(1:keyfrequency)); baryhigh=mean(Pxx1(1:keyfrequency));   %高低频段能量的分布重心
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %time-frequency feature index
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %结果写入到相应的文件中
    
%     if i==1 
%         fprintf(spindlemotor,'%9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f\n',diameter, cuttingspeed, feedspeed, rms, peaktime, wavefactor, kurt,state);
%         fclose(spindlemotor);
%     elseif i==2 
%          fprintf(servocontrol,'%9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f\n',diameter, cuttingspeed, feedspeed, rms, peaktime, wavefactor, kurt,state);
%         fclose(servocontrol);
%     elseif i==3
%          fprintf(servomotor,'%9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f   %9.4f\n',diameter, cuttingspeed, feedspeed, rms, peaktime, wavefactor, kurt,state);
%         fclose(servomotor);
%     end
end