disp(' ')
disp(' fourier_filter_Bessel.m   ver 1.1   June 27, 2007')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
disp(' This program calculates the Fourier transform of a time history')
disp(' which is pre-loaded into Matlab.  It then applies a Bessel ')
disp(' filter in the frequency domain. ')
disp(' ')
disp(' The time history should have less than 10000 points')
disp(' ')
disp(' The time history must have two columns: time(sec) & amplitude(units)')
disp(' ')
%
clear figure(1);
clear figure(2);
clear figure(3);
clear figure(4);
clear figure(5);
clear figure(6);
clear t;
clear y;
clear z;
clear tmx;
clear tmi;
clear n;
clear dt;
clear df;
clear sr;
clear f_real;
clear f_imag;
clear freq;
clear ff;
clear zz;
clear mu;
clear hw;
clear ms;
clear rms;
clear THM;
clear H;
clear MFT;
tpi=2.*pi;
%
%  need mean removal and hanning
%
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
%
t=double(THM(:,1));
y=double(THM(:,2));
%
figure(1);
plot(t,y);
xlabel(' Time(sec) ');
%
tmx=max(t);
tmi=min(t);
tfirst=tmi;
mu=mean(y);
%
n = length(y);
%
out3 = sprintf(' Number of samples  n = %d  ',n);
disp(out3)
%
if(n>10000)
    out1=sprintf('\n Warning: long calculation time required for high sample number.\n');
    disp(out1)
end
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
df=1./(tmx-tmi);
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)
%
out7 = sprintf('    df  = %8.4g Hz  \n',df);
disp(out7)
%
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)
%
%  Output options
%
disp(' ')
disp(' Select Fourier Transform output option ');
choice=input(' 1=plot only   2=plot & output text file  ' );
disp(' ')
%
if choice==2
   disp(' Enter the complex output filename ');
   filename1 = input(' ','s');
%   
   disp(' ')  
   disp(' Enter the one-sided, full amplitude filename ');
   filename2 = input(' ','s'); 
%
   fid1 = fopen(filename1,'w');
   fid2 = fopen(filename2,'w');   
%   
end
%
disp(' ')
disp(' mean removal? ');
mr_choice=input(' 1=yes   2=no  ' );
disp(' ')
%
if(mr_choice==1)
    y=y-mu;
end
disp(' ')
disp(' Hanning window? ');
hw_choice=input(' 1=yes   2=no  ' );
disp(' ')
%
   ae=((8/3)^(1/2));
