lpcauto

LPCAUTO  performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)
  Inputs:
     s(ns)   is the input signal
       p       is the order (default: 12)
       t(nf,3) specifies the frames size details: each row specifies
               up to three values per frame: [len anal skip] where:
                   len     is the length of the frame (default: length(s))
                   anal    is the analysis length (default: len)
                   skip    is the number of samples to skip at the beginning (default: 0)
               If t contains only one row, it will be used repeatedly
               until there are no more samples left in s.

 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 filter.
          k(nf,2)    gives the first and last sample of the analysis interval

0001 function [ar,e,k]=lpcauto(s,p,t)
0002 %LPCAUTO  performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)
0003 %  Inputs:
0004 %     s(ns)   is the input signal
0005 %       p       is the order (default: 12)
0006 %       t(nf,3) specifies the frames size details: each row specifies
0007 %               up to three values per frame: [len anal skip] where:
0008 %                   len     is the length of the frame (default: length(s))
0009 %                   anal    is the analysis length (default: len)
0010 %                   skip    is the number of samples to skip at the beginning (default: 0)
0011 %               If t contains only one row, it will be used repeatedly
0012 %               until there are no more samples left in s.
0013 %
0014 % Outputs:
0015 %          ar(nf,p+1) are the AR coefficients with ar(1) = 1
0016 %          e(nf)      is the energy in the residual.
0017 %                     sqrt(e) is often called the 'gain' of the filter.
0018 %          k(nf,2)    gives the first and last sample of the analysis interval
0019 
0020 % Notes:
0021 %
0022 % (1) The first frame always starts at sample s(1) and the analysis window starts at s(t(1,3)+1).
0023 % (2) The elements of t need not be integers.
0024 % (3) The analysis interval is always multiplied by a hamming window
0025 % (4) As an example, if p=3 and t=[10 5 2], then the illustration below shows
0026 %     successive frames labelled a, b, c, ... with capitals for the
0027 %     analysis regions. Note that the first frame starts at s(1)
0028 %
0029 %      a a A A A A A a a a b b B B B B B b b b c c C C C C C c c c d ...
0030 %
0031 %     For speech processing, it can be advantageous to restrict the analysis regions
0032 %     to time intervals when the glottis is closed.
0033 %
0034 % (5) Frames can overlap: e.g. t=[10 20] will use analysis frames of
0035 %     length 20 overlapped by 10 samples.
0036 % (6) For speech processing p should be at least 2*f*l/c where f is the sampling
0037 %     frequency, l the vocal tract length and c the speed of sound. For a typical
0038 %     male (l=17 cm) this gives f/1000.
0039 
0040 %      Copyright (C) Mike Brookes 1997
0041 %      Version: $Id: lpcauto.m 713 2011-10-16 14:45:43Z dmb $
0042 %
0043 %   VOICEBOX is a MATLAB toolbox for speech processing.
0044 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0045 %
0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0047 %   This program is free software; you can redistribute it and/or modify
0048 %   it under the terms of the GNU General Public License as published by
0049 %   the Free Software Foundation; either version 2 of the License, or
0050 %   (at your option) any later version.
0051 %
0052 %   This program is distributed in the hope that it will be useful,
0053 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0054 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0055 %   GNU General Public License for more details.
0056 %
0057 %   You can obtain a copy of the GNU General Public License from
0058 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0059 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0061 
0062 s = s(:); % make it a column vector
0063 if nargin < 2 p=12; end;
0064 if nargin < 3 t=length(s); end;
0065 %if nargin < 4 w='ham'; end;
0066 [nf,ng]=size(t);
0067 if ng<2 t=[t t]; end;
0068 if ng<3 t=[t zeros(nf,1)]; end;
0069 if nf==1
0070     nf=floor(1+(length(s)-t(2)-t(3))/t(1));
0071     tr=0;
0072 else
0073     tr=1;
0074 end;
0075 
0076 ar=zeros(nf,p+1);
0077 ar(:,1)=1;
0078 e=zeros(nf,1);
0079 
0080 t1=1;
0081 it=1;
0082 nw=-1;
0083 zp=zeros(1,p);
0084 r=(0:p);
0085 for jf=1:nf
0086     k(jf,1) = ceil(t1+t(it,3));
0087     k(jf,2) = ceil(t1+t(it,3)+t(it,2)-1);
0088     cs = (k(jf,1):k(jf,2)).';
0089     nc = length(cs);
0090     pp=min(p,nc);
0091     dd=s(cs);
0092     if nc~=nw
0093         % possibly we should have a window whose square integral equals unity
0094         ww=hamming(nc); nw=nc;
0095         y=zeros(1,nc+p);
0096         c=(1:nc)';
0097     end
0098     wd=dd(:).*ww;        % windowed data vector
0099     y(1:nc)=wd;          % data vector with p appended zeros
0100     z=zeros(nc,pp+1);    % data matrix
0101     %  was previously  z(:)=y(c(:,ones(1,pp+1))+r(ones(nc,1),1:pp+1));
0102     z(:)=y(repmat(c,1,pp+1)+repmat(r,nc,1));
0103     rr=wd'*z;
0104     rm=toeplitz(rr(1:pp));
0105     rk=rank(rm);
0106     if rk
0107         if rk<pp
0108             rm=rm(1:rk,1:rk);
0109         end
0110         ar(jf,2:rk+1)=-rr(2:rk+1)/rm;
0111     end
0112     e(jf)=rr*ar(jf,1:pp+1)';
0113     t1=t1+t(it,1);
0114     it=it+tr;
0115 end
0116