fconv.m

A MATLAB tool for analysis of Transient Otoacoustic Emission signals

function c = fconv(a, b)

%	c = fconv(a, b)
%
%	fconv	calculates the convolution of vectors a and b.
%
%   in:	a,b	vectors, polynomial coefficients, impulse responses etc.
%
%  out:	c	convolution of a and b

% (c) Heimo Ihalainen 25.7.1988 (14.2.1989 HI) (3.6.1989 HI) (23.8.1989 HI)
% (1.10.1989 HI new stb)
% (c) Timo Pyy 1.6.1989 (overlap-add)

% fconv uses FFT to achieve greater speed.

[na,n1] = size(a); ta = 0; if n1>na; a = a.'; na = n1; ta = 1; end
[nb,n1] = size(b); tb = 0; if n1>nb; b = b.'; nb = n1; tb = 1; end
t = ta; if nb>na; t = tb; end; nc = na+nb-1;
if t==0; c = zeros(nc,1); else c = zeros(1,nc); end; nm = min(na,nb);
%
if nm<100;					% direct without FFT
  if na>nb; c(:) = filter(b,1,[a;zeros(nc-na,1)]);
  else      c(:) = filter(a,1,[b;zeros(nc-nb,1)]);
  end
else; nf = 2^ceil(log(nc)/log(2)-10*eps);	% use FFT
  if (nm/nc>0.4)|(nc/nf>0.8)&(nf<4097)		% straight method
    cc = conj(fft(conj(fft([a;zeros(nf-na,1)]).*fft([b;zeros(nf-nb,1)]))));
  else						% overlap-add
    if na>=nb; nf = 2^ceil(log(nb)/log(2)+1-10*eps);
      B = fft([b;zeros(nf-nb,1)]); a(na+nf) = 0; cc = zeros(na+nf,1);
      inc = nf-nb; ina = 1:inc; ind = 1:nf; wc =zeros(nf,1); z = zeros(nb,1); 
      for i=1:ceil(na/inc);
	wc(:) = conj(fft(conj(fft([a(ina);z]).*B))); cc(ind) = cc(ind)+wc;
	ind(:) = ind+inc; ina(:) = ina+inc;
      end
    elseif nb>na; nf = 2^ceil(log(na)/log(2)+1-10*eps);
      A = fft([a;zeros(nf-na,1)]); b(nb+nf) = 0; cc = zeros(nb+nf,1);
      inc = nf-na; inb = 1:inc; ind = 1:nf; wc =zeros(nf,1); z = zeros(na,1);
      for i=1:ceil(nb/inc);
	wc(:) = conj(fft(conj(fft([b(inb);z]).*A))); cc(ind) = cc(ind)+wc;
	ind(:) = ind+inc; inb(:) = inb+inc;
      end
    end
  end
  if any(imag(a))|any(imag(b)); c(:) = cc(1:nc)/nf;
  else c(:) = real(cc(1:nc))/nf;
  end
end