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

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 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

Generated by m2html © 2003