%
% Function to compute a running periodogram of a uniformly sampled time series
% using FFTs of zero-padded (20x) segments. The output time-frequency 
% running periodogram is plotted as a contour plot
%
% INPUT PARAMETERS
% data--2 columns with independent,dependent variables of time series
% (time,value); time must have increasing values 
% dt--data sample rate (data must be uniformly sampled)
% window--window length for spectrogram (<< total data time) in time units
% step--step increment (<< window) in time units
% fmin--minimum output f (>= 0) in cycles/(time unit)
% fmax--maximum output f (<= Nyquist=1/(2dt)) in cycles/(time unit)
% unit--time unit of data, e.g. 'cm' or 'kyr'
% norm--1=normalize each spectrum; 2=no normalization (1 is recommended)
%
% CALCULATED PARAMETERS:
% kpts = 20 * mpts - padded FFT length of windows (frequency interpolation)
% df = 1/(kpts*dt) - output frequency bin size
% npts = length(xdata)
% nspec = (npts-mpts)/inc +1 - number of output spectra, must be a round number
% kfs = (fmax-fmin)/df+1 - number of points per output spectrum
%
% OUTPUT PARAMETERS:  
%
% s -- 2 columns, each row contains one spectrum (fmin thru fmax by df)
% and contour plot in 'jet' colormap created by pcolor.m
%
% 2002: written by Linda Hinnov and Florian Maurer  
% 2007: update Brian Zurek
% 2013: update by Nic Thibault 
% 2014: title added by Tim Herbert
% 2015: adjustments by Linda Hinnov for Kodama and Hinnov, 2015 book
% 2017: plot axis corrections and normalization options added by Linda Hinnov
%
%TEST:
%
% Load "evoffttest.mat", note the settings in Workspace, and in the
% Command Window enter:
%
% s=evofft(data,window,step,dt,fmin,fmax,unit,norm);
%------------------------------------------------------------------------------------
function s=evofft(data,window,step,dt,fmin,fmax,unit,norm)
x=[ ];
sx=[ ];
s=[ ];
taxis=[ ];
faxis=[ ];
time=data(:,1);
xdata=data(:,2);
mpts=window/dt;
mpts=floor(mpts)
inc=step/dt;
inc=floor(inc);
kpts=20*mpts;
df=1./(kpts*dt);
kfs=(fmax-fmin)/df+1;
kfs=floor(kfs);
kmin=(fmin/df)+1;
kmin=floor(kmin);
kmax=kmin+kfs-1;
npts=length(xdata);
% number of spectra to be calculated
nspec=floor(npts-mpts)/inc+1;
%Start loop through all spectra
for n=1:nspec
m1=inc*(n-1)+1;
m2=m1+mpts-1;
taxis(n)=(m1+(mpts/2)-1)*dt;
		for m=m1:m2
		j=m-m1+1;
		x(j)=xdata(m);
		end
x=detrend(x);
sx=fft(x,kpts);
		for k=kmin:kmax
            k1=k-kmin+1;
            s(n,k1)=sqrt(sx(k)*conj(sx(k)));
		end
end
s=real(s); % prepares spectrogram for plotting
% normalize each spectrum unto itself if norm=1
if norm==1
maxval=max(s');
for i=1:nspec;
    for k=kmin:kmax;
        s(i,k)=s(i,k)/maxval(i);
    end
end
end
% 2-D plot (color contour)
figure;
[nrow,ncol]=size(s);
ddt=time(2)-time(1);
sign=round(ddt/dt);
% Set up x and y axes grids
x_grid=fmin:df:fmax;
% For y axis, start first label at half the window point
igrid=floor(0.5*mpts);
ttt1=time(igrid+1);
%ttt2= ttt1+(nrow-1)*sign*step
ttt2=time(npts-mpts/2-1);
stepit=(ttt2-ttt1)/nrow;
y_grid=ttt1:stepit:ttt2-stepit;
pcolor(x_grid,y_grid,s)
% adjust color and add basic annotations 
colormap(jet)
shading interp 
str=sprintf('Window = %d %s',window,unit);
title(str);
%title(['Time-frequency Landscape', window,' = ', unit])
xlabel(['frequency ( cycles per ',unit,' )']) 
ylabel(['time ( in ',unit,' )'])
colorbar
% EOF
