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

cfft(varargin)
function Y = cfft(varargin)
%  CFFT Customized Discrete Fourier transform.
%     CFFT(X) is the discrete Fourier transform (DFT) of vector X.  For
%     matrices, the FFT operation is applied to each column. For N-D
%     arrays, the FFT operation operates on the first non-singleton
%     dimension.
%  
%     CFFT(X,N) is the N-point FFT, padded with zeros if X has less
%     than N points and truncated if it has more.
%
%     CFFT(X,N,P) will take the M-point FFT, where M is the next-largest
%     number to N that is factorable by prime numbers no larger than P.
%  
%     CFFT(X,[],[],DIM) or CFFT(X,N,[],DIM) or CFFT(X,N,P,DIM) applies the 
%     FFT operation across the dimension DIM.
%
%     
%     For length N input vector x, the DFT is a length N vector X,
%     with elements
%                      N
%        X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                     n=1
%     The inverse DFT (computed by IFFT) is given by
%                      N
%        x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
%                     k=1
%  
%     See also fft fft2, fftn, fftshift, fftw, ifft, ifft2, ifftn.
%
%
%
%    Christopher K. Turnes
%    Georgia Institute of Technology
%    Version:  1.0.0
%    Date: 9-Jul 2013
%
%

    
    % check for valid input
    if (length(varargin) < 1)
        error('MATLAB:minrhs', 'Not enough input arguments.');
    end
    
    % default to standard FFT if only two input arguments
    if (length(varargin) == 1)
        Y = fft(varargin{1});
        return;
    end
    if (length(varargin) == 2)
        Y = fft(varargin{1}, varargin{2});
        return;
    end
    
    % data is first argument, size is second
    X = varargin{1};
    N = varargin{2};
    if (isempty(N))
        N = size(X, 1);
    elseif (length(N) > 1)
        warning('MATLAB:cfftfcn:lengthNotNonNegInt', ...
            'FFT length must be a non-negative integer scalar.');
        N = N(1);
    else
        if (N < 0)
            error('MATLAB:cfftfcn:lengthNotNonNegInt', ...
                'FFT length must be a non-negative integer scalar.');
        end
    end

    % third argument 
    P = varargin{3};
    M = N;
    if (~isempty(P))
        
        if (length(P) > 1)
            warning('MATLAB:cfftfcn:primeBoundNotScalar', ...
                'Prime bound must be integer scalar greater than 1.');
            P = P(1);
        end
        if (P < 2)
            error('MATLAB:cfftfcn:primeBoundTooSmall', ...
                'Prime bound must be integer scalar greater than 1.');
        end
        
        % check factors
        f = factor(N);
        if (max(f) > P);
        
            % build quick lookup table
            T = 1;

            % go through all primes up to pmax
            for k = 2:P

                % check if prime
                if (~isprime(k))
                    continue;
                end

                % build viable range:
                km = ceil(log(N) / log(k));

                % product of powers with powers of existing primes
                T = T * k.^(0:km);

                % clip off anything above the minimal value above N
                % (saves pointless calculation)
                T = T(:);
                tm = min(T(T >= N));
                T = T(T <= tm);

            end

            % select minimum value above N
            M = min(T(T >= N));
            
        end
    end
    
    % call FFT at this point
    if (length(varargin) == 3)
        Y = fft(X,M);
        return;
    elseif (length(varargin) == 4)
        Y = fft(X,M,varargin{4});
    else
        error('MATLAB:maxrhs', 'Too many input arguments.');
    end
    
end

Contact us