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.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 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)=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