disp(' ')
disp(' waterfall_FFT.m ')
disp(' ver 2.6  January 21, 2011 ')
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 THM;
clear tmx;
clear tmi;
clear ts;
clear te;
clear n;
clear n1;
clear n2;
clear dt;
clear po;
clear as;
clear ts;
clear amp;
clear tim;
clear aa;
clear bb;
clear Y;
clear store;
%
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
%
size(THM)
tmx=max(THM(:,1));
tmi=min(THM(:,1));
n = max(size(THM(:,2)));
dt=(tmx-tmi)/(n-1);
%
out3 = sprintf('\n  number of samples = %d ',n);
disp(out3)
disp(' ')
disp(' Time Step ');
dtmin=min(diff(THM(:,1)));
dtmax=max(diff(THM(:,1)));
%
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)
%
out3 = sprintf('\n  start  = %g sec    end = %g sec \n',tmi,tmx);
disp(out3)
disp(' Enter processing option for signal duration');
disp(' 1=whole time history   2=segment ');
po=input('');
%
if(po==2)
    disp(' ');
    ts=input(' Enter segment start time  ');
    te=input(' Enter segment end   time  ');
    disp(' ');
%
    if(ts>tmx)
        ts=tmx;
    end
    if(te<tmi)
        ts=tmi;
    end
%
    n1=fix((ts-tmi)/dt);
    n2=fix((te-tmi)/dt);
%
    if(n1<1)
        n1=1;
    end
%
    if(n2>n)
        n2=n;
    end
%
    if(n1>n2)
        n2=n1;
    end
%
    as=THM(:,2);
%
    ts=THM(:,1);
%
    amp=as(n1:n2)';
    tim=ts(n1:n2)';
%
else
    amp=THM(:,2);
    tim=THM(:,1);
end
%
clear THM;
%
for(ijk=1:10)
%
clear Y;
clear store;
clear freq;
clear store_p;
clear freq_p;
clear time_a;
clear dt;
clear df;
clear NW;
clear mmm;
clear n;
clear max_a;
clear max_j;
%
tmx=max(tim);
tmi=min(tim);
%
n = max(size(amp));
dt=(tmx-tmi)/n;
%
disp(' ');
out4 = sprintf(' Time history length = %d ',n);
disp(out4)
%
for(i=1:22)
    if( 2^i > n )
        break;
    end
    N=2^i;
end
%
NC=18;
for(i=4:NC)
    ss(i) = 2^i;
    seg(i) =n/ss(i);
    i_seg(i) = fix(seg(i));
end
%
disp(' ');
out4 = sprintf(' Number of   Samples per   Time per    df');
out5 = sprintf('  Segments     Segment      Segment      ');
%
disp(out4)
disp(out5)
%        
for(i=4:NC)
    if( i_seg(i)>0)
        str = int2str(i_seg(i));
        tseg=dt*ss(i);
        ddf=1./tseg;
        out4 = sprintf(' \t  %s  \t  %d  \t   %6.3f  \t  %6.3f',str,ss(i),tseg,ddf);
        disp(out4)
    end
end
%
disp(' ')
NW = input(' Choose the Number of Segments:  ');
disp(' ')
%
mmm = 2^fix(log(n/NW)/log(2));
%
df=1./(mmm*dt);
%
minf = input(' select minimum output frequency (Hz) ');
disp(' ');
maxf = input(' select maximum output frequency (Hz) ');
%
disp(' ');
disp(' Select Overlap ');
disp(' 1=none  2=50% ');
io=input(' ');
%
for(i=1:mmm/2)
    freq(i)=(i-1)*df;
    if(freq(i)>maxf)
        break;
    end
end
mk=max(size(freq));
t1=tmi+(dt*mmm);
time_a(1)=t1;
%
if(io==1)
    for(i=2:NW)
        time_a(i)=time_a(i-1)+dt*mmm;
    end
else
    NW=2*NW-1;
    dt=dt/2;
    for(i=2:NW)
        time_a(i)=time_a(i-1)+dt*mmm;
    end    
end
%
jk=1;
for(ij=1:NW)
%
   clear sa;
%
   max_a(ij)=0.;
   max_f(ij)=0.;
   if(io==1)   
       for(k=1:mmm)
          sa(k)=amp(jk);
          jk=jk+1;
       end
   else
       for(k=1:mmm)
          sa(k)=amp(jk);
          jk=jk+1;
       end
       jk=jk-mmm/2;
   end
%
   clear Y;
   Y= fft(sa,mmm);
%
   store(ij,1) = abs(Y(1))/mmm; 
