clear;clc;close all;

[input,fs] = wavread('child_O.wav');
fn = fs/2;

%pre emphasis filter;
% num=[1 -.87];
% den=[1];
% input = filter(num,den,input);

%assuming males have 1 formant per 1kHz, females 1 per 1.5kHz and children
%1 per 2kHz
%p = ceil((fn)/1500)*2 + 2 + 4;
p=6;
%create a 25ms window
win = hamming(200);
sn = input(1:200) .* win;
autocorr = xcorr(sn);
rn = autocorr(200:399);
rn = rn/max(rn);

num=[1 -rn(2)/rn(1)];
den=[1];
sn = filter(num,den,sn);

autocorr = xcorr(sn);
rn = autocorr(200:399);
rn = rn/max(rn);

[matlabAlpha matlaberror matlabK] = levinson (rn, p);
matlabAlpha=matlabAlpha';

%first iteration of algorithm for i=1

Energy(1) = rn(1);
k(1) = rn(2)/Energy(1);
alpha(1,1) = k(1);

%loop for rest of iterations

for i=2:p
    summation = 0;
    for j=1:i-1
    summation = alpha(j,i-1)*rn((i-j)+1)+summation;    
    end   
    
Energy(i) = (1-k(i-1)^2) * Energy(i-1);
k(i) = (rn(i+1) - summation)/Energy(i);   
   
alpha(i,i) = k(i);

for j2 = 1:i-1
   alpha(j2,i) = alpha(j2,i-1)-k(i)*alpha(i-j2,i-1);
end

end

lpcoeffs = [1; -alpha(:,length(alpha))];
PARCOR = k;

alphavalues = alpha(:,length(alpha));
summation = 0;
for w = 1:length(alpha)
    
    summation = alphavalues(w)*rn(w+1);
        
end

gain = sqrt(rn(1)-summation)


[H W] = freqz(gain,lpcoeffs);
W = (W/pi)*fs/2;
plot(W,db(abs(H)))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Linear Prediction Spectral Envelope')

%get normalized cross sectional area values of vocal tract
A(length(PARCOR)) = 1;

for n=length(PARCOR):-1:2
   
    A(n-1) = A(n)*(1+PARCOR(n-1))/(1-PARCOR(n-1));
    
end

N_length_vector = [0 .019 .038 .057 .076 .095];
figure,stairs(N_length_vector,A/2)
hold on
stairs(N_length_vector,-A/2)
title('Estimated Vocal Tract Shape for 6 Tubes')
xlabel('Length (m)')
ylabel('Normalized Area')




