splitlp5.m

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

% FUNCTION: seperate the LP roots into a format and a balanced polynomial.
%   [Fpoly,Bpoly,FF,FB]=splitlp5(oripoly) gets the 13 roots of oripoly and 
%   splits them into Fpoly (constructed with formant roots) and Bpoly 
%   (unrelated roots).
%   In the function, we apply the pole-modificatoin algorithm to adjust
%   the radius of formant roots in order to get an accurate estimation
%
% INPUT: oripoly==LP polynomial, order of 13, the length is 14.
% OUTPUT: Fpoly==formant roots which are conjugated or on half sampling 
%               frequency, order is 10.
%        Bpoly==balanced poles which are real, order is 3.
%        FF == formant frequencies
%        FB == formant bandwidths

function [Fpoly,Bpoly,FF,FB]=splitlp5(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('Error in function splitlp5.');   
   disp('there is an error!  Less than 4 conjugated pairs are found.');
   disp(poles);
   disp(pole1);

end

%% 3. pole modification
%%%%%%%%%%%%%%%%%%%%%%%

% find the deleted poles
for i=1:length(pole1)
    if imag(pole1(i))~=0
       pole1=[pole1 conj(pole1(i))];
    end
end
ang1=angle(pole1);

% adjust the radius of the affected poles
npole=poles; % new poles 

for i=1:5    
    % correction for deleted poles
    tmp=angs(i); % angle of current selected pole
    fac5=ang1-tmp; % angle difference between poles and deleted poles
    idx3=1:length(pole1); % all pole1
    %idx3=find( (-1*pi/6)<fac5 & fac5<pi/6 ); % find pole1 in positive region
    fac5=fac5(idx3);

    dh=1;
    for ii=1:length(fac5)
        rj=abs( pole1(idx3(ii)) );  % radius of deleted pole
        %dh=(1+rj*rj-2*rj*cos( fac5(ii) ) )*2*dh; %=> under correction
        dh=(1+rj*rj-2*rj*cos( fac5(ii) ) )*dh; % => theoretic correction   
    end
 
    % correction for previous adjusted poles
    dhh=1;
    for k=1:(i-1)
        rk=abs( poles(k) ); % radius of previous pole
        delta=angs(k)-tmp; 
        rkp=abs( npole(k) ); % adjusted radius of previous pole
        dhh0=(1+rk*rk-2*rk*cos(delta));
        dhh1=(1+rkp*rkp-2*rkp*cos(delta));
        dhh=dhh*dhh0/dhh1;
    end

    ri=abs(poles(i)); % radius of affetced pole
    rn=1-sqrt(dh*dhh)*(1-ri); % radius of modified pole
    %[rj ri rn]
    npole(i)=rn*exp(j*tmp);
end

% 4. construct the formant polynomial and balanced polynomial
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FF=angle(npole)/pi*5000; % formant frequencies
rad=abs(npole);
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([npole conj(npole)]) );
Bpoly=real( poly(pole1) );

% plot the frequency responses
if nargout<2
   h0=freqz(1,oripoly,1000); % original LP spectrum
   h0=log(abs(h0));
   h1=freqz(1,real(poly([poles conj(poles)])),1000); % formant spectrum
   h1=log(abs(h1));
   h2=freqz(1,Fpoly,1000); % formant spectrum with pole modification
   h2=log(abs(h2));
   plot(1:5:5000,h0,'y-');hold on;
   plot(1:5:5000,h1,'y--');
   plot(1:5:5000,h2,'y.');hold off;
end