fbana.m

这是一个用于语音信号处理的工具箱

   [Froot,FF,FB]=smfrmt1(cofa);

end %% if Fmethod==1

    %%%-----------------------------%
    %%% smoothing the formant track %
    %%%-----------------------------%

 % find error frame
 for kf=1:nframe-1
     rcd(kf)=0; % rcd==0 --> no error ; rcd==-1 --> error
     if voicetype(kf)

        ff=FF(kf,:);
        fb=FB(kf,:);

        % the formant bandwidth is gibber than frequency
        if ~isempty(find( (ff-fb)<0 ))
           rcd(kf)=-1;
        % the second and third bandwidth to formant ratio should not exceed .5
        elseif fb(2)/ff(2)>.5 & fb(3)/ff(3)>.5
           rcd(kf)=-1;
        elseif ff(5)>4600
           rcd(kf)=-1;
        end
     end
 end

 %correct error frame
 for kf=1:nframe-1
    if voicetype(kf) & rcd(kf)==-1
       back=kf-1;
       while rcd(back)==-1
            back=back-1;
       end
       forward=kf+1;
       while rcd(forward)==-1
            forward=forward+1;
       end
       w1=kf-back;
       w2=forward-kf;
       cofb=real(poly(Froot(back,:)));
       coff=real(poly(Froot(forward,:)));
       cofa1=polystab(w1/(w1+w2)*cofb+w2/(w1+w2)*coff);

       rr=roots(cofa1);
       rr=rr(:)';
       rr=rr(imag(rr)>0);
       ff=angle(rr)/pi*5000;
       rad=abs(rr); % pole radius
       tmp=(4*rad-1-rad.^2)./(2*rad);
       fb=acos(tmp)/pi*10000; % formant bandwidth
       [ff,idx]=sort(ff);
       fb=fb(idx);

       if ~isempty(find( (ff-fb)<0 ))
            %zplane(cofa1,[]);pause(1);
            disp('Problem is found in the following frame!!');kf
            Froot(kf,:)=Froot(back,:);
            FF(kf,:)=FF(back,:);
            FB(kf,:)=FB(back,:);
       else
            Froot(kf,:)=[rr conj(rr)];
            FF(kf,:)=ff;
            FB(kf,:)=fb;
       end
    end
 end

     %%-------------------%%
     %% inverse filtering %%
     %%-------------------%%

froot=zeros(10,1);
kf=1;
      froot=Froot(kf,:)';
      fmpoly=real(poly(froot));
      order=length(fmpoly)-1;
 
      %****** frequency lifting ******%
      %cofa1=cofa(kf,:);
      %fmpoly=ad_lpc1(cofa1,fmpoly);

      coffa(kf,:)=fmpoly;

      sso=signal( (kf-1)*M_len+1:kf*M_len+order+O_lap );
      dgf1=filter(fmpoly, 1, sso);
      dgf1(1:order)=[];
      sso(1:10)=[];

      % Normalize th gain
      eng=sso*sso';
      eng1=dgf1*dgf1';
      dgf1=sqrt(eng/eng1)*dgf1;

      dgf( (kf-1)*M_len+order+1:kf*M_len+O_lap+order)=dgf1;
      dgf2=dgf1(M_len+1:F_len);

for kf=2:nframe-1
        froot=Froot(kf,:)';
        fmpoly=real(poly(froot));
        order=length(fmpoly)-1;

        %******* frequency lifting ******%
        %cofa1=cofa(kf,:);
        %fmpoly=ad_lpc1(cofa1,fmpoly);

        coffa(kf,:)=fmpoly;
        
        % inverse filtering
        sso=signal( (kf-1)*M_len+1:kf*M_len+O_lap+order );
        dgf1=filter(fmpoly, 1, sso);
        dgf1(1:order)=[];
        sso(1:10)=[];

        % Normalize th gain
        eng=sso*sso';
        eng1=dgf1*dgf1';
        dgf1=sqrt(eng/eng1)*dgf1;

        % overlap region
        dgfo=( dgf2.*(O_lap:-1:1)+dgf1(1:O_lap).*(1:O_lap) )/(O_lap+1); 
      dgf((kf-1)*M_len+order+1:kf*M_len+O_lap+order)=[dgfo dgf1(O_lap+1:F_len)];
        dgf2=dgf1(M_len+1:F_len);
