disp(' ')
disp(' Butterworth_filter.m ')
disp(' ver 2.5  July 1, 2013 ')
disp(' by Tom Irvine  Email: tom@vibrationdata.com ')
disp(' ')
disp(' Butterworth filter, sixth-order, infinite impulse response,');
disp(' cascade with refiltering option for phase correction ');                 
disp(' ')
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude ')
disp(' ')
%
close all;
%
clear THM;
clear tim;
clear y;
clear yt;
clear a;
clear b;
%
fig_num=1;
%
[tim,y,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history();
%
figure(fig_num);
fig_num=fig_num+1;
plot(tim,y);
grid on;
xlabel('Time(sec)');
%
ns = n;
%
if(ncontinue==1)
%
fl=0.;
fh=0.;
iband=0;  
%
iband= 0;
while((iband ~= 1) & (iband ~= 2) & (iband ~= 3))  
    disp(' ');
    disp(' Enter filter type ');
    disp(' 1=lowpass  2=highpass  3=bandpass ');
    iband=input('');
end
%
f=0.;
fl=0.;
fh=0.;
%
if(iband == 1)
          fl=input(' Enter lowpass frequency (Hz) ');
          f=fl;
end
if(iband == 2)
          fh=input(' Enter highpass frequency (Hz) ');      
          f=fh;
end
if(iband==1 | iband==2)
    if( f > 0.5*sr) 
	    disp('  error: cutoff frequency must be < Nyquist frequency ');    
	    iflag=999;
    end
end
if(iband == 3)
		  disp(' This bandpass filter is implemented as a ');
		  disp(' highpass filter and lowpass filter in series. ');
%    
		  disp(' Enter lower frequency (Hz) ');
		  fh=input(' '); 
		  disp(' Enter upper frequency (Hz) ');
		  fl=input(' ');
          f=fh;
%
    if( fl > 0.5*sr | fh > 0.5*sr) 
	    disp('  error: cutoff frequency must be < Nyquist frequency ');    
	    iflag=999;
    end
%
end
disp(' ');
disp('  Refiltering for phase correction? ');
disp(' 1=yes (appropriate for transient signal)');
disp(' 2=no  (appropriate for steady-state signal)');
iphase=input(' ');
%
if(iphase ~= 1 & iphase ~=2 )
    iphase=1;
end    
%
freq=f;
%
iflag=1;
%
%****** calculate coefficients *******
%
if(iflag ~= 999 & iband ~=3)
%
    [a,b,iflag] = filter_coefficients(f,dt,iband,iflag);
%    
end
%
if(iflag < 900 )
%
    if(iband == 1 | iband ==2)  % lowpass or highpass
%		 
         if(iphase==1)   % refiltering
			[y]=apply_filter(y,iphase,ns,a,b);
			[y]=apply_filter(y,iphase,ns,a,b); 
	      else
			[y]=apply_filter(y,iphase,ns,a,b);	
          end
%     
    end
%
    if(iband == 3)  % bandpass
 %     
		 f=fh;
		 freq=f;
%
		 disp(' Step 1');
		 iband=2;
         [a,b,iflag] = filter_coefficients(f,dt,iband,iflag);
 %		 
         if(iphase==1) % refiltering
			[y]=apply_filter(y,iphase,ns,a,b);
			[y]=apply_filter(y,iphase,ns,a,b);  
	      else
			[y]=apply_filter(y,iphase,ns,a,b);	
          end
%     
		 f=fl;
		 freq=f;
%
		 disp(' Step 2');
		 iband=1;
         [a,b,iflag] = filter_coefficients(f,dt,iband,iflag);
%		 
         if(iphase==1) % refiltering 
			[y]=apply_filter(y,iphase,ns,a,b);
			[y]=apply_filter(y,iphase,ns,a,b);  
	      else
			[y]=apply_filter(y,iphase,ns,a,b);	
          end
%          
    end
%
    mu=mean(y);
    sd=std(y);
    mx=max(y);
    mi=min(y);
    rms=sqrt(sd^2+mu^2);
    out1 = sprintf('\n Filtered Signal:  mean = %8.4g    std = %8.4g    rms = %8.4g ',mu,sd,rms);
    disp(out1);
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(tim,y);
    xlabel('Time (sec)');
    out5 = sprintf(' Filtered Signal   Overall Level = %8.3g ',rms);
    title(out5);
    grid;
    zoom on;
%
    disp(' ');
    disp(' Write filtered signal to ASCII file? ')
    pchoice=input(' 1=yes 2=no ');
    if pchoice == 1 
%
      [writefname, writepname] = uiputfile('*','Save filtered data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [tim, y ];
	  fid = fopen(writepfname,'w');
      clear length;
      n=length(tim);
      for i=1:n
         fprintf(fid,'%g \t %g \n',tim(i),y(i));
      end  
	  fclose(fid);
  end  
%
    clear filtered_data;
    filtered_data(:,1)=tim;
    filtered_data(:,2)=y;
    disp(' ');
    disp(' Output matrix is named: filtered_data ');
    clear y;
    clear yt;
%   
else
    disp(' ')
    disp('  Abnormal termination.  No output file generated. ');
end
end