disp(' ')
disp(' fourier_filter.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 Butterworth ')
disp(' sixth-order 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
%
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]=Butterworth_Transfer(1,fc,freq,nyq);
end
if(typ==2)
    [fH,H]=Butterworth_Transfer(2,fc,freq,nyq);
end
if(typ==3)
    clear G;
    fc=lpf;
    [fH,H]=Butterworth_Transfer(1,fc,freq,nyq);  % 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]=Butterworth_Transfer(2,fc,freq,nyq);  % 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(' Butterworth Lowpass Filter 6th order fc=%7.3g Hz ',lpf);
end
if(typ==2)
    out1=sprintf(' Butterworth Highpass Filter 6th order fc=%7.3g Hz ',hpf);
end
if(typ==3)
    out1=sprintf(' Butterworth Bandpass Filter 6th order %7.3g to %7.3g Hz ',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" ')