/*****************************************************************************
*                                                                            *
*   FCONV2D.SPL   Copyright (C) 2011 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    2D Convolution using FFT's                                  *
*                                                                            *
*   Revisions:   24 Feb 2011  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_FCONV2D

    FCONV2D

    Purpose: Computes the 2D convolution of two arrays using the FFT method.

    Syntax:  FCONV2D(a1, a2)

                 a1 - input array
                 a2 - input array

    Returns: An array, the convolved result.

    Example:
             W1: randn(30)
             W2: randn(30)
             W3: conv2d(w1, w2)
             W4: fconv2d(w1, w2)

             Performs the 2D convolution of two random arrays. W3 performs
             direct convolution in the time domain and W4 performs the
             same convolution using the FFT method.

    Remarks:
             FCONV2D performs convolution by computing the FFT's of the
             input arrays. This method is usually faster than CONV2D when
             both inputs are large arrays.

             See CONV2D for the time domain implementation.

    See Also:
             Conv
             Conv2d
             Fconv
#endif


/* frequency domain convolution */
fconv2d(s1, s2)
{
        local s, r1, c1, r2, c2;

        (r1, c1) = size(s1);
        (r2, c2) = size(s2);

        /* find best FFT length */
        rlen = bestpow2(2 * maxval(r1, r2));
        clen = bestpow2(2 * maxval(c1, c2));

        /* transform, multiply and inverse transform */
        s = ifft2(fft2(s1, rlen, clen) * fft2(s2, rlen, clen));

        /* resize */
        s = s[1..(r1 + r2 - 1), 1..(c1 + c2 - 1)];

        /* check if result should be purely real */
        if (isreal(s1) && isreal(s2))
        {
                s = real(s);
        }

        /* xoffset */
        setxoffset(s, xoffset(s1) + xoffset(s2));

        /* yoffset */
        setyoffset(s, yoffset(s1) + yoffset(s2));

        return(s);
}