%    
if(hw_choice==1)
    progressbar % Create figure and set starting time 
    for i=0:(n-1)
        progressbar(i/(n-1)) % Update figure  
        hw(i+1) =sin( (i*pi/n) )^2 ;
    end
     y=(y.*hw')*ae;
end
%
out1 = sprintf('\n  %d samples \n',n);
disp(out1)
%
out1 = sprintf(' SR  = %g samples/sec    dt = %g sec \n',sr,dt);
disp(out1)
%
iflag =0;
tpi=2.*pi;
%%%
ms=0;
jk=0;
kk=1;
%
disp(' ');
progressbar % Create figure and set starting time
%
if(n<=2)
    disp(' Error ');
    n
end
%
for  k=0:(n-1)
  progressbar(k/n) % Update figure    
%
  clear sum_c;
  clear sum_s;
  clear ccc;
  clear sss;
%
  f_real(k+1)=0.;
  f_imag(k+1)=0.;  
%  
  arg = linspace(0,tpi*k*(n-1)/n,n);
%
  ccc=y.*cos(arg)';
  sss=y.*sin(arg)';
%
  sum_c=sum(ccc);
  sum_s=sum(sss);
%
  f_real(k+1)=f_real(k+1)+sum_c;
  f_imag(k+1)=f_imag(k+1)+sum_s;
%
  freq(k+1)=k*df;
%       
  f_real(k+1)=f_real(k+1)/n;
  f_imag(k+1)=f_imag(k+1)/(-n);
%
  z(k+1)=((f_real(k+1)^2) + (f_imag(k+1)^2))^(1/2);
%  
  if choice==2
        fprintf(fid1,'%11.4e %11.4e %11.4e \n',freq(k+1), f_real(k+1), f_imag(k+1) );      
  end    
%		  
		  if(k<(n/2))  
%              
				 if(k > 0)			 
                    z(k+1)=2.*z(k+1);  

                 end
                 if choice==2
                    fprintf(fid2,'%11.4e %11.4e \n',freq(k+1), z(k+1) );
                 end
                 zz(kk)= z(k+1);
                 ff(kk)=freq(k+1);
                 ms=ms+(zz(kk))^2;     
                 kk=kk+1;
 %                
          end
  end   
  progressbar(1);
%
  if choice==2
       fclose(fid1);
       fclose(fid2);     
  end
%
if(length(freq)<=1)
    disp(' Error')
    freq
end
%
disp(' ');
disp(' Plot Fourier Transform ');
%
    figure(2); 
    plot(ff,zz);      
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');    
    out5 = sprintf(' Fourier Transform Magnitude ');
    title(out5);  
    grid on
    zoom on
%%%%
    rms=ms^(1/2);
%
    out1 = sprintf('\n  Fourier transform overall level = %g \n',rms);
    disp(out1)
%
    fmax=0;
    zmax=0;
    for(i=1:length(ff))
        if(zz(i)>zmax)
            zmax=zz(i);
            fmax=ff(i);
        end
    end
%
    out5 = sprintf('\n Peak occurs at %10.5g Hz ',fmax);
    disp(out5)
%
    out1=sprintf('\n Calculation complete. \n');
    disp(out1)
%%%%
end
%
clear FT;
clear complex_FT;
complex_FT=[freq', f_real', f_imag'];
FT=[freq',z'];
%
disp('   complex Fourier transform saved as: complex_FT ');
disp(' Fourier transform magnitude saved as: FT ');
%
%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select Filter Type ')
disp(' 1=lowpass   2=highpass  3=bandpass')
typ=input(' ');
disp(' ');
if(typ==2)
    hpf=input(' Enter highpass frequency(Hz) ');
    fc=hpf;
end
if(typ==1)
    lpf=input(' Enter lowpass frequency(Hz) ');
    fc=lpf;
end
if(typ==3)
    hpf=input(' Enter lower frequency(Hz) ');
    fc=hpf;
%
    lpf=input(' Enter upper frequency(Hz) ');
    fc=lpf;
end
%
L = input(' Enter Filter Order ');
%
n=length(freq)
%
a = n/2;
whole = floor(a);
part = a-whole;
%
if(part > 0) % odd
    oe_type=1;
    nyq=fix(n/2.)+1
else
    oe_type=2;
    nyq=(n/2)  
end
%
%
clear k;
clear H;
clear fH;
%
disp(' ')
disp(' Form Transfer Function ');
%%%
if(typ==1)
    [fH,H]=Bessel_Transfer(1,fc,freq,nyq,L);
end
if(typ==2)
    [fH,H]=Bessel_Transfer(2,fc,freq,nyq,L);
end
if(typ==3)
    clear G;
    fc=lpf;
    [fH,H]=Bessel_Transfer(1,fc,freq,nyq,L);  % lowpass filter
    G=H;
end
%
disp(' Implement Filter ');
clear G;
%
MFT=complex(f_real',f_imag');
for(i=1:nyq)
    Hfc=H(i)*MFT(i);
    f_real(i)=real(Hfc);
    f_imag(i)=imag(Hfc);
    G=H;
end  
if(typ==3)
    clear H;
    fc=hpf;
    [fH,H]=Bessel_Transfer(2,fc,freq,nyq,L);  % highpass filter
    MFT=complex(f_real',f_imag');
    for(i=1:nyq)
        Hfc=H(i)*MFT(i);
        f_real(i)=real(Hfc);
        f_imag(i)=imag(Hfc);
        H(i)=H(i)*G(i);
    end 
end
%
disp(' force symmetry ')
%
if(oe_type==1) % odd
    for(i=1:nyq-1)   
        f_real(n+1-i)= f_real(i+1);
        f_imag(n+1-i)=-f_imag(i+1);  
    end
end
if(oe_type==2)
    for(i=1:(nyq-1))   
        f_real(n+1-i)= f_real(i+1);
        f_imag(n+1-i)=-f_imag(i+1);  
    end
end
%
figure(3);
plot(fH,abs(H));
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
grid on;
xlabel(' Frequency (Hz) ');
ylabel(' Magnitude ');
if(typ==1)
    out1=sprintf(' Bessel Lowpass Filter Order=%d fc=%7.3g Hz ',L,lpf);
end
if(typ==2)
    out1=sprintf(' Bessel Highpass Filter Order=%d fc=%7.3g Hz ',L,hpf);
end
if(typ==3)
    out1=sprintf(' Bessel Bandpass Filter Order=%d  %7.3g to %7.3g Hz ',L,hpf,lpf);
end
title(out1);
%
MFT=complex(f_real',f_imag');
complex_FT=[freq', f_real', f_imag'];
figure(4);
plot(freq,abs(MFT));
xlabel(' Frequency (Hz) ');
ylabel(' Magnitude ');
title(' Fourier Magnitude with Filtering Applied ');
%
%%%%%%%%%%%%%%%%%%%%%%%
%
clear t;
clear y;
clear z;
clear tmx;
clear n;
clear dt;
clear df;
clear sr;
clear f;
clear t_real;
clear t_imag;
clear ff;
clear zz;
clear mu;
clear hw;
clear ms;
clear rms;
clear THM;
%
THM=complex_FT;
%
f=double(THM(:,1));
y_real=double(THM(:,2));
y_imag=double(THM(:,3));
%
k = length(f);
%
if(k>10000)
    out1=sprintf('\n Warning: long calculation time required for high sample number.\n');
    disp(out1)
end
%
dt=1./max(f);
%
%  Output options
%
disp(' ')
disp(' Select Inverse Fourier Transform output option ');
choice=input(' 1=plot only   2=plot & output text file  ' );
disp(' ')
%
if choice==2
   disp(' Enter the time history output filename ');
   filename1 = input(' ','s');
   fid1 = fopen(filename1,'w');
end
%    
iflag =0;
%
disp(' ');
progressbar % Create figure and set starting time
%
for  n=0:(k-1)
  progressbar(n/k) % Update figure    
%
  clear sum_c;
  clear sum_s;
  clear ccc;
  clear sss;
%
  t_real(n+1)=0.;
  t_imag(n+1)=0.;  
%  
  arg = linspace(0,tpi*n*(k-1)/k,k);   
%%
  ccc=y_real.*cos(arg)'-y_imag.*sin(arg)';
  sss=y_real.*sin(arg)'+y_imag.*cos(arg)';
%
  sum_c=sum(ccc);
  sum_s=sum(sss);
%
  t_real(n+1)=t_real(n+1)+sum_c;
  t_imag(n+1)=t_imag(n+1)+sum_s;
%
  t(n+1)=n*dt;
%
  if choice==2
        fprintf(fid1,'%11.4e %11.4e %11.4e \n',t(n+1), t_real(n+1), t_imag(n+1) );      
  end    
%		 
  end   
  progressbar(1);
%
end
if choice==2
    fclose(fid1);     
end
%
disp(' ');
disp(' Plot Time History ');
%
    figure(5); 
    plot(t,t_imag);      
    xlabel('Time (sec)');
    ylabel('Amplitude');    
    out6 = sprintf(' Inverse Fourier Transform Imaginary');
    title(out6);  
%
    figure(6); 
    plot(t,t_real);      
    xlabel('Time (sec)');
    ylabel('Amplitude');    
    out5 = sprintf(' Inverse Fourier Transform Real');
    title(out5);  
    zoom on
%
t=t+tfirst;
inv_real=[t',t_real'];
inv_imag=[t',t_imag'];
filtered_data=inv_real;
disp(' ')
disp(' Real time history renamed as "filtered data" ')