disp(' ')
disp(' mean_filter.m   ver 1.2   March 15, 2013')
disp(' by Tom Irvine   Email: tom@vibrationdata.com')
disp(' ')
disp(' This script performs a mean filter on a time history. ');
disp(' ')
disp(' The input time history must have two columns: time(sec) & amplitude')
disp(' ')
%
%
clear raw_minus_mf;
clear mean_filtered;
clear THM;
clear t;
clear y;
clear yf;
clear a;
clear w;
clear q;
%
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
disp('   3=Excel file ');
file_choice = input('');
%
if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
end
if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
t=THM(:,1);
y=THM(:,2);
%
tmx=max(t);
tmi=min(t);
n = length(y);
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
disp(' ')
disp(' Time Step ');
dtmin=min(diff(t));
dtmax=max(diff(t));
%
out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
out5 = sprintf(' dt     = %8.4g sec  ',dt);
out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
disp(out4)
disp(out5)
disp(out6)
%
disp(' ')
disp(' Sample Rate ');
out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
disp(out4)
disp(out5)
disp(out6)
%
sr=1/dt;
%
            figure(1);
            plot(t,y);
            title(' Raw Signal ');
            grid on;       
%
ncontinue=1;
if(((dtmax-dtmin)/dt)>0.01)
    disp(' ')
    disp(' Warning:  time step is not constant.  Continue calculation? 1=yes 2=no ')
    ncontinue=input(' ');
end
if(ncontinue==1)
%%
%      The number of passes and window size should be determined via
%      trial-and-error.
%%
        a=y;
%
        disp(' ');
        disp(' Select method ');
        disp('  1=specify number of passes & window size');
        disp('  2=specify -3 dB filter frequency. ');
        im=input(' ');        
%
        disp(' ')
%
        maxw=round(double(n)/2);
%
        if(im==1)
            np=input(' Enter number of passes ');
            disp(' ');
            out1=sprintf(' Enter the window size.  Must be odd number. < %d ',maxw);
            disp(out1);
            w=input(' ');
        else
            np=1;
            disp(' Enter primary filter goal ');
            disp('  1=lowpass ');
            disp('  2=highpass ');
            ip=input(' ');
            disp(' ');
            if(ip==1)
              fc=input(' Enter lowpass frequency (Hz) ');
              q=0.400*(sr/fc);
            else
              fc=input(' Enter highpass frequency (Hz) ');
              q=0.728*(sr/fc);
            end
            w= 2*floor(q/2)+1;
            if(w>maxw)
                w= 2*floor(q/2)-1; 
            end
            if(w<3)
                w=3;
            end            
            out1=sprintf('\n  window size=%d \n',w);
            disp(out1);
        end
%
        k=fix(double(w-1)/2.);
%
        last=n;
        for(m=1:np)
%
            for(i=1:last)
%
                ave=0.;
                n=0;
%
                for(j=(i-k):(i+k))
%
                    if(j>=1 && j<=last )
%
                        ave=ave+a(j);
                        n=n+1;
                    end	
                end
                if(n>1)
                    a(i)=ave/double(n);
                end
            end
  
        end
%
            yf=y-a;
%
            vvmax=max(yf);
            vvmin=min(yf);
            disp(' ');
            disp(' Raw Signal with Mean Filtered Signal Removed');          
            out1=sprintf(' max=%8.4g  min=%8.4g ',vvmax,vvmin);
            disp(out1);    
 %
            figure(2);
            plot(t,yf);
            title(' Raw Signal with Mean Filtered Signal Removed, mf highpass ');
            grid on;
            raw_minus_mf(:,1)=t';
            raw_minus_mf(:,2)=yf;        
 %           
            figure(3);
            plot(t,a);
            title(' Mean Filtered signal, mf lowpass');
            grid on;  
            mean_filtered(:,1)=t';
            mean_filtered(:,2)=a';   
%
            vvmax=max(a);
            vvmin=min(a);
            disp(' ');
            disp(' Mean Filtered Signal');          
            out1=sprintf(' max=%8.4g  min=%8.4g ',vvmax,vvmin);
            disp(out1);
            disp(' ');
%
            clear mf_lowpass;
            clear mf_highpass;
            mf_lowpass=mean_filtered;
            mf_highpass=raw_minus_mf;
%
            disp(' Output Files: ');
            disp('   mf_lowpass  ');
            disp('   mf_highpass ');
%
end