Code covered by the BSD License  

Highlights from
Quick Convolution

Quick Convolution

by

 

10 Jul 2013 (Updated )

Qconv is a faster alternative to the "conv" function of the SP Toolbox when the inputs are long

qconv(X, Y, P)
function Z = qconv(X, Y, P)
%  QCONV Quick linear convolution
%     Z = CONV(X, Y) convolves vectors X and Y.  The resulting vector is 
%     length LENGTH(X)+LENGTH(Y)-1, and is computed with the FFT.
%
%     Z = CONV(X,Y,P) convolves vectors X and Y through the FFT, where the
%     FFT size is chosen to be the minimum number with prime factors no
%     larger than P.  Choosing the FFT size in this way can dramatically
%     reduce the computation time.
%
%     If X and Y are vectors of polynomial coefficients, convolving
%     them is equivalent to multiplying the two polynomials.
%  
%     See also cfft conv
%
%
%
%    Christopher K. Turnes
%    Georgia Institute of Technology
%    Version:  1.0.0
%    Date: 9-Jul 2013
%
%

    % check for valid input
    if (nargin < 2)
        error('MATLAB:minrhs', 'Not enough input arguments.');
    end
    
    % get vector lengths
    N1 = length(X);
    N2 = length(Y);
    if (N1 ~= numel(X) || N2 ~= numel(Y))
        error('MATLAB:qconvfcn:AorBNotVector', 'X and Y must be vectors.');
    end
    if ((N1 == 0) || (N2 == 0))
        error('MATLAB:qconvfcn:vector', 'X and Y must be vectors.');
    end
    
    % check if min prime factor is supplied, start convolution
    if (nargin < 3)
        X = cfft(X, N1+N2-1);
    else
        
        if (isempty(P))
            warning('MATLAB:qconvfcn:noPrimeBound', ...
                'Third argument is empty; setting default prime bound to 5');
            P = 5;
        else
            if (length(P) > 1)
                warning('MATLAB:qconvfcn:primeBoundNotScalar', ...
                    'Prime bound must be integer scalar greater than 1.');
                P = P(1);
            end
            if (P < 2)
                error('MATLAB:qconvfcn:primeBoundTooSmall', ...
                    'Prime bound must be integer scalar greater than 1.');
            end
        end
        X = cfft(X, N1+N2-1, P);
    end
    
    % perform convolution
    Z = ifft(cfft(Y, length(X)) .* X);
    Z = Z(1:(N1+N2-1));

end

Contact us