splitlp3.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 83 行

% FUNCTION: seperate the LP roots into a format and a balanced polynomial.
%   [Fpoly,Bpoly,FF,FB]=splitlp3(oripoly) gets the 12 roots of oripoly and 
%   splits them into Fpoly (conjugated) and Bpoly (real roots).
%
% INPUT: oripoly==LP polynomial, order of 12, the length is 13.
% OUTPUT: Fpoly==formant roots which are conjugated or on half sampling 
%               frequency, flexible number.
%        Bpoly==balanced poles which are real, flexible number.
%        FF == formant frequencies
%        FB == formant bandwidths

function [Fpoly,Bpoly,FF,FB]=splitlp3(oripoly);

%%% 0. calculate the roots(poles) from the LP polynomials
%%%    and sort it according to their angles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
poles=roots(oripoly);
poles=poles(:)';
poles=poles( imag(poles)>=0 ); % remove the conjugate part
angs=angle(poles);
if find(angs==-pi)~=[]
   xindx=find(angs==-pi);
   angs(xindx)=-1*angs(xindx);
end
[angs,a_sort]=sort(angs);
poles=poles(a_sort);

%%% 1. put poles with imagnary part==0 into pole1; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
indx1=find( imag(poles)==0 );
pole1=poles(indx1);
poles(indx1)=[];
angs(indx1)=[];

%% 2. Check if the number of poles in poles
%%    >5  => find the biggest bandwidth to frequency ratio
%%    <5  => issue an error message
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
npole=length(poles);

if npole>5
   ff=angs/pi*5000; % formant frquency
   rad=abs(poles); % pole radius
   tmp=(4*rad-1-rad.^2)./(2*rad);
   fb=acos(tmp)/pi*10000; % formant bandwidth

   ratio=fb./ff; % bandwidth to frequency ratio
   [dum,indx2]=sort(ratio);
   unwnt=indx2(6:npole); % the extra pole with big ratio
   pole1=[pole1 poles(unwnt)];
   poles(unwnt)=[];
   angs(unwnt)=[];

elseif npole==4
   [dum1,indx2]=min(pole1);
   pole1(indx2)=[];
   [dum2,indx2]=min(pole1);
   pole1(indx2)=[];
   sudo=sqrt((1+dum1)*(1+dum2))-1;
   sudo=sudo*exp(j*3.11); % create a complex pair of roots
   poles=[poles sudo];
   angs(5)=3.11;

elseif npole<4   
   disp('there is an error!  Less than 4 conjugated pairs are found');
   disp(poles);
   disp(pole1);

end

% 3. construct the formant polynomial and balanced polynomial
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FF=angle(poles)/pi*5000;
rad=abs(poles);
tmp=(4*rad-1-rad.^2)./(2*rad);
FB=acos(tmp)/pi*10000; % formant bandwidth
[FF,idx]=sort(FF);
FB=FB(idx);

Fpoly=real( poly([poles conj(poles)]) );
Bpoly=real( poly(pole1) );