www.pudn.com > voicebox2.zip > LPCCOVAR.M

function [ar,e,dc]=lpccovar(s,p,t,w) 
%LPCCOVAR performs covariance LPC analysis [AR,E,DC]=(S,P,T) 
% 
%  Inputs:  S(NS)    is the input signal 
%           P        is the order (default: 12) 
%           T(NF,:)  specifies the frames size details: each row specifies one frame 
%                    T can be a cell array if rows have unequal numbers of values 
%                       T(:,1) gives the start of the analysis interval: must be >P 
%                       T(:,2) gives the end of the anaylsis interval [default: t(:+1,1)-1] 
%                       subsequent pairs can be used to specify multiple disjoint segments 
%                    If T is omitted, T(1,1)=P+1, T(1,2)=NS; 
%                    The elements of t need not be integers. 
%           W(NS)    The error at each sample is weighted by W^2 (default: 1) 
% 
% Outputs:  AR(NF,P+1)  are the AR coefficients with AR(:,1) = 1 
%           E(NF)       is the energy in the residual. 
%                       sqrt(E) is often called the 'gain' of the LPC filter. 
%           DC          is the DC component of the signal S. If this output is included, 
%                       the LPC equations are modified to include a DC offset. 
 
% Notes: 
% 
% (1) For speech processing P should be at least 2*F*L/C where F is the sampling 
%     frequency, L the vocal tract length and C the speed of sound. For a typical 
%     male (l=17 cm) this gives f/1000. 
% 
% (2) Each analysis frame should contain at least 2P samples. If note (1) is followed 
%     this implies at least 2 ms of speech signal per frame. 
% 
% (3) It can be advantageous to restrict the analysis regions to time intervals 
%     when the glottis is closed (closed-phase analysis). This can be achieved by 
%     setting the T input parameter appropriately. If the closed-phase is shorter than 
%     2 ms then two or more successive closed-phases should be used by defining 4 or more 
%     elements in the corresponding row of T. 
% 
% (4) A previous version of this routine allowed T() to have a single row which would 
%     be replicated for the entire file length. This has be removed because it gave rise 
%     to an ambiguity. 
 
%  Bugs: sould really detect a singular matrix and reduce the order accordingly 
 
%	Copyright (C) Mike Brookes 1995 
% 
%      Last modified Fri Sep 24 17:51:43 1999 
% 
%   VOICEBOX is a MATLAB toolbox for speech processing. Home page is at 
%   http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   This program is free software; you can redistribute it and/or modify 
%   it under the terms of the GNU General Public License as published by 
%   the Free Software Foundation; either version 2 of the License, or 
%   (at your option) any later version. 
% 
%   This program is distributed in the hope that it will be useful, 
%   but WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
%   GNU General Public License for more details. 
% 
%   You can obtain a copy of the GNU General Public License from 
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to 
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
s = s(:); % make it a column vector 
if nargin < 2 p=12; end; 
if nargin < 3 t=[p+1 length(s)]; end; 
wq = nargin>3; 
[nf,ng]=size(t); 
if iscell(t) 
   t{nf+1}=length(s)+1; 
else 
   if rem(ng,2) 
      t(:,end+1)=[t(2:nf,1)-1; length(s)]; 
   end 
end 
ar=zeros(nf,p+1); 
ar(:,1)=1; 
e=zeros(nf,1); 
dc=zeros(nf,1); 
d0=nargout >2; 
 
rs=(1:p); 
 
for jf=1:nf 
   if iscell(t) 
      tj=t{jf}; 
      if rem(ntj,2) 
         tj(end+1)=t{jf+1}(1)-1; 
      end 
   else 
      tj=t(jf,:); 
   end 
    
   ta = ceil(tj(1)); 
   tb = floor(tj(2)); 
   cs = (ta:tb).'; 
   for js=3:2:length(tj) 
      ta = ceil(tj(js)); 
      tb = floor(tj(js+1)); 
      cs = [cs; (ta:tb).']; 
   end 
   %disp(cs([logical(1); (cs(2:end-1)~=cs(1:end-2)+1)|(cs(2:end-1)~=cs(3:end)-1); logical(1)])'); 
   nc = length(cs); 
   pp=min(p,nc-d0); 
   dm=zeros(nc,pp);	% predefine shape 
   dm(:) = s(cs(:,ones(1,pp))-rs(ones(nc,1),1:pp)); 
   if nargout>2 
      if wq 
         dm = [ones(nc,1) dm].*w(cs(:,ones(1,1+pp))); 
         sc=(s(cs).*w(cs)); 
         aa = (dm\sc).'; 
      else 
         dm = [ones(nc,1) dm]; 
          sc=s(cs); 
        aa = (dm\sc).'; 
      end 
      ar(jf,2:pp+1) = -aa(2:pp+1); 
      e(jf)=sc.'*(sc - dm*aa.'); 
      dc(jf)=aa(1)/sum(ar(jf,:)); 
   else 
      if wq 
         dm = dm.*w(cs(:,ones(1,pp))); 
         aa = (dm\(s(cs).*w(cs))).'; 
      else 
         aa = (dm\s(cs)).'; 
      end; 
      ar(jf,2:pp+1) = -aa; 
      if nargout>1 
         e(jf)=sc.'*(sc - dm*aa.'); 
      end 
   end 
end 
if ~nargout 
   lpcar2db(ar,127); 
end