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);
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 General Public License as published by 0063 % the Free Software Foundation; either version 2 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 General Public License for more details. 0070 % 0071 % You can obtain a copy of the GNU General Public License from 0072 % http://www.gnu.org/copyleft/gpl.html or by writing to 0073 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 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