%   
   j=1;
   for(k=2:mk) 
      store(ij,k) =2.*abs(Y(k))/mmm;   
      if(freq(k)>=minf)
            store_p(ij,j)=store(ij,k);
            freq_p(j)=freq(k);
            if(store_p(ij,j)>max_a(ij))
                max_a(ij)=store_p(ij,j);
                max_f(ij)=freq_p(j);
            end
            j=j+1;
      end
   end    
end
%
disp(' ')
disp(' Choose color scheme ')
color_choice = input(' 1=default 2=red 3=black ');
%
figure(1);    
if(color_choice ==1)
     colormap(hsv(128));
end  
if(color_choice ==2)
     colormap(hsv(1));
end    
if(color_choice ==3)
     for(i=1:64)
          map(i,:)=[0 0 0];
          colormap(map);
     end      
end 
figure(1);
set(gcf,'renderer','openGL');    
waterfall(freq_p,time_a,store_p);  
xlabel(' Frequency (Hz)');
ylabel(' Time (sec)'); 
zlabel(' Magnitude'); 
view([-12 15]);   
%
disp(' ');
disp(' Plot as surface plot? ')
schoice = input(' 1=yes 2=no ');
%
while(schoice==1)
    disp(' ')
    disp(' Choose Magnitude Format');
    disp(' 1=linear 2=log ');
    imag=input(' ');
    clear store_pp;
    if(imag==1)
        store_pp=store_p;
    else
        store_pp=log10(store_p);
    %
        disp(' ');
        disp(' Select Y-axis format ');
        iyf=input(' 1=default 2=set number of decades ');
        if(iyf==2)
            disp(' ');
            disp(' Select number of decades for Y-axis ');
            idec=input(' ');
            maxP=round(max(max(store_pp))+0.5);
            minP=maxP-idec+0.0001;
            sz=size(store_pp);
            for(i=1:sz(1))
                for(j=1:sz(2))
                    if(store_pp(i,j)<minP)
                        store_pp(i,j)=minP;
                    end
                end
            end
        end
    %
    end
    disp(' ')
    disp(' Selected Shading ');
    ie=input(' 1=flat 2=faceted 3=interp ');
%
    figure(2);
    colormap(hsv(128))
    set(gcf,'renderer','openGL');
    surf(freq_p,time_a,store_pp);
%
    if(ie==1)
        shading flat;
    end
    if(ie==2)
        shading faceted;
    end
    if(ie==3)
        shading interp;
    end
% 
    xlabel(' Frequency (Hz)');
    ylabel(' Time (sec)'); 
    zlabel(' Magnitude');   
%  axis([xmin,xmax,ymin,ymax,zmin,zmax]);
    disp(' ');
    disp(' Replot surface with different options? ');
    disp(' 1=yes 2=no ');
    ire=input(' ');
    if(ire==2)
        break;
    end
end
%
disp(' ');
disp(' Plot the Spectrogram? ')
choice = input(' 1=yes 2=no ');
if(choice==1)
   figure(3);
   colormap(hsv(128));
   surf(freq_p,time_a,store_p,'edgecolor','none')
   colormap(jet); axis tight;
   view(0,90);
   ylabel('Time (sec)'); xlabel('Frequency (Hz)');    
end
%
disp(' ');
disp(' Plot the Contour? ')
c_choice = input(' 1=yes 2=no ');
if(c_choice==1)
    while(1)
        disp(' ');
        disp(' Enter the number of contour levels ');
        n=input(' ');
        figure(4);
        contour(freq_p,time_a,store_p,n);
        ylabel('Time (sec)'); xlabel('Frequency (Hz)');  
        grid on;
        disp(' ');
        disp(' Replot contour with different options? ');
        disp(' 1=yes 2=no ');
        ire=input(' ');
        if(ire==2)
            break;
        end        
    end
end
%
disp(' ')
disp(' Peak Values ');
disp(' Time(sec)  Freq(Hz)  Amplitude ');
for(ij=1:NW)
    out4 = sprintf(' \t  %6.3f  \t  %6.3f  \t  %6.3f',time_a(ij),max_f(ij),max_a(ij));
    disp(out4) 
end
%
%disp(' ')
%disp(' Repeat analysis with new parameters? ')
%repeat = input(' 1=yes  2=no  ');
%
%if(repeat==2)
%    break
%end
break
end
end
if(file_choice==1)
%
%  Rename matlab matrices with descriptive names
%
    signal=[tim amp];
    sz=signal;
    if(sz(2)>sz(1))
        signal=signal';
    end
    out5 = sprintf('\n The time history is renamed as "signal" \n');
    disp(out5);      
end