# v_convfft

## PURPOSE

V_CONFFT 1-D convolution or correlation using FFT

## SYNOPSIS

function z=v_convfft(x,h,d,m,h0,x1,x2)

## DESCRIPTION

    V_CONFFT 1-D convolution or correlation using FFT

      Usage: (1) z=v_convfft(x,h,d1,'',1,1,size(x,d)+length(h)-1);            % equivalent to conv(x,h)

             (2) hh=v_convfft(length(x),h,2,'z',1,1,length(x)+length(h)-1);   % precalculate fft(h) for multiple calls
                 z=v_convfft(x,hh);                                           % also equivalent to conv(x,h)

             (3) z=v_convfft(x,h);                                            % equivalent to filter(h,1,x)
                 z=v_convfft(x,h,d1,'',1,1,size(x,d));                        % equivalent to filter(h,1,x)

             (4) z=v_convfft(x,h,d,'X',1,2-max(size(x,d),length(h)),max(size(x,d),length(h)));  % equivalent to xcorr(x,h)

             (5) z=v_convfft(x,h,d,'x',floor((1+length(h))/2),1,size(x,d));   % equivalent to imfilter(x,h)

      Inputs: x     input vector or array (or size(x,d) for the 'z' option)
              h     impulse response (or z output from previous call with the 'z' option)
              d     dimension of x to do convolution along [first non-singleton]
              m     mode options (see below)
              h0    origin sample number in h (can be outside the range 1:length(h)) [default: 1]
              x1,x2 range of rows/columns in x to align with origin of h (can be outside the range 1:size(x,d)) [default: 1,size(x,d)]

     Outputs: z     output from convolution/correlation. The same size and shape as x except that size(z,d)=x2-x1+1
                    If m='z' is specified, then z is a structure that can be used as h in a subsequent call

     Mode is any sensible combination of the following
              'x'   perform real correlation rather than convolution (reflects h around sample h0)
              'X'   perform complex correlation rather than convolution (reflects and conjugates h around sample h0)
              'z'   Precalculate fft(h) for efficiency. Input d must be given explicitly and input x equals size(x,d).

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function z=v_convfft(x,h,d,m,h0,x1,x2)
    0002 %V_CONFFT 1-D convolution or correlation using FFT
    0003 %
    0004 %  Usage: (1) z=v_convfft(x,h,d1,'',1,1,size(x,d)+length(h)-1);            % equivalent to conv(x,h)
    0005 %
    0006 %         (2) hh=v_convfft(length(x),h,2,'z',1,1,length(x)+length(h)-1);   % precalculate fft(h) for multiple calls
    0007 %             z=v_convfft(x,hh);                                           % also equivalent to conv(x,h)
    0008 %
    0009 %         (3) z=v_convfft(x,h);                                            % equivalent to filter(h,1,x)
    0010 %             z=v_convfft(x,h,d1,'',1,1,size(x,d));                        % equivalent to filter(h,1,x)
    0011 %
    0012 %         (4) z=v_convfft(x,h,d,'X',1,2-max(size(x,d),length(h)),max(size(x,d),length(h)));  % equivalent to xcorr(x,h)
    0013 %
    0014 %         (5) z=v_convfft(x,h,d,'x',floor((1+length(h))/2),1,size(x,d));   % equivalent to imfilter(x,h)
    0015 %
    0016 %  Inputs: x     input vector or array (or size(x,d) for the 'z' option)
    0017 %          h     impulse response (or z output from previous call with the 'z' option)
    0018 %          d     dimension of x to do convolution along [first non-singleton]
    0019 %          m     mode options (see below)
    0020 %          h0    origin sample number in h (can be outside the range 1:length(h)) [default: 1]
    0021 %          x1,x2 range of rows/columns in x to align with origin of h (can be outside the range 1:size(x,d)) [default: 1,size(x,d)]
    0022 %
    0023 % Outputs: z     output from convolution/correlation. The same size and shape as x except that size(z,d)=x2-x1+1
    0024 %                If m='z' is specified, then z is a structure that can be used as h in a subsequent call
    0025 %
    0026 % Mode is any sensible combination of the following
    0027 %          'x'   perform real correlation rather than convolution (reflects h around sample h0)
    0028 %          'X'   perform complex correlation rather than convolution (reflects and conjugates h around sample h0)
    0029 %          'z'   Precalculate fft(h) for efficiency. Input d must be given explicitly and input x equals size(x,d).
    0030
    0031 %      Copyright (C) Mike Brookes 2016-2017
    0032 %      Version: $Id: v_convfft.m 10865 2018-09-21 17:22:45Z dmb $
    0033 %
    0034 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0035 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0036 %
    0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0038 %   This program is free software; you can redistribute it and/or modify
    0039 %   it under the terms of the GNU Lesser General Public License as published by
    0040 %   the Free Software Foundation; either version 3 of the License, or
    0041 %   (at your option) any later version.
    0042 %
    0043 %   This program is distributed in the hope that it will be useful,
    0044 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0045 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0046 %   GNU Lesser General Public License for more details.
    0047 %
    0048 %   You can obtain a copy of the GNU Lesser General Public License from
    0049 %   https://www.gnu.org/licenses/ .
    0050 %
    0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0052 if ~isstruct(h) % normal input calling sequence
    0053     if nargin<4 || isempty(m)
    0054         m='';                       % set default mode string
    0055     end
    0056     if nargin<5 || isempty(h0)
    0057         h0=1;                       % set default h origin
    0058     end
    0059     s=size(x);                      % get the structure of x
    0060     ps=numel(x);                      % number of elements in x
    0061     ns=length(s);                   % number of dimensions in x
    0062     if any(m=='z')      % output pre-computed structure instead of convolution
    0063         if nargin<3 || isempty(d)   % if d unspecified
    0064             error('d must be specified explicitly');
    0065         end
    0066         if numel(x)~=1              % x contains nx for a future call
    0067             error('x must equal size(*,d)');
    0068         end
    0069         nx=x;                       % save x as the value of nx
    0070     else
    0071         if nargin<3 || isempty(d)       % if d unspecified
    0072             if ps<2
    0073                 d=1;                    % if d is a singleton or is empty
    0074             else
    0075                 d=find(s>1,1);          % d = first nonsingleton dimension
    0076             end
    0077         end
    0078         nx=s(d);                        % length in correlation dimension
    0079     end
    0080     k=ps/nx;                        % total size of all other dimensions
    0081     if nargin<6 || isempty(x1)
    0082         x1=1;                       % default initial lag
    0083     end
    0084     if nargin<7 || isempty(x2)
    0085         x2=nx;                      % default final lag
    0086     end
    0087     nz=x2-x1+1;                     % number of output lags required
    0088     h=h(:);                         % force h to be a column
    0089     nh=length(h);
    0090     if any(m=='X')                  % do complex correlation rather than convolution
    0091         h=conj(h(nh:-1:1));         % reflect and conjugate h
    0092         h0=nh+1-h0;                 % reflect the position of h0
    0093     elseif any(m=='x')              % do real correlation
    0094         h=h(nh:-1:1);               % reflect h
    0095         h0=nh+1-h0;                 % reflect the position of h0
    0096     end
    0097     hmin=h0+x1-nx;                  % smallest h index ever used
    0098     hmax=h0+x2-1;                   % largest h index ever used
    0099     xmin=x1+h0-nh;                  % smallest x index ever used
    0100     xmax=x2+h0-1;                   % largest x index ever used
    0101     if hmin>1 || hmax<nh            % we can delete some unused h values
    0102         hmin=max(hmin,1);
    0103         h=h(hmin:min(hmax,nh));     % trim h if possible
    0104         nh=length(h);
    0105         h0=h0-hmin+1;               % update h0 to new position
    0106     end
    0107     if xmin>1 || xmax<nx            % we can delete some unused x values
    0108         vmin=max(xmin,1);           % we will trim v to v(vmin:vmax)
    0109         vmax=min(xmax,nx);
    0110         x1=x1-vmin+h0;              % update x1,x2 to new positions assuming h0=0
    0111         x2=x2-vmin+h0;
    0112     else
    0113         vmin=1;
    0114         vmax=nx;
    0115         x1=x1+h0-1;                 % update x1,x2 to new positions assuming h0=0
    0116         x2=x2+h0-1;
    0117     end
    0118     nv=vmax-vmin+1;                 % number of elements of v to retain
    0119     nxz=min(max(max(nh-x1,0),max(x2-nv,0)),nh-1); % number of zeros to add to v is the larger of the number needed at each end
    0120     [fnx,enx]=log2(nv+nxz);         % round up length of zero-padded v to next power of 2
    0121     nf=pow2(1,enx-(fnx==0.5));      % actual length of dft is a power of 2
    0122     fmin=max(x1,1);                 % range of indices to extract from circular convolution
    0123     fmax=min(x2,min(nf,nx+nh-1));
    0124     zmin=max(1,2-x1);               % range of indices for non-zero output values
    0125     zmax=zmin+fmax-fmin;
    0126     fh=fft(h,nf,1);
    0127 else                % h is the z output of a previous call with the 'z' option
    0128     d=h.d;          % x dimension to convolve along
    0129     nx=h.nx;        % original size(x,d)
    0130     ns=h.ns;        % number of dimensions of x
    0131     nv=h.nv;        % number x values needed (might be < nx)
    0132     vmin=h.vmin;    % first x value needed (might be > 1)
    0133     vmax=h.vmax;    % last x value needed (might be < nx)
    0134     nf=h.nf;        % dft length
    0135     fh=h.fh;        % dft of impulse response h
    0136     fmin=h.fmin;    % first fh value needed
    0137     fmax=h.fmax;    % last fh value needed
    0138     nz=h.nz;        % output size(z,d)
    0139     zmin=h.zmin;    % first non-zero z value
    0140     zmax=h.zmax;    % last non-zero z value
    0141     s=size(x);      % get the structure of x
    0142     k=numel(x)/nx;  % total size of the rest of x
    0143     if size(x,d)~=nx || length(s)~=ns
    0144         error('input x has incompatible dimensions');
    0145     end
    0146 end
    0147 if  ~isstruct(h) && any(m=='z')    % save required information as a structure for next call
    0148     z.d=d;          % x dimension to convolve along
    0149     z.nx=nx;        % original size(x,d)
    0150     z.ns=ns;        % number of dimensions of x
    0151     z.nv=nv;        % number x values needed (might be < nx)
    0152     z.vmin=vmin;    % first x value needed (might be > 1)
    0153     z.vmax=vmax;    % last x value needed (might be < nx)
    0154     z.nf=nf;        % dft length
    0155     z.fh=fh;        % dft of impulse response h
    0156     z.fmin=fmin;    % first fh value needed
    0157     z.fmax=fmax;    % last fh value needed
    0158     z.nz=nz;        % output size(z,d)
    0159     z.zmin=zmin;    % first non-zero z value
    0160     z.zmax=zmax;    % last non-zero z value
    0161 elseif ~isempty(x)  % there is data to convolve
    0162     if d==1         % reshape x if necessary
    0163         v=reshape(x,nx,k);
    0164     else
    0165         v=reshape(permute(x,[d:ns 1:d-1]),nx,k);
    0166     end
    0167     if nv<nx        % we can delete some unused v values
    0168         v=v(vmin:vmax,:);
    0169     end
    0170     zz=ifft(fft(v,nf,1).*repmat(fh,1,k));    % do the convolution
    0171     z=zeros(nz,k);                  % reserve space for the output
    0172     z(zmin:zmax,:)=zz(fmin:fmax,:); % extract values from the convolution
    0173     s(d)=nz;                        % update the size of the output
    0174     if d==1
    0175         z=reshape(z,s);
    0176     else
    0177         z=permute(reshape(z,s([d:ns 1:d-1])),[ns+2-d:ns 1:ns+1-d]);
    0178     end
    0179 else
    0180     z=[];           % no input data
    0181 end
    0182
    0183
    0184