/*****************************************************************************
*                                                                            *
*   FXCORR.SPL   Copyright (C) 2000 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Cross-correlation using FFT's                               *
*                                                                            *
*   Revisions:    2 May 2000  RRR  Creation - from FREQ.MAC                  *
*                                                                            *
*****************************************************************************/

#include 


#if @HELP_FXCORR

    FXCORR

    Purpose: Cross correlation using the FFT method

    Syntax:  FXCORR(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: fxcorr(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: fxcorr(w1, w1, 1)
             W4: fxcorr(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. FXCORR performs correlation by
             computing the FFT's 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 XCORR for the time domain implementation.

    See Also:
             Acorr
             Conv
             Facorr
             Fconv
             Xcorr
#endif


/* frequency domain cross correlation */
fxcorr(s1, s2, norm)
{
        local len, maxlen, f, xoff1, xoff2, dx;

        if (argc < 3)
        {
                if (argc < 2) 
                {
                        if (argc < 1) 
                        {
                                error("fxcorr - input series required");
                        }
                        
                        s2 = refseries(s1);
                }
                
                if (isarray(s2))
                {
                        norm = 0;
                }
                else
                {
                        norm = s2;
                        s2   = refseries(s1);
                }
        }

        /* get best fft length */
        maxlen = maxval(max(collength(s1)), max(collength(s2)));
        len    = bestpow2(maxlen * 2);

        /* perform correlation in the frequency domain */
        f = ifft(fft(s1, len) * fft(conj(rev(s2)), len));

        if (isreal(s1) && isreal(s2)) 
        {
                f = real(f);
        }

        /* extract proper output length */
        f = extract(f, 1, length(s1) + length(s2) - 1);

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

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

                case 2: /* biased */
                        maxlen = maxval(collength(s1), collength(s2));
                        f /= maxlen;
                        break;

                case 3: /* unbiased */
                        dx = deltax(f);
                        f /= ( {1..maxlen, (maxlen - 1)..-1..1});
                        
                        /* preserve original deltax */
                        setdeltax(f, dx);
                        break;

        }

        /* set proper xoffset */
        xoff1 = xoffset(s1);
        
        /* reversed offset */
        xoff2 = -xoffset(s2) - deltax(s2) * (length(s2) - 1);
        setxoffset(f, xoff1 + xoff2);
        
        return(f);
}