Home > voicebox > convfft.m

convfft

PURPOSE ^

CONFFT 1-D convolution or correlation using FFT

SYNOPSIS ^

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

DESCRIPTION ^

CONFFT 1-D convolution or correlation using FFT

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

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

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

         (4) z=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=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=convfft(x,h,d,m,h0,x1,x2)
0002 %CONFFT 1-D convolution or correlation using FFT
0003 %
0004 %  Usage: (1) z=convfft(x,h,d1,'',1,1,size(x,d)+length(h)-1);            % equivalent to conv(x,h)
0005 %
0006 %         (2) hh=convfft(length(x),h,2,'z',1,1,length(x)+length(h)-1);   % precalculate fft(h) for multiple calls
0007 %             z=convfft(x,hh);                                           % also equivalent to conv(x,h)
0008 %
0009 %         (3) z=convfft(x,h);                                            % equivalent to filter(h,1,x)
0010 %             z=convfft(x,h,d1,'',1,1,size(x,d));                        % equivalent to filter(h,1,x)
0011 %
0012 %         (4) z=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=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: convfft.m 10118 2017-09-17 19:45:51Z 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 General Public License as published by
0040 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0047 %
0048 %   You can obtain a copy of the GNU General Public License from
0049 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0050 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003