rt2lp.m

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

% FUNCTION: convert formant frequencies and bandwidths to a formant
%           polynomial.
%  fpoly=rt2lp(nFF,FF,FB) generate an filter polynomial from formant frequencies
%  and bandwidths.  In the function, we apply the pole-modification algorithm 
%  to adjust the radius of the formant roots in order to obtain a closer spectrum.
%
% INPUT: nFF == new formant frequencies
%        FF == original formant frequencies
%        FB == original formant bandwidths
% OUTPUT: fpoly == formant polynomial with modification
%         opoly == formant polynomial without modification
%         nFB == adjusted formant bandwidth
% See Also: splitlp5.m

function [fpoly,opoly,nFB]=rt2lp(nFF,FF,FB);

% 0. convert frequencies and bandwidths into roots
%------------------------------------------------

   numff=length(FF);
   [FF,indx]=sort(FF);
   FB=FB(indx);

   % here we assume the sampling frequency is 10k Hz
   phi=FF*pi/5000;  % original angle of the root
   tmp=cos( FB*pi/10000 );
   for k=1:numff
        tr=roots([1 2*tmp(k)-4 1]);
        tr=tr( tr<1 & tr>0 );
        if length(tr)~=1
           tr=exp(-1.0*FB(k)*pi/10000);
        end
        rdis(k)=tr;  % original radius of the root   
        frt(k)=tr*exp(j*phi(k));  % original formant root
   end

   % original formant polynomial before modification
   opoly=poly(frt); % poly without modification

% 1. estimate the spectral density for the new formant frequencies
%-----------------------------------------------------------------

   h0=freqz(1,opoly,phi);
   h0=log(abs(h0));

   h0=sort(h0);  %%==> create spectral tilt
   h0=h0(numff:-1:1);

   theta=nFF*pi/5000;  % new angle of the root

   x0=2*phi(1)-phi(2);
   x6=2*phi(numff)-phi(numff-1);
   h0_0=2*h0(1)-h0(2);
   h0_6=2*h0(numff)-h0(numff-1);
   tmp_h1=interp1([x0 phi x6],[h0_0 h0 h0_6],theta(2:numff-1) );

   h1(1)=h0(1);
   h1(numff)=h0(numff);
   h1(2:numff-1)=tmp_h1;
   h1=exp(h1);

   %for k=1:numff
   %    if theta(k)==phi(k)
   %        h1(k)=exp(h0(k));
   %    end
   %end

% 2. adjust the radius
%---------------------

Nrdis=rdis;
npole=zeros(1,numff);

for loop=1:2
    for k=1:numff

        dh=h1(k)*h1(k); % original spectral response
        Ndh=1;  % new spectral response

        % effect of other roots
        for ii=1:numff
           if ii~=k
              rj=Nrdis(ii);
              dthe=theta(ii)-theta(k);
              Ndh=(1+rj*rj-2*rj*cos(dthe) )*Ndh;
           end
        end

        rn=1-1/sqrt(dh*Ndh);

        if rn<0.5
          rn=max(0.5,rn*1.2);
        elseif rn>1
          rn=1-(rn-1);
        end

        Nrdis(k)=rn;  % new formant radius    
        npole(k)=rn*exp( j*theta(k) ); % new formant root
        tmp=(4*rn-1-rn*rn)/(2*rn);
        nFB(k)=acos(tmp)/pi*10000;  % new formant bandwidth
     end
end

% 2. construct the formant polynomial
%------------------------------------

opoly=real( poly([frt conj(frt)]) ); % poly without modification
fpoly=real( poly([npole conj(npole)]) ); % poly with modification