disp(' ')
disp(' mean_filter.m   ver 1.1   July 15, 2010')
disp(' by Tom Irvine   Email: tomirvine@aol.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) & accel(G)')
disp(' ')
%
%
clear raw_minus_mf;
clear mean_filtered;
clear THM;
clear t;
clear y;
clear yf;
clear a;
%
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)
%
            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(' ')
        np=input(' Enter number of passes ');
%
        disp(' ');
        out1=sprintf(' Enter the window size.  Must be odd number. < %d ',round(double(n)/2.));
        disp(out1);
        w=input(' ');
%
        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 ');
            grid on;
            raw_minus_mf(:,1)=t';
            raw_minus_mf(:,2)=yf;        
 %           
            figure(3);
            plot(t,a);
            title(' Mean Filtered signal ');
            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(' ');
%
            disp(' Output Files: ');
            disp('   mean_filtered ');
            disp('   raw_minus_mf ');
%
end