function [pitchtrack]=pitch(spec,samprate)
% [pitchtrack]=pitch(spec,samprate)
% 
% PITCH input parameters: spec (MxN matrix returned from specgram()
%                         samprate (scalar holding sample rate)
%
% PITCH output:           pitchtrack (fundamental freq vs. time vector)
%
% PITCH description:      PITCH takes in a specgram and corresponding sample
%                         rate to create a pitch track, a plot of the
%                         fundamental frequency vs. time.  It then eliminates
%                         false pitch estimates by zeroing the areas of low
%                         energy (because no one is talking or a fricative
%                         occurs).
%
% ELEC 301 Project - "Mars Lander the Theologian", aka, "Speaker Verification"
% Sara MacAlpine, Aamir Virani, Nipul Bharani, JP Slavinsky
% WRC (OCEE) Class of 2001 - Fall '99

b           = abs(spec);                     %need frequency magnitudes
bsqr        = b.^2;                          %square magnitudes
energy      = sum(bsqr);                     %Parseval: energy vs time
[numf,numt] = size(b);                       %get num of div in time and freq
f1200       = ceil(1200*numf/(samprate/2));  %row that holds 1200Hz value
b           = b((1:f1200),:);                %fund freq lies b/w 0-1200Hz

%get list of three largest peaks at each time
dmat=diff(b);                    %get diff b/w freq
[numFdiff,numTdiff]=size(dmat);  %get size of diff matrix
for i=1:numTdiff,                %go thru time columns
   threes = zeros(3,2);          %a new vector for each time t col1=n col2=mag
   for n=2:numFdiff,             %go thru freq diffs
      if ((dmat(n,i)<0) & (dmat(n-1,i)>0)) %if slope before = + and after = -
         if (b(n,i) > min(threes(:,2)))    %if peak mag is one of three largest
            if (threes(1,2) == min(threes(:,2))) %then record in this part
               threes(1,1) = n;
               threes(1,2) = b(n,i);
            elseif (threes(2,2) == min(threes(:,2)))
               threes(2,1) = n;
	       threes(2,2) = b(n,i);
            else
	       threes(3,1) = n;
               threes(3,2) = b(n,i);
            end
         end
      end
   end
   freqlist(:,i) = threes(:,1); %gives column three largest vals at time i
end
 
%create initial pitch track
for i=1:numTdiff,                              %go thru times
    sorted        = sort(freqlist(:,i));       %sort three n at each time
    dvals         = diff(sorted);              %to get d1 and d2
    n             = dvals(1)/dvals(2);         %n for first pitch estimate
    F0            = (dvals(1)+dvals(2))/(n+1); %pitch estimate
    pitchtrack(i) = (F0/numf)*(samprate/2);    %return actual freq values
end

%zero points in peak track with low energy
noise      = energy(1:10);
avg        = mean(noise);
dev        = std(noise);
threshold  = avg+200*dev;
killzone   = energy>threshold;
pitchtrack = pitchtrack.*killzone;