Dawn@4: function[formant]=formant_tracker(input,dimension,step,window_size,p,threshold)
Dawn@4: % %Formant_tracker based on LPC Analysis
Dawn@4: % format:[formants]=formant_tracker(input,dimension,step,window_size,p,threshold)
Dawn@4: % where,
Dawn@4: % input:input speech waveform;
Dawn@4: % dimension:input size, in number of samples.
Dawn@4: % step:step size , in number of samples;
Dawn@4: % window_size: frame/window size, in number of samples;
Dawn@4: % p: number of LP coeffcients;
Dawn@4: % threshold: the threshold value for pole magnitude.Poles having higher magnitude are candidates for being formant
Dawn@4: % frequencies.
Dawn@4: % 
Dawn@4: % Author:Narsimh Kamath.
Dawn@4: % NITK, India.Thanks!
Dawn@4: num=1;
Dawn@4: dimension-window_size;
Dawn@4:     
Dawn@4: for i=1:step:(dimension/50)-window_size
Dawn@4:     
Dawn@4:     frame=input(num*step-window_size/2:num*step+window_size/2);
Dawn@4:      
Dawn@4:  
Dawn@4: 
Dawn@4: for k=1:p+1
Dawn@4:     autocorrelation(k)=0;
Dawn@4:     for l=1:window_size-k
Dawn@4:         autocorrelation(k)=autocorrelation(k)+frame(l)*frame(l+k-1);
Dawn@4:     end 
Dawn@4:     autocorrelation(k)=autocorrelation(k)/window_size;
Dawn@4: end
Dawn@4: 
Dawn@4: t=window_size/2;
Dawn@4: for k=20:t
Dawn@4:     acf(k)=0;
Dawn@4:     for l=1:window_size-k
Dawn@4:         acf(k)=acf(k)+frame(l)*frame(l+k-1);
Dawn@4:     end 
Dawn@4:     acf(k)=acf(k)/window_size;
Dawn@4: end
Dawn@4: [rr,ii]=max(acf(10:t));
Dawn@4: tp=10000/ii;
Dawn@4: 
Dawn@4: alpha=zeros(p,p);
Dawn@4: a=zeros(p,1);
Dawn@4: E=zeros(p,1);
Dawn@4: E(1)=autocorrelation(1);
Dawn@4: for i=1:p
Dawn@4:     sum=0;
Dawn@4:     for j=1:i-1
Dawn@4:         if(i)>1
Dawn@4:         sum=sum+alpha(j,(i-1))*autocorrelation(abs(i-j)+1);
Dawn@4:         end
Dawn@4:     end
Dawn@4:     
Dawn@4:         k(i)=(autocorrelation(i+1)-sum)/E(i);
Dawn@4:         alpha(i,i)=k(i);
Dawn@4:         if i>1
Dawn@4:             for j=1:i-1
Dawn@4:                 alpha(j,i)=alpha(j,i-1)-k(i)*alpha(i-j,i-1);
Dawn@4:             end
Dawn@4:         end
Dawn@4:         E(i+1)=(1-k(i)^2)*E(i);
Dawn@4:         var=E(i+1);
Dawn@4:     end
Dawn@4:     for j=1:p
Dawn@4:         a(j)=alpha(j,p);
Dawn@4:     end
Dawn@4:     f=[1 -a'];
Dawn@4:   gain=0;
Dawn@4:   p;
Dawn@4:     for i =1:p+1
Dawn@4:         gain=gain+f(i)*autocorrelation(i);
Dawn@4:     end
Dawn@4:    gain;
Dawn@4:     gain=gain^0.5;
Dawn@4:    
Dawn@4: 
Dawn@4: root1=roots([1,-a']);
Dawn@4: 
Dawn@4:     mag_root=abs(root1);
Dawn@4:     arg_root=angle(root1);
Dawn@4:     my_roots(:,num) = mag_root(:);
Dawn@4: %     if( max(frame) > 0.001 )
Dawn@4: %         subplot(311); plot(frame);
Dawn@4: %         subplot(312); plot(abs(fft(frame)));
Dawn@4: %         subplot(313); plot(mag_root);
Dawn@4: %     end
Dawn@4:     k=1;
Dawn@4:     for j=1:p
Dawn@4:         if mag_root(j)>threshold
Dawn@4:             if arg_root(j)>0  &arg_root(j)<pi
Dawn@4:             formant(num,k)=arg_root(j)/pi*5000;
Dawn@4:             
Dawn@4:             num;
Dawn@4:             k;
Dawn@4:             k=k+1;
Dawn@4:             end
Dawn@4:         end
Dawn@4:     end
Dawn@4:    
Dawn@4:     
Dawn@4:     num=num+1;
Dawn@4: end
Dawn@4:  s=size(formant);
Dawn@4:     for i=1:length(formant)
Dawn@4:         for j=1:s(2)
Dawn@4:             if formant(i,j)==0
Dawn@4:                 formant(i,j)=NaN;
Dawn@4:             end
Dawn@4:         end
Dawn@4:     end
Dawn@4:    
Dawn@4:     plot(formant,'.b');
Dawn@4:     return;
Dawn@4: end
Dawn@4: