function [instfreq,t]=installen(x,t,wlen,olap)
% [instfreq,t]=INSTALLEN(x,t,wlen,olap)
%
% Computes the predominant frequency of a signal sampled at regular
% intervals. Following a method suggested by Allen and Kanamori (2003). 
%
% INPUT:
%
% x           The real signal
% t           The times at which the signal is sampled, in s
% wlen        Window length, in seconds (default: 3)
% olap        Window overlap, in percent (default: 95)
%
% OUTPUT:
%
% instfreq    The instantaneous (not angular) frequency
% t           The times at which this is evaluated
%
% EXAMPLE: (Careful with Nyquist frequency and window length!)
%
% t=linspace(0,8*pi,500); f=sin(2*pi*3*t);
% subplot(211); plot(t,f);
% subplot(212); [in,t]=installen(f,t); plot(t,in);
% ylim([0 5])
%
% Last modified by fjsimons-at-alum.mit.edu, 14.04.2006

defval('wlen',3)
defval('olap',95)

% Computes sampling interval
dt=t(2)-t(1);

% Window length, in samples
wlen=round(wlen/dt);

% Overlap, in samples
olap=floor(olap/100*wlen);

% It's a boxcar window in this case
dwin=ones(wlen,1);

% Determine window parameters and truncation of data set
% Number of points in the signal
npts=length(x);
% Effective length of the window
eflen=wlen-olap;
% Check this works
checkit=(npts-olap)/eflen;
nwin=floor(checkit);
disp(sprintf('Number of overlapping data segments: %i',nwin))   
if nwin~=checkit
  disp(sprintf(...
      'Number of  segments is not  integer: %i / %i points truncated',...
      npts-(nwin*eflen+olap),npts))
end
if nwin<=0; error('Data sequence not long enough'); end
% First, divide the signal in overlapping windowed segments
rows=[1:wlen];
cols=[0:(nwin-1)]*eflen;
xsd=detrend(...
    x(repmat(rows(:),1,nwin)+...
      repmat(cols,wlen,1)));
% This is the overlapping windowed data set
% Do not normalize this window... this is not the spectrogram
xsdw=xsd.*repmat(dwin,1,nwin);

% Now do algorithm from Allen and Kanamori
% Don't do the cumlative sum, just return the result at the end of the  
% window. Note that this is unstable due to the division by a small number!
in=sum([repmat(0,1,nwin) ; diff(xsdw,1)/dt].^2,1);
xsdw=sum(xsdw.^2,1);
in=sqrt(in(:)./xsdw(:));

instfreq=in/2/pi;
t=cols(:)*dt;