view Code/Test/formant_tracker.m @ 4:92ca03a8fa99 tip

Update to ICASSP 2013 benchmark

function[formant]=formant_tracker(input,dimension,step,window_size,p,threshold)
% %Formant_tracker based on LPC Analysis
% format:[formants]=formant_tracker(input,dimension,step,window_size,p,threshold)
% where,
% input:input speech waveform;
% dimension:input size, in number of samples.
% step:step size , in number of samples;
% window_size: frame/window size, in number of samples;
% p: number of LP coeffcients;
% threshold: the threshold value for pole magnitude.Poles having higher magnitude are candidates for being formant
% frequencies.
% 
% Author:Narsimh Kamath.
% NITK, India.Thanks!
num=1;
dimension-window_size;
    
for i=1:step:(dimension/50)-window_size
    
    frame=input(num*step-window_size/2:num*step+window_size/2);
     
 

for k=1:p+1
    autocorrelation(k)=0;
    for l=1:window_size-k
        autocorrelation(k)=autocorrelation(k)+frame(l)*frame(l+k-1);
    end 
    autocorrelation(k)=autocorrelation(k)/window_size;
end

t=window_size/2;
for k=20:t
    acf(k)=0;
    for l=1:window_size-k
        acf(k)=acf(k)+frame(l)*frame(l+k-1);
    end 
    acf(k)=acf(k)/window_size;
end
[rr,ii]=max(acf(10:t));
tp=10000/ii;

alpha=zeros(p,p);
a=zeros(p,1);
E=zeros(p,1);
E(1)=autocorrelation(1);
for i=1:p
    sum=0;
    for j=1:i-1
        if(i)>1
        sum=sum+alpha(j,(i-1))*autocorrelation(abs(i-j)+1);
        end
    end
    
        k(i)=(autocorrelation(i+1)-sum)/E(i);
        alpha(i,i)=k(i);
        if i>1
            for j=1:i-1
                alpha(j,i)=alpha(j,i-1)-k(i)*alpha(i-j,i-1);
            end
        end
        E(i+1)=(1-k(i)^2)*E(i);
        var=E(i+1);
    end
    for j=1:p
        a(j)=alpha(j,p);
    end
    f=[1 -a'];
  gain=0;
  p;
    for i =1:p+1
        gain=gain+f(i)*autocorrelation(i);
    end
   gain;
    gain=gain^0.5;
   

root1=roots([1,-a']);

    mag_root=abs(root1);
    arg_root=angle(root1);
    my_roots(:,num) = mag_root(:);
%     if( max(frame) > 0.001 )
%         subplot(311); plot(frame);
%         subplot(312); plot(abs(fft(frame)));
%         subplot(313); plot(mag_root);
%     end
    k=1;
    for j=1:p
        if mag_root(j)>threshold
            if arg_root(j)>0  &arg_root(j)<pi
            formant(num,k)=arg_root(j)/pi*5000;
            
            num;
            k;
            k=k+1;
            end
        end
    end
   
    
    num=num+1;
end
 s=size(formant);
    for i=1:length(formant)
        for j=1:s(2)
            if formant(i,j)==0
                formant(i,j)=NaN;
            end
        end
    end
   
    plot(formant,'.b');
    return;
end