function onsts = onsets(y, fs, frame)

% find onsets in a piece of sound by looking for derivatives in the
% rows of a spectrogram that many rows share.  Onsets is a matrix
% with two columns, the first is the onsets, the second is the
% offsets.

if(nargin < 3) frame = 512; end
hop = frame/2;

% $$$ % Do something like Fletcher-Munsen
% $$$ wn = [500 4000];   % Stop bands of ear, in Hz
% $$$ [b,a] = butter(2, wn * 2/fs);
% $$$ y = filter(b,a,y);

S = specgram(y, frame, fs);
S = abs(S);

% Normalize spectral bands to be around the same level each
% $$$ S = S.^(1/3);
% $$$ S = diag(1./mean(S')) * S;
% $$$ S = 20*log10(abs(S));
S = sqrt(S);
good = find(isfinite(sum(S)));
Sg = S(:,good);

% Find onsets
[onsts(:,1), D, os, osts] = ons2(Sg, good, hop, fs);

% find offsets
onsts(:,2) = offs(Sg, onsts, good, hop, fs);


% Draw Graphs
subplot 311, imagesc(log(Sg)), axis xy;
subplot 312, imagesc(os), axis xy;
subplot 313, plot([1:size(os,2)], sum(os), [1:size(osts,2)], osts, '*')
  xlim([1 size(D,2)]);
subplot 111


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [o, D, os, osts] = ons3(Sg, good, hop, fs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find onsets in the spectrogram Sg.  All of the elements of Sg are
% finite, and GOOD lists where they came from in an original
% spectrogram with possibly non-finite elements.  HOP is the number
% of samples between STFT frames and FS is the sampling rate of the
% original signal.  The returned vector O is a column of onset
% values in terms of samples in the original time series.

% Longest time window over which to look for onsets is 40ms in each
% direction
time_thresh = .08 / (hop/fs)
% Threshold for filtered contents of a spectral band to be
% considered for being part of an onset
df_thresh = 000 * hop/fs
% threshold for onset as stddevs above the mean filter response.
% Things look pretty gaussian => P_fa = 4.5%
sim_thresh = 2

% Use zero crossings of bandpass filters at different scales to
% find onsets.  Bandpass is Laplacian of Gaussian (c.f. Marr).
i=1;
D = zeros(size(Sg));
while i<time_thresh
  G = 1/sqrt(2*pi*i.^2)*exp(-[-3*i:3*i].^2/(2*i));
  LapG = diff(G,2);
  t1 = filter(LapG, 1, [Sg zeros(size(Sg,1), 3*i-1)]')';
  t1 = t1(:,3*i:end);
    t2 = zeroCrossings(t1, 'neg');
  D = D + t2*i;   % Weight larger scales more
  i = i*2;
end
  
  
os = (D > df_thresh) .* D;
%os = localmax(os', round(time_thresh))';  % local maxes per band
Sos = sum(os);
mulm = mean(Sos);
stdm = std(Sos);
sim_thresh*stdm+mulm
osts(good) = localmax(Sos.*(Sos>sim_thresh*stdm+mulm), round(time_thresh));

% convert back to timescale of original y
o = find(osts) * hop;
o = o(:);  % Make sure o is a column


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [o, D, os, osts] = ons2(Sg, good, hop, fs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find onsets in the spectrogram Sg.  All of the elements of Sg are
% finite, and GOOD lists where they came from in an original
% spectrogram with possibly non-finite elements.  HOP is the number
% of samples between STFT frames and FS is the sampling rate of the
% original signal.  The returned vector O is a column of onset
% values in terms of samples in the original time series.

% Longest time window over which to look for onsets is 40ms in each
% direction
time_thresh = .08 / (hop/fs)
% Threshold for filtered contents of a spectral band to be
% considered for being part of an onset
df_thresh = 0000 * hop/fs
% threshold for onset as a fraction of the mean filter response
sim_thresh = 1.3

% Use ndiff and rectification to extract onsets
i=1;
D = zeros(size(Sg));
while i<time_thresh
  t = ndiff(Sg, i);
  D = D + t .* (t > df_thresh);
  i = i*2;
end
  
  
os = (D > df_thresh) .* D;
%os = localmax(os', round(time_thresh))';  % local maxes per band
Sos = sum(os);
mulm = mean(Sos)
osts(good) = localmax(Sos .* (Sos > sim_thresh*mulm), round(time_thresh));

% convert back to timescale of original y
o = find(osts) * hop;
o = o(:);  % Make sure o is a column


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [o, D, os, osts] = ons(Sg, good, hop, fs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find onsets in the spectrogram Sg.  All of the elements of Sg are
% finite, and GOOD lists where they came from in an original
% spectrogram with possibly non-finite elements.  HOP is the number
% of samples between STFT frames and FS is the sampling rate of the
% original signal.  The returned vector O is a column of onset
% values in terms of samples in the original time series.

% Longest time window over which to look for onsets is 40ms in each
% direction
time_thresh = .08 / (hop/fs)
% Threshold for filtered contents of a spectral band to be
% considered for being part of an onset
df_thresh = 0000 * hop/fs
% threshold for onset as a fraction of the mean filter response
sim_thresh = 1.5

% Bandpass each channel to extract onsets
wn = [3 50];               % Stopbands of bandpass filter in Hz
wn = min(wn*2*hop/fs, .9)  % Convert to frame frequency, limit to <=.9
[b,a] = butter(2, wn);
D = filter(b,a,Sg')';

os = (D > df_thresh) .* D;
%os = localmax(os', round(time_thresh))';  % local maxes per band
Sos = sum(os);
mulm = mean(Sos)
osts(good) = localmax(Sos .* (Sos > sim_thresh*mulm), round(time_thresh));

% convert back to timescale of original y
o = find(osts) * hop;
o = o(:);  % Make sure o is a column


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o = offs(Sg, onsts, good, hop, fs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find offsets in the spectrogram Sg.  Finds one offset for each
% onset in ONSETS, the rest of the arguments are the same as ons(.)
% above.  Returned vector O is a column of offsets in samples of
% the original time series.

% Maximum length of an event
maxlen = .5*fs;

o = min([onsts(2:end); (size(Sg,2)-1)*hop], onsts+maxlen);
o = o(:);  % Make sure o is a column


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = localmax(X, k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Keep local maxima the same, set all other elements of the vector
% X to zero.  K is the window size to look around each element
% (passed to maxfilt).

Xmf = maxfilt(X, k);
y = X .* (X == Xmf);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = orderStat(X, frac)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the frac*size(X,1)-th largest element in each column of X.
% E.g. orderStat(X,1) = max(X), orderStat(X,0) = min(X),
% orderStat(X, .5) = median(X).

k = round(frac*size(X,1));
X = sort(X);
y = X(k,:);



