disp(' ');
disp(' remove_noise_frequency.m  version 1.1   May 27, 2009');
disp(' By Tom Irvine   Email:  tomirvine@aol.com ');
disp(' ');
disp(' This program removes noise frequencies from data. ');
%
clear t;
clear tt;
clear a;
clear amp;
clear cleaned;
clear damp;
clear delay;
clear f;
clear p;
clear fc;
clear fl;
clear fu;
clear THM;
clear dt;
clear sr;
clear syn;
clear num2;
clear s_amp_origin;
hold off;
clear figure(1);
%    
tp=2.*pi;    
%%%
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude ')
disp(' ')
%
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
file_choice = input('');
%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);     
%    disp(' Enter the input filename ');
%    filename = input(' ','s');
    fid = fopen(filename,'r');
    THM = fscanf(fid,'%g %g',[2 inf]);
    THM=THM';
else
    THM = input(' Enter the matrix name:  ');
end
%
s_amp=THM(:,2);
s_amp_origin=s_amp;
s_tim=THM(:,1);
%
clear progressbar;
%
figure(1);
plot(s_tim,s_amp);
title(' Raw Data ');
xlabel(' Time (sec) ');
grid on;
%
n = length(s_amp); 
clear syn;
syn=zeros(1,n);
%
tmx=max(s_tim);
tmi=min(s_tim);
%
out1 = sprintf('\n  %d samples \n',n);
disp(out1)
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
out1 = sprintf(' SR  = %g samples/sec    dt = %g sec \n',sr,dt);
out3 = sprintf(' start  = %g sec    end = %g sec            ',tmi,tmx);
disp(out1)
disp(out3)
%
n1=1;
T=1./60.;
T=T*20;
ndelta=round(T/dt);
ndelta_old=ndelta;
% ndelta=10000;
%
nlast=fix(n/ndelta);
%
residual=n-fix(nlast*ndelta);
%
iflag=0;
%
if(nlast>1 && residual>1)
    nlast=nlast+1;
    iflag=1;
end
%
if(nlast==0)
    nlast=1;
    ndelta=n;
end
%
fmax=0.;
%
disp(' ');
nt=input(' Enter the number of trials per frequency ');
disp(' ');
%
disp(' Enter the noise frequency (Hz)');
fc(1)=input(' ');
np=input(' Enter the number of passes ');
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for(kjk=1:np)
%
clear a;
clear t;
clear amp;
clear tim;
clear temp;
%
n1=1;
%
j=0;
%
fu=fc*1.04;
fl=fc*0.96;
%
nfr=max(size(fl));
%
fmax=max(fu);
%
if(kjk>1)
    s_amp=cleaned;
end
%
progressbar % Create figure and set starting time
%
for(ijk=1:nlast)
%
    progressbar(ijk/nlast) % Update figure    
%  
    ndelta=ndelta_old;
    if(iflag==1 && ijk==nlast)
        ndelta=residual;
    end
%
n2=n1+ndelta-1;
amp=s_amp(n1:n2)';
tim=s_tim(n1:n2)';
np=length(tim);
n=np;
%
tt=tim-tim(1);
%
%% out1 = sprintf(' np  = %d  \n',np);
%% disp(out1)
%
    t=tim;
    a=amp;
    dur=max(t)-min(t);
    num2=length(t);
%
    clear THM;
%
    drate=fix(sr/(fmax*16.));
%
    if(drate <2)
        drate=1;
    end    
    if(drate > 10)
        drate=10;
    end    
%
    for(ie=1:nfr)
%
        if(ijk==1)
            x1r(ie)=0.;
            x2r(ie)=0.;
            x3r(ie)=0.;
        end
%
        out4 = sprintf('\n frequency case %d \n',ie);
        disp(out4)    
%
        flow=fl(ie);
        fup=fu(ie);
%
   	    [a,x1r,x2r,x3r,error_rth,syna]=sfa_engine_60(a,tim,drate,num2,flow,fup,ie,nt,x1r,x2r,x3r,ijk,nlast);  
%
    end
%
    disp(' ');
    disp(' Results');
    disp(' ');
    disp(' Case   Amplitude   fn(Hz)   Phase(rad)  ');
%
    for(ie=1:nfr)
        out4 = sprintf('  %d  %10.4f  %10.4f  %10.4f  ',ie,x1r(ie),x2r(ie)/tp,x3r(ie));
        disp(out4)      
    end    
%
    clear tt;
    tt=tim-tim(1);
%
    ssaa=size(amp);
    if(ssaa(2)>ssaa(1))
        amp=amp';
    end
%
    for(ie=1:nfr)
%	    
   		syn(n1:n2)=syn(n1:n2) + x1r(ie)*sin(x2r(ie)*tt(1:(n2-n1)+1)-x3r(ie));
%       
    end
%
    n1=n2+1;
%%%%%%%%%%%%%%%
end
cleaned=s_amp_origin-syn';    
end
%%%%%%%%%%%%%%%
%
figure(2)
%
tt=tt+t(1);
%
tq=s_tim(1:max(size(cleaned)));
plot(tq,cleaned);
%
xlabel(' Time (sec) ');
title(' Cleaned Data ');
grid on;
%
clear cleaned_data;
cleaned_data=[tq,cleaned];
%
noise=[tq,syn'];
%
figure(3)
plot(noise(:,1),noise(:,2));
title(' Noise ');
xlabel(' Time (sec) ');