# v_lpcauto

## PURPOSE

V_LPCAUTO performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)

## SYNOPSIS

function [ar,e,k]=v_lpcauto(s,p,t,w,m)

## DESCRIPTION

    V_LPCAUTO  performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)
     Usage: (1) [ar,e]=v_lpcauto(x,p,[],'r','e'); same as lpc(x,p);

## CROSS-REFERENCE INFORMATION

This function calls:

- v_windows V_WINDOWS Generate a standard windowing function (TYPE,N,MODE,P,OV)

This function is called by:

- v_dypsa V_DYPSA Derive glottal closure instances from speech [gci,goi] = (s,fs)
- v_fxrapt V_FXRAPT RAPT pitch tracker [FX,VUV]=(S,FS,M,Q)

## SOURCE CODE

    0001 function [ar,e,k]=v_lpcauto(s,p,t,w,m)
    0002 %V_LPCAUTO  performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)
    0003 % Usage: (1) [ar,e]=v_lpcauto(x,p,[],'r','e'); same as lpc(x,p);
    0004
    0005 %  Inputs:
    0006 %     s(ns,nch)  is the input signal (ns samples, nch channels)
    0007 %      p          is the order (default: 12)
    0008 %      t(nf,3)    specifies the frames size details: each row specifies
    0009 %                 up to three values per frame: [hop anal skip] where:
    0010 %                   hop     is the length of the frame (default: length(s))
    0011 %                   anal    is the analysis length (default: hop)
    0012 %                   skip    is the number of samples to skip at the beginning (default: 0)
    0013 %                 If t contains only one row, it will be used repeatedly
    0014 %                 until there are no more samples left in s.
    0015 %     w          window or window type (see v_windows)
    0016 %                  Code Window     Sidelobe  3dB-BW  6dB-BW
    0017 %                  'c'   cos        -23dB     1.2     1.6  used in MP3
    0018 %                  'k'   kaiser(5)  -23dB     1.2     1.7  used in AAC
    0019 %                  'm'   hamming    -43dB     1.3     1.8  low sidelobes [default]
    0020 %                  'n'   hanning    -31dB     1.4     2.0  rapid falloff
    0021 %                  'r'   rectangle  -13dB     0.87    1.2  narrow main lobe
    0022 %                  'w'   rsqvorbis  -26dB     1.1     1.5  sqrt Vorbis
    0023 %                  's'   hamming    -24dB     1.1     1.5  sqrt Hamming
    0024 %                  'v'   vorbis     -21dB     1.3     1.8  used in Vorbis
    0025 %     m          set of mode options:
    0026 %                  'e'   scale window to make sum of squares = 1
    0027 %                  'j'   find a single set of LPC coefficients for all channels
    0028 %
    0029 % Outputs:
    0030 %     ar(nf,p+1,nch) are the AR coefficients with ar(:,1,:) = 1
    0031 %     e(nf,nch)      is the energy in the residual.
    0032 %                      sqrt(e) is often called the 'gain' of the filter.
    0033 %     k(nf,2)        gives the first and last sample of the analysis intervals
    0034
    0035 % Notes:
    0036 %
    0037 % (1) The first frame always starts at sample s(1) and the analysis window starts at s(t(1,3)+1).
    0038 % (2) The elements of t need not be integers; window parameters will be rounded to the nearest integers.
    0039 % (3) The analysis interval is always multiplied by a hamming window
    0040 % (4) As an example, if p=3 and t=[10 5 2], then the illustration below shows
    0041 %     successive frames labelled a, b, c, ... with capitals for the
    0042 %     analysis regions.
    0043 %
    0044 %      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 ...
    0045 %
    0046 %     The first frame starts at s(1) and the first analysis interval is t(1,3)+(1:t(1,2)).
    0047 %
    0048 % (5) Frames can overlap: e.g. t=[5 20] will use analysis frames of
    0049 %     length 20 overlapped by 15 samples.
    0050 % (6) For speech processing p should be at least 2*f*l/c where f is the sampling
    0051 %     frequency, l the vocal tract length and c the speed of sound. For a typical
    0052 %     male (l=17 cm) this gives f/1000.
    0053
    0054 %      Copyright (C) Mike Brookes 1997-2018
    0055 %      Version: $Id: v_lpcauto.m 10865 2018-09-21 17:22:45Z dmb $
    0056 %
    0057 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0058 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0059 %
    0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0061 %   This program is free software; you can redistribute it and/or modify
    0062 %   it under the terms of the GNU Lesser General Public License as published by
    0063 %   the Free Software Foundation; either version 3 of the License, or
    0064 %   (at your option) any later version.
    0065 %
    0066 %   This program is distributed in the hope that it will be useful,
    0067 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0068 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0069 %   GNU Lesser General Public License for more details.
    0070 %
    0071 %   You can obtain a copy of the GNU Lesser General Public License from
    0072 %   https://www.gnu.org/licenses/ .
    0073 %
    0074 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0075 persistent wnam wtypes
    0076 if isempty(wnam)
    0077     wnam={'c','k','m','n','r','w','s','v'};
    0078     wtypes=[10 11 3 2 1 8 3 7];
    0079 end
    0080 [ns,nch]=size(s);
    0081 if ns==1
    0082     s = s(:); % transpose if a row vector
    0083     [ns,nch]=size(s);
    0084 end
    0085 if nargin < 2 || isempty(p)
    0086     p=12;
    0087 end
    0088 if nargin < 3 || isempty(t)
    0089     t=ns; % default frame is entire signal
    0090 end
    0091 if nargin <4 || isempty(w)
    0092     w='m'; % default to Hamming window
    0093 end
    0094 if nargin<5 || isempty(m)
    0095     m='';
    0096 end
    0097 modee=any(m=='e'); % normalize window to unit energy
    0098 modej=any(m=='j'); % do LPC jointly for all channels
    0099 wch=ischar(w);
    0100 if wch
    0101     if modee
    0102         wopt='e';
    0103     else
    0104         wopt='';
    0105     end
    0106     wch=find(strcmp(w,wnam),1);
    0107 else
    0108     w=w(:);
    0109     if modee
    0110         w=w/sqrt(w'*w); % make sum of squares equal to 1
    0111     end
    0112 end
    0113 [nf,ng]=size(t);
    0114 if ng==1
    0115     t=[t t zeros(nf,1)];
    0116 elseif ng==2
    0117     t=[t zeros(nf,1)];
    0118 end
    0119 if nf==1
    0120     nf=floor(1+(ns-t(2)-t(3))/t(1));
    0121     t=repmat(t,nf,1); % repeat t for each frame
    0122 end
    0123 k=round(repmat(cumsum([0;t(1:nf-1,1)])+t(:,3),1,2)+[ones(nf,1) t(:,2)]); % integer start and end of each analysis frame
    0124 kd=k*[-1; 1]+1; % frame lengths
    0125 ku=unique(kd);
    0126 nk=length(ku); % number of unique frame lengths
    0127 if modej % do jointly over all channels
    0128     ar=zeros(p+1,nf); % space for LPC coefficients
    0129     e=zeros(nf,1); % space for residual energy
    0130     for ik=1:nk % loop for each unique frame length
    0131         kui=ku(ik); % analysis frame length
    0132         if wch
    0133             w = v_windows(wtypes(wch),kui,wopt,5)'; % only Kaiser needs a parameter; hence always = 5
    0134         end
    0135         pk=min(p,kui); % actual LPC order for these frames
    0136         km=kd==kui; % mask of frames with this length
    0137         nkm=sum(km); % number of frame with this length
    0138         sk=zeros(kui,nkm,nch); % space for data
    0139         sk(:)=s(repmat(repmat(k(km,1)',kui,1)+repmat((0:kui-1)',1,nkm),[1,1,nch])+reshape(repmat(ns*(0:nch-1),kui*nkm,1),[kui,nkm,nch])).*repmat(w(1:kui),[1,nkm,nch]);
    0140         rr=zeros(pk+1,nkm); % space for autocorrelation coefs
    0141         ark=rr;
    0142         rr(1,:)=sum(sum(sk.^2,1),3); % zero-lag autocorrelation
    0143         rr(2,:)=sum(sum(sk(1:kui-1,:,:).*sk(2:kui,:,:),1),3);
    0144         ark(1,:)=1; % first LPC coefficient is always 1
    0145         ark(2,:)=-rr(2,:)./rr(1,:);
    0146         ek=rr(1,:).*(ark(2,:).^2-1); % **** maybe force non-negative
    0147         for n=2:pk
    0148             rr(n+1,:)=sum(sum(sk(1:kui-n,:,:).*sk(n+1:kui,:,:),1),3); % calculate new autocorrelation term
    0149             ka=(rr(n+1,:)+sum(rr(n:-1:2,:).*ark(2:n,:),1))./ek; % ***** what if ek is zero, could limit to +-1
    0150             % mka=abs(ka)>=1; ka(mka)=sign(ka(mka));
    0151             ark(2:n,:)=ark(2:n,:)+repmat(ka,n-1,1).*ark(n:-1:2,:);
    0152             ark(n+1,:)=ka;
    0153             ek=ek.*(1-ka.^2);
    0154         end
    0155         ar(1:pk+1,km)=ark;
    0156         e(km)=-ek;
    0157     end
    0158     ar=permute(ar,[2 1 3]);
    0159 else                            % do each channel independently
    0160     ar=zeros(p+1,nf,nch);       % space for LPC coefficients
    0161     e=zeros(nf,nch);            % space for residual energy
    0162     for ik=1:nk                 % loop for each unique frame length
    0163         kui=ku(ik);             % analysis frame length
    0164         if wch
    0165             w = v_windows(wtypes(wch),kui,wopt,5)'; % only Kaiser needs a parameter; hence always = 5
    0166         end
    0167         pk=min(p,kui);          % actual LPC order for these frames
    0168         km=kd==kui;             % mask of frames with this length
    0169         nkm=sum(km);            % number of frame with this length
    0170         sk=zeros(kui,nkm*nch);  % space for data
    0171         sk(:)=s(repmat(repmat(k(km,1)',kui,1)+repmat((0:kui-1)',1,nkm),1,nch)+reshape(repmat(ns*(0:nch-1),kui*nkm,1),kui,nkm*nch)).*repmat(w(1:kui),1,nkm*nch);
    0172         rr=zeros(pk+1,nkm*nch);    % space for autocorrelation coefs
    0173         ark=rr;
    0174         rr(1,:)=sum(sk.^2,1);   % zero-lag autocorrelation
    0175         rr(2,:)=sum(sk(1:kui-1,:).*sk(2:kui,:),1);
    0176         ark(1,:)=1;             % first LPC coefficient is always 1
    0177         ark(2,:)=-rr(2,:)./rr(1,:);
    0178         ek=rr(1,:).*(ark(2,:).^2-1); % **** maybe force non-negative
    0179         for n=2:pk
    0180             rr(n+1,:)=sum(sk(1:kui-n,:).*sk(n+1:kui,:),1);      % calculate new autocorrelation term
    0181             ka=(rr(n+1,:)+sum(rr(n:-1:2,:).*ark(2:n,:),1))./ek;
    0182             % could limit to +-1 by doing: mka=abs(ka)>=1; ka(mka)=sign(ka(mka));
    0183             ark(2:n,:)=ark(2:n,:)+repmat(ka,n-1,1).*ark(n:-1:2,:);
    0184             ark(n+1,:)=ka;
    0185             ek=ek.*(1-ka.^2);
    0186         end
    0187         ar(1:pk+1,km,:)=reshape(ark,[pk+1,nkm,nch]);    % insert into output array
    0188         e(km,:)=reshape(-ek,nkm,nch);                   % insert into output array
    0189     end
    0190     ar=permute(ar,[2 1 3]);
    0191 end