end

%disp('Inverse filtering is complete!');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                           %%%%%%
%%%%%%%  Find glottal coefficient %%%%%%
%%%%%%%  6-order polynomial model %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % This step will create the following variables :
 %
 %	gpcof = glottal codeword for each pitch period.

gpcof=zeros(nframe,7);

dgf=filter([1 1],[1 0.98],dgf); % remove high-frequency noise 
%dgf=integ(dgf);  % integrate the differential glottal flow in a whole 
%gf=filtfilt([1 -1],[1 -0.99],dgf); % gf== glottal flow

          %***** Select Glottal Source model
          %***** Smodel==1 --> 6 order polynomial model  sounds better ******%
          %***** Smodel==2 --> LF model                                ******%

Smodel=basic(3);

for kf=2:nframe-1
	if ploc1(kf,1)
	    startp=M_len*(kf-1)+Order;
	    self=ploc1(kf,:);
	    self=self(self~=0);
	    if ploc1(kf+1,1)
		after=ploc1(kf+1,1)+M_len;
	    elseif voicetype(kf+1) & kf <= nframe-2 
	      	if ploc1(kf+2,1)
		    after=ploc1(kf+2,1)+2*M_len;
		else
		    after=[];
	        end;
	    else
		after=[];
	    end;
	    period=[self after]+startp;
	    lenp=length(period);
	    if lenp>1
		gcof=[];
		for k=1:lenp-1
                        
			seggp=dgf(period(k):period(k+1));
                        seggp=seggp-mean(seggp);

                        % Polynomial model when Smodel==1 ***
                        if Smodel==1
			   gcof(k,:)=polym1( seggp );
                           %  poly1 <--> gen_dgf1

                        % LF model when Smodel==2 ***
                        elseif Smodel==2
                           seggp=integ(seggp);
                           [Rd,Ee]=lfmodel(seggp);
                           Ra=(-1+4.8*Rd)/100;
                           Rk=(22.4+11.8*Rd)/100;
                           Rg=Rk/4/( 0.11*Rd/(0.5+1.2*Rk)-Ra);
                           T0=length(seggp)-1;
                           Ta=T0*Ra;
                           Tp=T0/2/Rg;
                           Te=floor(Tp*Rk+Tp);
                           Tc=T0;
                           gcof(k,1:7)=[Tp/T0 Te/T0 Ta/T0 Tc/T0 Ee 0 0];
                        end
		end;
		if lenp<=2
			avegcof=gcof;
		else
			avegcof=mean(gcof);
		end;
		gpcof(kf,:)=avegcof;

            end % if lenp>1
	end;
end;
gcof=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  energy encoding for voiced speech  %%
%%%%%%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This step will create the following variables :
%
%	gm = signal power for each pitch period.

% encoding the energy for pitch period
numgci=length(gci);
gm=zeros(1,numgci);
if numgci>1
   lenss=gci(2)-gci(1);
