%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