/*****************************************************************************
*                                                                            *
*   XCORR.SPL   Copyright (C) 2000-2017 DSP Development Corporation          *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:     Randy Race                                                   *
*                                                                            *
*   Synopsis:   Cross-correlation using convolution                          *
*                                                                            *
*   Revisions:   2 May 2000  RRR  Creation - from FREQ.MAC                   *
*               29 Nov 2009  RRR  conv argument reversal                     *
*               26 Apr 2017  RRR  matrix support                             *
*                                                                            *
*****************************************************************************/


#if @HELP_XCORR

    XCORR

    Purpose: Cross correlation using the convolution method

    Syntax:  XCORR(s1, s2, norm)

                 s1 - input series

                 s2 - input series

               norm - optional integer, normalization method,

                      0: None,
                      1: Unity (-1 to 1)
                      2: Biased
                      3: Unbiased

                      defaults to 0: None

    Returns: A series

    Example:
             W1: gsin(1000, .001, 4)
             W2: gsin(1000, .001, 4)
             W3: xcorr(w1, w2)

             Performs the cross correlation of two sinewaves. The
             peaks of the result indicate the two waveforms are
             very similar at the time intervals where the peaks
             occur.

    Example:
             W1: gsin(1000, .001, 4)
             W2: gnorm(1000, .001)
             W3: xcorr(w1, w1, 1)
             W4: xcorr(w1, w2, 1)

             W3 displays the cross correlation of a sinewave normalized
             to -1 and 1. W4 shows the cross correlation between a sinewave
             and random noise. The normalized maximum of W3 is 1.0, indicating
             perfect correlation at time t == 0. Although the waveform of
             W4 displays some peaks, the normalized maximum is roughly 0.04
             indicating little correlation between W1 and W2. For
             a graphical representation, overplot W4 in W3.

    Remarks:
             The cross-correlation is used to determine how similar two
             series are to each other. XCORR performs correlation by
             computing the direct convolution of the input series.

             The output length L is:

                      L = length(s1) + length(s2) - 1

             For series of equal lengths and offsets, the zeroth lag
             component is the mid point of the series.

             The BIASED normalization divides the result by M, the
             maximum length of the input series.

             The UNBIASED normalization divides the result by

                              M - abs(M - i - 1) + 1

             where i is the index of the result.

             See FXCORR for the frequency domain implementation.


    See Also:
             Acorr
             Conv
             Facorr
             Fconv
             Fxcorr
#endif


/* time domain cross correlation */
xcorr(s1, s2, norm)
{
        local len, maxlen, t, dx, j, s3, u, snorm;

        if (argc < 3)
        {
                if (argc < 2) 
                {
                        if (argc < 1) 
                        {
                                error("xcorr - input series required");
                        }

                        s2 = refseries(s1);
                }

                if (isarray(s2))
                {
                        norm = 0;
                }
                else
                {
                        norm = s2;
                        s2   = refseries(s1);
                }
        }

        if (isstring(norm))
        {
                switch(tolower(norm))
                {
                        case "coeff":
                        case "unity":
                                norm = 1;
                                break;

                        case "biased":
                                norm = 2;
                                break;

                        case "unbiased":
                                norm = 3;
                                break;

                        default:
                                norm = 0;
                                break;
                }
        }

        if (numcols(s1) > 1 && numcols(s2) > 1)
        {
                /* all correlations column-wise */
                t  = {};
                s3 = conj(rev(s2));

                if (norm == 1)
                {
                        /* unity normalization for each column */
                        snorm = colsum(s2 * s2);

                        loop (j = 1.. numcols(s1))
                        {
                                u = conv(col(s1, j), s3);
                                u = u / sqrt(colsum(col(s1, j)^2) * snorm);

                                t = ravel(t, u);
                        }
                }
                else
                {
                        loop (j = 1..numcols(s1))
                        {
                                t = ravel(t, conv(col(s1, j), s3));
                        }

                        t = xcorr_norm(t, s1, s2, norm);
                }
        }
        else
        {
                /* perform correlation in the time domain */
                t = conv(s1, conj(rev(s2)));

                /* normalization */
                t = xcorr_norm(t, s1, s2, norm);
        }

        /* all done */
        return(t);
}


xcorr_norm(t, s1, s2, norm)
{
        local dx, maxlen;

        /* normalization */
        switch (norm)
        {
                case 0:        /* none */
                default:
                        break;

                case 1: /* -1 to +1 */
                        t = t / sqrt(colsum(s1 * s1) * colsum(s2 * s2));
                        break;

                case 2: /* biased */
                        maxlen = max({collength(s1), collength(s2)});
                        t = t / maxlen;
                        break;

                case 3: /* unbiased */
                        dx = deltax(t);
                        maxlen = max({length(s1), length(s2)});
                        t = t / ({1..maxlen, (maxlen - 1)..-1..1});

                        /* preserve original deltax */
                        setdeltax(t, dx);
                        break;
        }

        return(t);
}