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;