function onsets = findOnsets(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 = 20*log10(abs(S));
S = abs(S);
good = find(isfinite(sum(S)));
Sg = S(:,good);

% Normalize spectral bands to be around the same level each
Sg = sqrt(Sg);
%Sg = Sg.^(1/3);
%Sg = diag(1./mean(Sg')) * Sg;

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

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


% Draw Graphs
subplot 311, imagesc(Sg), axis xy;
subplot 312, imagesc(os), axis xy;
i = [1:size(D,2)];
subplot 313, plot(i, sum(os), i, 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 = 0000 * hop/fs
% threshold for onset as a fraction of the mean filter response
sim_thresh = 1.3

% Use zero crossings of bandpass filters at different scales to
% find onsets
i=1;
D = zeros(size(Sg));
wn = [3 6];  % Stopbands of first bandpass filter in Hz
wn = wn*2*hop/fs  % Convert to frame frequency
while i<time_thresh
  % Bandpass each channel to extract onsets
  [b,a] = butter(4, wn);
  t1 = filter(b,a,Sg')';
  t2 = zeroCrossings(t1, 'pos');
  subplot 211, imagesc(t1)
  subplot 212, imagesc(t2), pause(0)
  D = D + t2;
  i = i*2;
  wn = min(wn*2, .9)  % Take next octave, limit high edge to .9
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] = 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, onsets, 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([onsets(2:end); size(Sg,2)*hop], onsets+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,:);