disp(' ')
disp(' interp_fourier.m   ver 1.0   September 3, 2011')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
disp(' This program interplates a time history via its Fourier transform.')
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 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 PSD;
%
%  need mean removal and hanning
%
disp(' ')
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=double(THM(:,1));
y=double(THM(:,2));
%
tmx=max(t);
tmi=min(t);
mu=mean(y);
%
n = length(y)
%
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)
%
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 output option ');
%% choice=input(' 1=plot only   2=plot & output text file  ' );
%% disp(' ')
choice=1;
%
disp(' ')
disp(' mean removal? ');
mr_choice=input(' 1=yes   2=no  ' );
disp(' ')
%
if(mr_choice==1)
    y=y-mu;
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
%
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 \t %11.4e \t %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 \t %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
%
disp(' ');
disp(' Plot Fourier Transform ');
%
    figure(1); 
    plot(ff,zz);
    grid on;
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');    
    out5 = sprintf(' Fourier Transform Magnitude ');
    title(out5);  
    zoom on
%
    clear FT;
    FT=[ff',zz'];
%%%%
    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 complex_FT;
complex_FT=[freq', f_real', f_imag'];
%
ipp=2;
%
   disp('   complex Fourier transform saved as: complex_FT ');
   disp(' Fourier transform magnitude saved as: FT ');
%
disp(' ');
mi=input(' Enter sample rate multiplier integer ');
mi=round(mi);
nsr=mi*sr;
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear int_real;
clear int_imag;
clear length;
%
clear t;
no=length(f_real);
nnn=(mi*no);
dt=dt/mi;
for(i=1:nnn)
    t(i)=dt*(i-1);
end
%
%% out1=sprintf(' nnn=%d  mi=%d  no=%d  \n',nnn,mi,no);
%% disp(out1);
%
int_real=zeros(nnn,1);
int_imag=zeros(nnn,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Render Symmetric
%
nd=nnn;
%
if(mod(no,2)==0)
% even
    disp(' even ');
    mid=no/2;
    midp1=mid+1;
    midm1=mid-1;
%
    int_real(1:midp1)=f_real(1:midp1);
    int_imag(1:midp1)=f_imag(1:midp1);
    int_real(midp1)=int_real(midp1)/2
    int_imag(midp1)=int_imag(midp1)/2
%
    for(i=2:midp1)
        int_real(nd-i+2)= int_real(i);
        int_imag(nd-i+2)=-int_imag(i);
    end
else
% odd
    disp(' odd ');
    nx=(no-1)/2;
    nxp1=nx+1;
    int_real(1:nxp1)=f_real(1:nxp1);
    int_imag(1:nxp1)=f_imag(1:nxp1);    
    for(i=2:nxp1)
        int_real(nd-i+2)= int_real(i);
        int_imag(nd-i+2)=-int_imag(i);
    end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear f;
clear y_real;
clear y_imag;
clear t_real;
clear t_imag;
clear sum_c;
clear sum_s;
clear ccc;
clear sss;
%
y_real=int_real;
y_imag=int_imag;
%
k = nnn;
%
%
disp(' ')
%
iflag =0;
tpi=2.*pi;
%
disp(' ');
progressbar % Create figure and set starting time
%
for  n=0:(k-1)
  progressbar(n/k) % Update figure    
%
  t_real(n+1)=0.;
  t_imag(n+1)=0.;  
%  
%  arg = linspace(0,tpi*n,k);
  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;
%
  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);
%
if choice==2
    fclose(fid1);     
end
%
disp(' ');
disp(' Plot Time History ');
%
    figure(2); 
    plot(t,t_real);      
    xlabel('Time (sec)');
    ylabel('Amplitude');    
    out5 = sprintf(' Inverse Fourier Transform Real');
    title(out5);  
%
    figure(3); 
    plot(t,t_imag);      
    xlabel('Time (sec)');
    ylabel('Amplitude');    
    out6 = sprintf(' Inverse Fourier Transform Imaginary');
    title(out6);  
%
    zoom on
%
    figure(4); 
    plot(t,t_real,THM(:,1),THM(:,2));  
    title('Interpolated Time History');
    legend ('interpolated','raw');  
    xlabel('Time (sec)');
    ylabel('Amplitude')
    grid on;
%%%%
%
sz=size(t);
if(sz(2)>sz(1))
    t=t';
end
%
sz=size(t_real);
if(sz(2)>sz(1))
    t_real=t_real';
end
clear interpolated_th;
interpolated_th=[t t_real];
disp(' ');
disp(' Output array = interpolated_th ');