# v_lpccovar

## PURPOSE

V_LPCCOVAR performs covariance LPC analysis [AR,E,DC]=(S,P,T)

## SYNOPSIS

function [ar,e,dc]=v_lpccovar(s,p,t,w)

## DESCRIPTION

    V_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,4)     each row is [Er Es Pr Ps] and gives the energy ("E") and power ("P")
                           in the input signal window ("s") and in the LPC residual "r".
                           The 'gain' of the LPC filter is g=sqrt(Pr); x=filter(g,ar,randn(:,1)) will
                           generate noise with approximately the same power spectrum as the input s.
               DC          is the DC component of the signal S. If this output is included,
                           the LPC equations are modified to include a DC offset.

## CROSS-REFERENCE INFORMATION

This function calls:

- v_lpcar2ff V_LPCAR2FF LPC: Convert AR coefs to complex spectrum FF=(AR,NP)

This function is called by:

## SOURCE CODE

    0001 function [ar,e,dc]=v_lpccovar(s,p,t,w)
    0002 %V_LPCCOVAR performs covariance LPC analysis [AR,E,DC]=(S,P,T)
    0003 %
    0004 %  Inputs:  S(NS)    is the input signal
    0005 %           P        is the order (default: 12)
    0006 %           T(NF,:)  specifies the frames size details: each row specifies one frame
    0007 %                    T can be a cell array if rows have unequal numbers of values
    0008 %                       T(:,1) gives the start of the analysis interval: must be >P
    0009 %                       T(:,2) gives the end of the anaylsis interval [default: t(:+1,1)-1]
    0010 %                       subsequent pairs can be used to specify multiple disjoint segments
    0011 %                    If T is omitted, T(1,1)=P+1, T(1,2)=NS;
    0012 %                    The elements of t need not be integers.
    0013 %           W(NS)    The error at each sample is weighted by W^2 (default: 1)
    0014 %
    0015 % Outputs:  AR(NF,P+1)  are the AR coefficients with AR(:,1) = 1
    0016 %           E(NF,4)     each row is [Er Es Pr Ps] and gives the energy ("E") and power ("P")
    0017 %                       in the input signal window ("s") and in the LPC residual "r".
    0018 %                       The 'gain' of the LPC filter is g=sqrt(Pr); x=filter(g,ar,randn(:,1)) will
    0019 %                       generate noise with approximately the same power spectrum as the input s.
    0020 %           DC          is the DC component of the signal S. If this output is included,
    0021 %                       the LPC equations are modified to include a DC offset.
    0022
    0023 % Notes:
    0024 %
    0025 % (1a) If no DC output is specified AR(j,:)*S(n-(0:P)) ~ 0 or, equivalently,
    0026 %      S(n) ~ -AR(j,2:P)*S(n-(1:P)) where T(j,1) <= n <= T(j,2).
    0027 % (1b) If a DC output is specified AR(j,:)*(S(n-(0:P))-DC) ~ 0 or, equivalently,
    0028 %      S(n) ~ DC - AR(j,2:P)*(S(n-(1:P))-DC) = DC*sum(AR,j,:)) - AR(j,2:P)*S(n-(1:P))
    0029 %      where T(j,1) <= n <= T(j,2).
    0030 %
    0031 % (2) For speech processing P should be at least 2*F*L/C where F is the sampling
    0032 %     frequency, L the vocal tract length and C the speed of sound. For a typical
    0033 %     male (l=17 cm) this gives f/1000.
    0034 %
    0035 % (3) Each analysis frame should contain at least 2P samples. If note (1) is followed
    0036 %     this implies at least 2 ms of speech signal per frame.
    0037 %
    0038 % (4) It can be advantageous to restrict the analysis regions to time intervals
    0039 %     when the glottis is closed (closed-phase analysis). This can be achieved by
    0040 %     setting the T input parameter appropriately. If the closed-phase is shorter than
    0041 %     2 ms then two or more successive closed-phases should be used by defining 4 or more
    0042 %     elements in the corresponding row of T.
    0043 %
    0044 % (5) A previous version of this routine allowed T() to have a single row which would
    0045 %     be replicated for the entire file length. This has been removed because it gave rise
    0046 %     to an ambiguity.
    0047
    0048 %  Bugs: should really detect a singular matrix and reduce the order accordingly
    0049
    0050 %       Copyright (C) Mike Brookes 1995
    0051 %      Version: $Id: v_lpccovar.m 10865 2018-09-21 17:22:45Z dmb $
    0052 %
    0053 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0054 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0055 %
    0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0057 %   This program is free software; you can redistribute it and/or modify
    0058 %   it under the terms of the GNU Lesser General Public License as published by
    0059 %   the Free Software Foundation; either version 3 of the License, or
    0060 %   (at your option) any later version.
    0061 %
    0062 %   This program is distributed in the hope that it will be useful,
    0063 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0064 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0065 %   GNU Lesser General Public License for more details.
    0066 %
    0067 %   You can obtain a copy of the GNU Lesser General Public License from
    0068 %   https://www.gnu.org/licenses/ .
    0069 %
    0070 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0071 s = s(:); % make it a column vector
    0072 if nargin < 2 p=12; end;
    0073 if nargin < 3 t=[p+1 length(s)]; end;
    0074 wq = nargin>3;
    0075 [nf,ng]=size(t);
    0076 if iscell(t)
    0077     t{nf+1}=length(s)+1;
    0078 else
    0079     if rem(ng,2)
    0080         t(:,end+1)=[t(2:nf,1)-1; length(s)];
    0081     end
    0082 end
    0083 ar=zeros(nf,p+1);
    0084 ar(:,1)=1;
    0085 e=zeros(nf,4);
    0086 dc=zeros(nf,1);
    0087 d0=nargout >2;
    0088 rs=(1:p);
    0089 for jf=1:nf
    0090     if iscell(t)
    0091         tj=t{jf};
    0092         if rem(length(tj),2)
    0093             tj(end+1)=t{jf+1}(1)-1;
    0094         end
    0095     else
    0096         tj=t(jf,:);
    0097     end
    0098
    0099     ta = ceil(tj(1));
    0100     tb = floor(tj(2));
    0101     cs = (ta:tb).';
    0102     for js=3:2:length(tj)
    0103         ta = ceil(tj(js));
    0104         tb = floor(tj(js+1));
    0105         cs = [cs; (ta:tb).'];
    0106     end
    0107     %disp(cs([logical(1); (cs(2:end-1)~=cs(1:end-2)+1)|(cs(2:end-1)~=cs(3:end)-1); logical(1)])');
    0108     nc = length(cs);
    0109     pp=min(p,nc-d0);
    0110     dm=zeros(nc,pp);    % predefine shape
    0111     dm(:) = s(cs(:,ones(1,pp))-rs(ones(nc,1),1:pp));
    0112     if nargout>2
    0113         if wq
    0114             dm = [ones(nc,1) dm].*w(cs(:,ones(1,1+pp)));
    0115             sc=(s(cs).*w(cs));
    0116             aa = (dm\sc).';
    0117         else
    0118             dm = [ones(nc,1) dm];
    0119             sc=s(cs);
    0120             aa = (dm\sc).';
    0121         end
    0122         ar(jf,2:pp+1) = -aa(2:pp+1);
    0123         e(jf,1)=sc.'*(sc - dm*aa.');
    0124         e(jf,2)=sc.'*sc;
    0125         e(jf,3:4)=e(jf,1:2)/nc;
    0126         dc(jf)=aa(1)/sum(ar(jf,:));
    0127     else
    0128         if wq
    0129             dm = dm.*w(cs(:,ones(1,pp)));
    0130             sc=(s(cs).*w(cs));
    0131             aa = (dm\sc).';
    0132         else
    0133             sc=s(cs);
    0134             aa = (dm\sc).';
    0135         end;
    0136         ar(jf,2:pp+1) = -aa;
    0137         if nargout~=1
    0138             e(jf,1)=real(sc'*(sc - dm*aa.'));
    0139             e(jf,2)=real(sc'*sc);
    0140             e(jf,3:4)=e(jf,1:2)/nc;
    0141         end
    0142     end
    0143 end
    0144 if ~nargout
    0145     v_lpcar2ff(repmat(sqrt(e(:,3).^(-1)),1,p+1).*ar,255);
    0146     ylabel('Power (dB)');
    0147 end
    0148