conv2 performance

[Lukas Reichlin]

> I would happily use a
> free software library with a GPL-compatible license to implement this
> function, but I don't know whether one is available.

I don't know whether DE01OD or DE01PD do exactly what you want, but since 
SLICOT (www.slicot.org) is GPL 2, these two convolution routines might be of 
interest.

hth
Lukas


DE01OD
Convolution or deconvolution of two real signals
[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To compute the convolution or deconvolution of two real signals
  A and B.


Specification

      SUBROUTINE DE01OD( CONV, N, A, B, INFO )
C     .. Scalar Arguments ..
      CHARACTER         CONV
      INTEGER           INFO, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(*), B(*)


Arguments

Mode Parameters

  CONV    CHARACTER*1
          Indicates whether convolution or deconvolution is to be
          performed as follows:
          = 'C':  Convolution;
          = 'D':  Deconvolution.


Input/Output Parameters

  N       (input) INTEGER
          The number of samples.  N must be a power of 2.  N >= 2.

  A       (input/output) DOUBLE PRECISION array, dimension (N)
          On entry, this array must contain the first signal.
          On exit, this array contains the convolution (if
          CONV = 'C') or deconvolution (if CONV = 'D') of the two
          signals.

  B       (input) DOUBLE PRECISION array, dimension (N)
          On entry, this array must contain the second signal.
          NOTE that this array is overwritten.


Error Indicator

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value.


Method

  This routine computes the convolution or deconvolution of two real
  signals A and B using an FFT algorithm (SLICOT Library routine
  DG01MD).


References

  [1] Rabiner, L.R. and Rader, C.M.
      Digital Signal Processing.
      IEEE Press, 1972.


Numerical Aspects

  The algorithm requires 0( N*log(N) ) operations.






DE01PD
Convolution or deconvolution of two real signals using Hartley transform
[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To compute the convolution or deconvolution of two real signals
  A and B using the Hartley transform.


Specification

      SUBROUTINE DE01PD( CONV, WGHT, N, A, B, W, INFO )
C     .. Scalar Arguments ..
      CHARACTER         CONV, WGHT
      INTEGER           INFO, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(*), B(*), W(*)


Arguments

Mode Parameters

  CONV    CHARACTER*1
          Indicates whether convolution or deconvolution is to be
          performed as follows:
          = 'C':  Convolution;
          = 'D':  Deconvolution.

  WGHT    CHARACTER*1
          Indicates whether the precomputed weights are available
          or not, as follows:
          = 'A':  available;
          = 'N':  not available.
          Note that if N > 1 and WGHT = 'N' on entry, then WGHT is
          set to 'A' on exit.


Input/Output Parameters

  N       (input) INTEGER
          The number of samples.  N must be a power of 2.  N >= 0.

  A       (input/output) DOUBLE PRECISION array, dimension (N)
          On entry, this array must contain the first signal.
          On exit, this array contains the convolution (if
          CONV = 'C') or deconvolution (if CONV = 'D') of the two
          signals.

  B       (input) DOUBLE PRECISION array, dimension (N)
          On entry, this array must contain the second signal.
          NOTE that this array is overwritten.

  W       (input/output) DOUBLE PRECISION array,
                         dimension (N - LOG2(N))
          On entry with WGHT = 'A', this array must contain the long
          weight vector computed by a previous call of this routine
          or of the SLICOT Library routine DG01OD.f, with the same
          value of N. If WGHT = 'N', the contents of this array on
          entry is ignored.
          On exit, this array contains the long weight vector.


Error Indicator

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value.


Method

  This routine computes the convolution or deconvolution of two
  real signals A and B using three scrambled Hartley transforms
  (SLICOT Library routine DG01OD).


References

  [1] Van Loan, Charles.
      Computational frameworks for the fast Fourier transform.
      SIAM, 1992.


Numerical Aspects

  The algorithm requires O(N log(N)) floating point operations.