Home > voicebox > lpccovar.m

lpccovar

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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: This function is called by:

SOURCE CODE ^

0001 function [ar,e,dc]=lpccovar(s,p,t,w)
0002 %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: lpccovar.m 8211 2016-07-20 20:59:16Z 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 General Public License as published by
0059 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0066 %
0067 %   You can obtain a copy of the GNU General Public License from
0068 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0069 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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)=sc.'*(sc - dm*aa.');
0139             e(jf,2)=sc.'*sc;
0140             e(jf,3:4)=e(jf,1:2)/nc;
0141         end
0142     end
0143 end
0144 if ~nargout
0145     lpcar2ff(repmat(sqrt(e(:,3).^(-1)),1,p+1).*ar,255);
0146     ylabel('Power (dB)');
0147 end
0148

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003