function y=sconv(x,h)
%        FAST LINEAR CONVOLUTION OF A LONG AND SHORT SEQUENCE USING FFT
%                    (algorithm of overlap-add sectioning)
% function y=sconv(x,h)
% x ... long  input sequence - data to be filtered (vector)
% h ... short input sequence - finite impulse response of a filter (vector)
% y ... output sequence - filtered result as a linear convolution (y=x*h)
%       (column vector)

% (C) 1996, V. Vesely, Masaryk University of Brno (Czech Republic)
%
% Constants for the optimal buffer length of the overlapp-add sectioning method
% from:
% E.O.Brigham:"The Fast Fourier Transform and Its Applications",
% Prentice-Hall, 1988, p.217
OPTBUFCON = [1,11,20,30,50,100,200,300,600,1000,2000,4000];
LEN = length(OPTBUFCON);

error(nargchk(2,2,nargin));
m  = length(h);
r = find( m>=OPTBUFCON(1:LEN-1) & m<OPTBUFCON(2:LEN) );
if isempty(r),
   error(['The length of the short sequence ''h'' exceeded the upper limit ',...
          num2str(OPTBUFCON(LEN)-1)]);
end

N  = 2^(r+4);            % optimal buffer length
s  = N - m + 1;          % length of the input/output segment
n0 = length(x);
n  = ceil(n0/s)*s;       % the nearest multiple of the segment length
x  = [reshape(x,n0,1);zeros(n-n0,1)];     % pad x with zeros to the length n
H0 = fft([reshape(h,m,1);zeros(s-1,1)]);  % FFT of prolonged impulse response h
y  = zeros(n+m-1,1);
buf= zeros(N,1);         % clear buffer
realdata = isreal(x) & isreal(h);
for sk=1:s:n             % loop through segments
    X0 = fft([x(sk:sk+s-1);zeros(m-1,1)]); % FFT of prolonged input segment
    X0 = ifft(X0.*H0);                     % Convolution theorem
    if realdata, X0 = real(X0); end;
    buf = buf + X0;
    y(sk:sk+s-1) = buf(1:s);
    buf = [buf(s+1:N);zeros(s,1)];
end
y(n+1:n+m-1) = buf(1:m-1);
y = y(1:n0+m-1);         % discard extra zeros