end
for njk=1:numgci

       startp=gci(njk)+1;  % starting point of the pitch period
       kf=fix( (startp-Order)/M_len+1 ); % frame number
       if njk==numgci
          endp=startp+lenss-1;
       else
          endp=gci(njk+1);  % ending point of the pitch period
          olenx=lenss;
          lenss=gci(njk+1)-gci(njk);

          if lenss>1.5*olenx
             if voicetype(kf)==1 & voicetype(kf+1)==0
                lenss=olenx;
                endp=startp+lenss-1;
             end
          end
       end

       sp=signal(startp:endp);
       gm(njk)=(sp*sp')/lenss;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                                        %%
%%%%%%%  energy encoding for the whole speech  %%
%%%%%%%                                        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This step will create the following variables :
%
%    ncgm = signal power for Xlen period; size(ncgm)=[nframe,4]

ncgm=zeros(nframe,4);
Xlen=M_len/4;
for kf=1:nframe
    for kk=1:4
        startp=(kf-1)*M_len+Order+(kk-1)*Xlen+1;
        ss1=signal(startp:startp+Xlen-1);
        ncgm(kf,kk)=(ss1*ss1')/Xlen;
    end
end

%disp('Glottal codeword searching is completed!'); toc;
%disp('');
%disp('It takes a while to perform codeword searching for unvoiced speech.');
%disp('Please wait! Please wait! Please wai--------------------------------t.');
%disp('');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                            %%
%%%%%%%  Find unvoiced codeword    %%
%%%%%%%                            %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % This step will create the following variables :
 %
 %	ncgm = segmental power for one pitch period.
 %	ncidx = glottal codeword.

% Function : preprocess before performing codeword searching.  This includes
%		1. removal of filter memory contribution,
%		2. spectral weighting.

load nc;  % Load the stochastic codebook
Xlen=M_len/4;
Ziy1=zeros(1,Order);  % Reset filter memory at the beginning
for kf=1:nframe
	if kf==1 | kf==nframe
		proc='nc';
	elseif sum(voicetype(kf-1:kf))<2 | sum(voicetype(kf:kf+1))<2
		proc='nc';
	else
		proc='mp';
	end;
	if proc=='nc'  % Proceed if it is an unvoiced segment
	    ocofa=cofa(kf,:);
            ncofa=ocofa.*((0.8).^(0:length(ocofa)-1));
	    if kf>1
	    	if voicetype(kf-1)
			Ziy1=signal((kf-1)*M_len+Order:-1:(kf-1)*M_len+1);
			% Note : filter memory is regarded as the past samples
			% for the Direct-I implementation
		end;
	    end;
	    if kf==1 
		so=signal((kf-1)*M_len+Order+1:kf*M_len+Order);
	    else
		so=signal((kf-1)*M_len+Order+1:kf*M_len+Order)-...
		    filt(1,ocofa,zeros(1,M_len),Ziy1);
		% Remove filter memory carried over from the previous frame
	    end
	    ss=filter(ocofa,ncofa,so);  % spectral weighting 
	    % PURPOSE: De-emphasize errors at the formant regions
	    noise1=[];
	    Ziy=Ziy1;  % filter memory

	    for kk=1:4
		cmp1=ss(1+(kk-1)*Xlen:Xlen*kk);
		sub1=filt(1,ncofa,zeros(1,Xlen),Ziy);
		cmp=cmp1-sub1;  % eliminate memory contribution
		ncidx(kf,kk)=fdncb0(cmp,ncofa,nc0,Rnc0);

                % track the filter memory
                eng1=ncgm(kf,kk);
		noise=nc0(ncidx(kf,kk),:);
                if Xlen~=50
                   noise=interpft(noise,Xlen);
                end
                ss2=filt(1,ocofa,noise,Ziy);
                eng2=sum(ss2.^2);
                amp=sqrt(eng1/eng2);

		% Concatenate the excitation in every subframe as a whole
		noise1=[noise1 amp*noise];
	    	nsyn=amp*ss2;
		Ziy=nsyn(Xlen:-1:Xlen-Order+1);
	    end
	    ss1=filt(1,ocofa,noise1,Ziy1);
	    Ziy1=ss1(M_len:-1:M_len-Order+1);  % Track the filter memory
	end	
end

clear Rnc0 nc0 ss1 ss2 noise noise1;

%disp('Analysis is OK! ');toc;