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