%
%   integrate_th_mf.m  ver 1.0  March 15, 2013
%
function[v,figure_number]=integrate_th_mf(nba,ist,t,v,figure_number)   
%
    disp(' ');
    if(nba==1)
%%       disp(' remove trend from input signal via mean filtering?  1=yes  2=no ');
        isk=1;
    else
        disp(' remove trend from output signal via mean filtering?  1=yes  2=no ');
        isk=input(' ');
    end
%
    n=max(size(v));
    sr=n/(t(n)-t(1));
%
    if(isk==1)
%%
%      The number of passes and window size should be determined via
%      trial-and-error.
%%
        a=v;
%%
%
%
        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;
%
            fc=input(' Enter highpass frequency (Hz) ');
            q=0.728*(sr/fc);
%
            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
%
        vv=v-a;
        v=vv;
        figure(figure_number);
        figure_number=figure_number+1;
%        
        plot(t,v);     
 %       
        if(nba==1)
            if(ist==1)
                title(' Acceleration Mean Filtered Signal prior to Integration');
            else
                title(' Velocity Mean Filtered Signal prior to Integration');                
            end
        else
            if(ist==1)
                title(' Velocity Mean Filtered Signal after Integration');
            else
                title(' Displacement Mean Filtered Signal after Integration');                
            end 
        end 
%
        grid on;
        xlabel('Time (sec)');
    end