%[figNum,classification] = ...
%    calc_zeros(figNum,B,F,T,percentile,max_frequency,do_output)
%
%Classifies a time slices in a spectrogram as being solo or non-solo.
%
%Inputs:
%  figNum - first figure number for plots
%  B - the spectrogram
%  F - the frequencies associated with this spectrogram
%  T - the times associated with this spectrogram
%  percentile - the noise cut off.  Pick something between 5 and 10
%  max_frequency - the maximum frequency (in Hz) to look at in the
%    spectrogram
%  do_output - set to true if you want this function to make
%    plots and false otherwise
%
%Outputs:
%  figNum - the next available figure number.  (One more than the
%    last figure number this function used.)
%  classification - the classification of the spectrogram.  The
%    same length as T.  Ones indicate solo music and zeros indicate
%    non-solo music.
%
% Dec-14-2006 ces2130@columbia.edu
function [figNum,classification] = ...
    calc_zeros(figNum,B,F,T,percentile,max_frequency,do_output)

dB = 20*log10(abs(B));

pwr = sum(2*dB,1);

% we don't want to label silence as solo.  Figure out
% approximately where the silence is.
silence_cutoff = prctile(pwr,percentile*3);
silent_times = logical(pwr < silence_cutoff);

if do_output
    %plot what we have
    temp = silence_cutoff*ones(1,length(T));
    figure(figNum);
    figNum = figNum+1;
    plot(T,pwr);
    hold on
    plot(T,temp,'g');
    hold off;
    title('Where is the silence?');
    xlabel('Time (s)');
    legend('dB Power','Silence cutoff (below = silence, above = non-silence)');
end

[nrow, ncol] = size(B);
% if norm == 1
%     % normalize the power
%     dB = dB./meshgrid(pwr,[1:nrow]);    
% end

% figure out where our noise cutoffs are.  In other words, find 
% "zero".
cutoffs = prctile(dB',percentile)';

if do_output
    figure(figNum);
    figNum = figNum+1;
    
    inds = logical(F<=max_frequency);
    
    subplot(211);
    imagesc(T,F(inds),dB(inds,:));
    axis xy
%     if norm == 1
%         title('dB Normalized Spectrgram');
%     else
%         title('dB Spectrogram');
%     end
    title('dB Spectrogram');
    %colorbar;
    xlabel('Time (s)');
    ylabel('Frequency (Hz)');
    
    subplot(212);
    plot(cutoffs(inds),F(inds));
    title('Noise floor cutoffs');
    ylabel('Frequency (Hz)');
    xlabel('Cutoff (dB)');
    
end

cutoffs_repeat = repmat(cutoffs,1,ncol);
small_vals = logical(dB <= cutoffs_repeat);
f_vals = logical(F<=max_frequency);




%now find the proportions for each column
prop = zeros(1,ncol);

for i=1:ncol
    prop(i) = length(find(small_vals(f_vals,i)));
end

% set the no-sound areas to zero so that we don't pick
% up on those
prop(silent_times) = 0;

%do some smoothing on the raw data
smooth_prop = smooth(prop,10);

% find out what passed the test
pass = find(smooth_prop>=mean(smooth_prop));
classification = zeros(length(smooth_prop),1);
classification(pass) = 1;

if do_output
    
    figure(figNum);
    figNum = figNum+1;
    
    subplot(211);
    plot(T,prop);
    title(sprintf('Number of frequencies close to zero below %d Hz',...
        max_frequency));
    axis([min(T) max(T) min(prop) max(prop)]);
    xlabel('Time (s)');
    ylabel('Number of zeros');
    
    subplot(212);
    plot(T,smooth_prop);
    hold on
    plot(T,mean(smooth_prop)*classification,'g');
    hold off
    axis([min(T) max(T) min(smooth_prop) max(smooth_prop)]);
    title('Smoothed zero-frequency count and classification');
    xlabel('Time (s)');
    ylabel('Smoothed zero-frequency count');
    legend('zero-frequency count',...
        'classification (0 = non-solo, non-zero = solo)');
end



