disp(' ')
disp(' fourier.m   ver 4.0   November 13, 2012')
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.  ')
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(' ')
%
close all;
%
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 sum;
clear THM;
clear PSD;
%
fig_num=1;
%
%  need mean removal and hanning
%
[t,y,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history();
%
x_label=sprintf('Time(sec)');
y_label=sprintf('Amplitude');
t_string=sprintf('Time History');
[fig_num]=plot_TH(fig_num,x_label,y_label,t_string,[t y]);
%
dur=tmx-tmi;
df=1/dur;
mu=mean(y);
%
if(ncontinue==1)
%
%  Output options
%
    disp(' ')
    disp(' Select 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(' ')
%
    disp(' ')
    disp(' Select Window ')
    hw_choice = input(' 1=Rectangular 2=Hanning  ');
    disp(' ')
%
    if(mr_choice==1 || hw_choice==1)
        [y]=mean_removal_Hanning(y,mr_choice,hw_choice);
    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;
%%%
    jk=0;
%
    disp(' ');
    progressbar % Create figure and set starting time
%
    nhalf=floor(n/2);
%
    f_real=zeros(1,n);
    f_imag=zeros(1,n);
    z=zeros(1,n);
    freq=zeros(1,n);
    zz=zeros(1,nhalf);
%
    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    
%		  
    end   
    progressbar(1);
%
    zz(1)=z(1);
    zz(2:nhalf)=2*z(2:nhalf);
    ms=0.5*sum(zz(1:nhalf).^2);  % include peak to rms conversion
%    
    ff=freq(1:nhalf);
%
    if choice==2
        for k=1:nhalf
            fprintf(fid2,'%11.4e \t %11.4e \n',freq(k), zz(k) );
        end  
    end
%
    if choice==2
       fclose(fid1);
       fclose(fid2);     
    end
%
    disp(' ');
    disp(' Plot Fourier Transform ');
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(ff,zz);
    grid on;
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');    
    out5 = sprintf(' Fourier Transform Magnitude ');
    title(out5);  
    zoom on
%
    clear FT;
    ff=fix_size(ff);
    zz=fix_size(zz);    
    FT=[ff zz];
%
    rms=sqrt(ms);
%
    out1 = sprintf('\n  Fourier transform overall level = %g \n',rms);
    disp(out1)
%
    [xmax,fmax]=find_max(FT);
%
    out5 = sprintf('\n Peak occurs at %10.5g Hz ',fmax);
    disp(out5)
%
    out1=sprintf('\n Calculation complete. \n');
    disp(out1);
%%%%
end
%
freq=fix_size(freq);
f_real=fix_size(f_real);
f_imag=fix_size(f_imag);
clear complex_FT;
complex_FT=[freq, f_real, f_imag];
%
disp(' ');
disp(' Plot PSD? ')
disp(' 1=yes 2=no');
ipp=input(' ');
%
if(ipp==1)
    clear PS;
    zz=zz/sqrt(2);
    PS=zz.*zz/df;
    rms=sqrt(sum(PS)*df);
    PS=fix_size(PS);
    PSD=[ff,PS];
%
    power_spectral_density=PSD;
    clear PSD;
    %
    t_string=...
    sprintf(' Power Spectral Density  Overall Level = %6.3g RMS ',rms); 
    x_label=sprintf(' Frequency (Hz)');
    ylab=sprintf(' PSD (amplitude^2/Hz)');
    [fig_num]=plot_PSD(fig_num,x_label,ylab,t_string,power_spectral_density);
%    
end
%
   disp('   complex Fourier transform saved as: complex_FT ');
   disp(' Fourier transform magnitude saved as: FT ');
if(ipp==1)
   disp('      Power spectral density saved as: power_spectral_density ');
end