function F = calcFeatures(y, onsets_t, fs)

% Calculate features for the sound Y at each of the ONSETS provided.
% FS is the sampling rate, used for some feature calculations.  F is a
% matrix with one row per onset and one column per feature.

% Eventually use harmonicity, centroid, f0, noisiness, etc. or shazam
% features.  For now, just unwrap 1/2 second of spectrogram.

frame = 32;
hop = frame/2;
halfSec = round(.5 * fs/hop);

S = abs(specgram(y, frame, fs));
onsets_f = round((onsets_t - 1)/hop + 1);

for i=1:length(onsets_t)
  chunk_t = y(onsets_t(i,1):onsets_t(i,2));
  empty = [];
  chunk_S = S(:,onsets_f(i,1):onsets_f(i,2));
  if(numel(chunk_S) <= 0) continue; end

  F(i,:) = [env(empty, chunk_S, fs) avg_spec(empty, chunk_S, fs)];
  % $$$   F(i,:) = mean(S(:,onsets(i,1):onsets(i,2))');
end

% Normalize features to zero mean, unit var
infinite = find(~isfinite(F));
% $$$ F = F.*F./F;  % Convert inf to NaN
F(infinite) = NaN;
F = F - repmat(nanmean(F), size(F,1), 1);
F = F * diag(1./nanstd(F));


%%%%%%%%%%
% Interface to feature functions: 
%   f = feature(timeform, specgram, fs)
%
% F is a row vector (or scalar), timeform is the waveform, specgram
% is the magnitude of the spectrogram on a linear scale.
%%%%%%%%%%


function f = env(chunk_t, chunk_S, fs)
% Envelope contour at N points, not starting at the beginning, but
% the beginning is normalized to be of magnitude 1.
N = 16;
mag = sum(chunk_S);
% $$$ rs = resample(mag, N+1, length(mag));
% $$$ f = rs(2:end)/rs(1);
% $$$ f = resample(mag, N, length(mag));
if(length(mag) < N)
  f = [mag zeros(1,length(mag)-N)];
elseif(length(mag) > N)
  f = mag(1:N);
end


function f = avg_spec(chunk_t, chunk_S, fs)
% Average spectrum
f = mean(chunk_S, 2)';


function f = entrpy(chunk_t, chunk_S, fs)
% Entropy of the fft of a whole chunk.
f = mean(entropy(chunk_S));


function f = centroid(chunk_t, chunk_S, fs)
% Average spectral centroid across spectrogram
N = (size(chunk_S,1)-1)*2;
freq = [N/2:-1:0]/N * fs;
f = mean(freq * chunk_S);