disp(' ')
disp(' full_FFT.m ')
disp(' ver 2.4 July 9, 2009 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the one-sided, full amplitude FFT ')
disp(' of a time history ')
disp(' ')
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude ')
disp(' ')
%
clear n;
clear Y;
clear amp;
clear tim;
clear N;
clear df;
clear dt;
clear Mag;
clear full;
clear freq;
clear complex_FFT;
clear THM;
clear mu;
clear mean;
%
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
%
%N=4096;
%
amp=double(THM(:,2));
tim=double(THM(:,1));
n = max(size(amp));
%
for(i=1:18)
    if( 2^i > n )
        break;
    end
    N=2^i;
end
%
out4 = sprintf(' time history length = %d ',n);
disp(out4)
disp(' ');
out5 = sprintf(' samples used for FFT = %d',N);
disp(out5)
%
%disp(' mean values ')
%
mu=mean(amp);
sd=std(amp);
mx=max(amp);
mi=min(amp);
rms=sqrt(sd^2+mu^2);
tmx=max(tim);
tmi=min(tim);
dt=(tmx-tmi)/(n-1);
%
sr=1/dt;
%
disp(' ')
disp(' Time Step ');
dtmin=min(diff(tim));
dtmax=max(diff(tim));
%
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)
%
disp(' ')
disp(' Select Window ')
h_choice = input(' 1=Rectangular 2=Hanning  ');
%
if(h_choice==2)
%
    ae=sqrt(8./3.);
	disp(' Hanning window ');
	for(i=1:N)
	    amp(i)=amp(i)*( ae*(sin( (i*pi/N) )^2) );       
    end
end
%
disp(' ')
disp(' begin FFT ')
disp(' ')
Y = fft(amp,N);
%
clear complex_FFT;
complex_FFT=zeros(N,3);
clear FF;
%
df=1./(N*dt);
out4 = sprintf(' df  = %8.4g Hz  ',df);
disp(out4)
%
FF=linspace(0,df*(N-1),N);
%
complex_FFT(:,1)=FF';
complex_FFT(:,2)=real(Y)/N;
complex_FFT(:,3)=imag(Y)/N;
%
freq(1)=0.;
%
disp(' scale for full FFT ')
freq(1)=0;
for(i=2:N/2)
   aa=real(Y(i));
   bb=imag(Y(i));
   Mag(i) =sqrt(aa^2+bb^2);
%
   phase(i) =atan2(bb,aa);
   freq(i)=(i-1)*df;
end
%
clear magnitude_FFT;
magnitude_FFT=[freq',Mag'];
%
full=2.*Mag/N;
full(1)=0.;
phase = phase*180./pi; 
%
fmax=0;
zmax=0;
for(i=1:max(size(freq)))
        if(full(i)>zmax)
            zmax=full(i);
            fmax=freq(i);
        end
end
%
out5 = sprintf('\n Peak occurs at %10.5g Hz \n',fmax);
disp(out5)
%
disp(' Plot the FFT magnitude? ')
choice = input(' 1=yes 2=no ');
%
if(choice == 1)
    figure(1);     
    plot(freq,full)
    xlabel(' Frequency (Hz)');
    ylabel(' Magnitude');
    title(' FFT Magnitude ');     
    zoom on;
end
disp(' ')
disp(' Plot the FFT phase? ')
choice = input(' 1=yes 2=no ');
%
if(choice == 1)
    figure(2);      
    plot(freq,phase)
    xlabel(' Frequency (Hz)');
    ylabel(' Phase (deg)'); 
    title(' FFT Phase '); 
end
%
disp(' ');
disp(' The complex FFT name is  : complex_FFT ' );
disp(' The FFT magnitude name is: magnitude_FFT ' );
%
disp(' ');
disp(' Output full FFT Magnitude file to ASCII text? ');
choice=input(' 1=yes  2=no  ' );
disp(' ')
%
if choice == 1 
        disp(' Enter output filename ');
        g_filename = input(' ','s');
%
        fid = fopen(g_filename,'w');
        for j=1:length(freq);
            fprintf(fid,'%12.5e \t %10.4e \n',freq(j),full(j));
        end
        fclose(fid);
end
%
disp(' ')
%% disp(' Plot the spectrogram? ')
%% disp(' (Requires Signal Processing Toolbox) ')
%% choice = input(' 1=yes 2=no ');
%%
%% if(choice == 1)
%%    figure(3);     
%%    fs=1/dt;
%%    specgram(amp,[],fs)
%%    ylabel(' Frequency (Hz)');
%%    xlabel(' Time (sec)'); 
%% end
end