clear all
%step 1: pitch of human is located between 50 Hz and 500 hz so filter all frequencies below 50 hz and above 600 hz
%step 2: Find a vowel in the sentence because only vowel will contain the pitch
%step 3: Vowel can be determined by its high energy
%step 4: divide the sentence to segments of size 50 msec
%step 5: calculate the energy for each segment from its spectral components
%step 6: choose the segment with the highest energy
%step 7: Find the FFT and count cepstrum of the vowel segment
%step 8:
[sigOrg, fs] = wavread('Aaron1.wav');
[bl,al]=butter(4,[50 500]/(fs/2)); %creating the coefficients for a bandpass filter with low freq=50 and high freq = 600 Hz
sig=filter(bl,al,sigOrg); %filtering out frequencies below 50 and above 600 Hz from the signal
sigLngth= length(sig);
WindSize = round((20/1000)*fs);     %size of the window is 50 msec to look for the vowel
NumSegm = floor(sigLngth/WindSize); %Number of segments in the entire sentence
for i=1:NumSegm                   %This FOR loop will calculate the energy in every segment
    windStrt = (i-1)*WindSize+1;
    windEnd = windStrt+WindSize;
    EnrgySegm(i)=sum(abs(sig(windStrt:windEnd)));
end
[maxSpec, indx]=max(EnrgySegm);  %segment with highest energy represent a vowel
% WindSize=WindSize*3;

x = sig((indx-1)*WindSize:(indx+1)*WindSize);
N = length(x);
Frs = fs/N;
x = x(:); % assure column vector
plot(x);
figure
%To find the pitch in
%the time domain you need to find the period T between
%two consecutive peaks in the autocorrelation of the vowel segment and the pitch will be 1/T
sigCorr = xcorr(x);
plot(sigCorr)

%To find the pitch in the frequency domain you need to find the frequency F0 between
%two consecutive peaks in the FFT of the vowel segment and the pitch will
%be F0=(number of samples between consecutive peaks)*Frs


window =  hamming(N);
% Find the fourier transform of a the vowel segment
x = x(:).* window(:);
y = fft(x, N);
yabs = abs(y);
L = length(yabs);
ypos = yabs(1:round(500/Frs));
figure
plot(ypos)

% [A,K] = findpeaks(x);
[A,K] = findpeaks(ypos);
Umax = max(K);
Umin = min(K);
Sd = Umax - Umin;
Fo = Frs*Sd/(length(K)-1)