# v_lpcifilt

## PURPOSE

V_LPCIFILT Apply inverse filter to speech signal U=(S,AR,T,DC,FADE)

## SYNOPSIS

function u=v_lpcifilt(s,ar,t,dc,fade)

## DESCRIPTION

    V_LPCIFILT Apply inverse filter to speech signal U=(S,AR,T,DC,FADE)

     Inputs:    S             speech signal
                    AR(nf,p+1)    array of ar coefficients; one row per frame
                    T             column vector giving the index of the first sample in each frame
               DC(nf)        optional dc component will be subtracted from the signal
               FADE          AR coefficients will be linearly interpolated for FADE samples
                             either side of frame boundaries

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_dypsa V_DYPSA Derive glottal closure instances from speech [gci,goi] = (s,fs)

## SOURCE CODE

    0001 function u=v_lpcifilt(s,ar,t,dc,fade)
    0002 %V_LPCIFILT Apply inverse filter to speech signal U=(S,AR,T,DC,FADE)
    0003 %
    0004 % Inputs:    S             speech signal
    0005 %                AR(nf,p+1)    array of ar coefficients; one row per frame
    0006 %                T             column vector giving the index of the first sample in each frame
    0007 %           DC(nf)        optional dc component will be subtracted from the signal
    0008 %           FADE          AR coefficients will be linearly interpolated for FADE samples
    0009 %                         either side of frame boundaries
    0010 %
    0011
    0012 % Example usage: generate an inverse filtered waveform
    0013 %
    0014 %      [sp,fs]=v_readwav('infile.wav');
    0015 %      lpcord=2+floor(fs/1000);
    0016 %      spp=filter([1 -exp(-2*pi*50/fs)],1,sp);                    % preemphasis zero is at 50 Hz
    0017 %      [lpar,lpe,lpk]=v_lpcauto(spp,lpcord,floor([0.01 0.02]*fs));  % 20ms frame with 10ms frame increment
    0018 %      overlap=lpk(1,2)-lpk(2,1)+1;                               % overlap between adjacent frames
    0019 %      u=v_lpcifilt(sp,lpar,lpk(:,1)+floor(overlap/2),0,overlap/4); % do inverse filtering
    0020 %      v_writewav(u,fs,'outfile.wav');
    0021
    0022 %      Copyright (C) Mike Brookes 1997
    0023 %      Version: $Id: v_lpcifilt.m 10865 2018-09-21 17:22:45Z dmb $
    0024 %
    0025 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0026 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0027 %
    0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0029 %   This program is free software; you can redistribute it and/or modify
    0030 %   it under the terms of the GNU Lesser General Public License as published by
    0031 %   the Free Software Foundation; either version 3 of the License, or
    0032 %   (at your option) any later version.
    0033 %
    0034 %   This program is distributed in the hope that it will be useful,
    0035 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0036 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0037 %   GNU Lesser General Public License for more details.
    0038 %
    0039 %   You can obtain a copy of the GNU Lesser General Public License from
    0040 %   https://www.gnu.org/licenses/ .
    0041 %
    0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0043
    0044 [nf,p1]=size(ar);
    0045 if nargin<4 | isempty(dc)
    0046    dc=zeros(nf,1);
    0047 elseif length(dc)==1
    0048    dc=dc(ones(nf,1));
    0049 end
    0050 if nf==1
    0051    u=filter(ar,1,s-dc);
    0052 else
    0053    p=p1-1;
    0054    ns=length(s);
    0055    if nargin<5 | isempty(fade)
    0056       fade=0;
    0057       if nargin<3 | isempty(t)
    0058          t=p1+(0:nf-1)*(ns-p)/nf;
    0059       end
    0060    end
    0061    u=zeros(ns,1);
    0062    if fade<1
    0063       x0=(max(1,ceil(t(nf)-p)):ns)';
    0064       u(x0)=filter(ar(nf,:),1,s(x0)-dc(nf));
    0065       for i=nf-1:-1:2
    0066          x0=(max(1,ceil(t(i)-p)):ceil(t(i+1)-1))';
    0067          u(x0)=filter(ar(i,:),1,s(x0)-dc(i));
    0068       end
    0069       x0=(1:ceil(t(2)-1))';
    0070       u(x0)=filter(ar(1,:),1,s(x0)-dc(1));
    0071
    0072    else
    0073       tb=min(t(nf)+fade,(t(nf)+ns)/2);
    0074       ta=max(t(nf)-fade,(t(nf-1)+t(nf))/2);
    0075       t0=max(1,ceil(ta)-p);
    0076       x0=(t0:ns)';
    0077       xb=x0-tb;
    0078       u(x0)=filter(ar(nf,:),1,s(x0)-dc(nf)).*max((1+(xb-abs(xb))/(2*(tb-ta))),0);
    0079       for i=nf-1:-1:2
    0080          td=tb;
    0081          tc=ta;
    0082          tb=min(t(i)+fade,ta);
    0083          ta=max(t(i)-fade,(t(i-1)+t(i))/2);
    0084          t1=floor(td);
    0085          t0=max(1,ceil(ta)-p);
    0086          x0=(t0:t1)';
    0087          xb=x0-tb;
    0088          xc=x0-tc;
    0089          u(x0)=u(x0)+filter(ar(i,:),1,s(x0)-dc(i)).*max((1+(xb-abs(xb))/(2*(tb-ta))-(xc+abs(xc))/(2*(td-tc))),0);
    0090       end
    0091       t1=floor(tb);
    0092       x0=(1:t1)';
    0093       xc=x0-ta;
    0094       u(x0)=u(x0)+filter(ar(1,:),1,s(x0)-dc(1)).*(1-(xc+abs(xc))/(2*(tb-ta)));
    0095
    0096    end
    0097 end
    0098