/*****************************************************************************
*                                                                            *
*   DECILP.SPL   Copyright (C) 1998-2009 Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Low pass filters and decimates a series                     *
*                                                                            *
*   Revisions:    9 Jun 1998  RRR  Creation                                  *
*                22 Dec 2009  RRR  Use series end padding                    *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_DECILP

    DECILP

    Purpose: Low pass filters and decimates a series

    Syntax:  DECILP(s, n, edge, attn)

              s    - the input series
              n    - an optional integer, decimation factor (default 2)
              edge - an optional real, the transition band of the low pass
                     filter (default 0.8125)
              attn - an optional real, the attenuation in dB of the
                     low pass filter (default 70.0 db)


    Returns: A series or array

    Example:
             W1: integ(gnorm(1000, 1/1000))
             W2: decilp(W1, 5)
             W3: spectrum(W2, 8*1024)
             W4: spectrum(W1, 8*1024);overp(W3, lred);setylog(1, 0, 1)

             The spectrum in W3 is a virtually the first 1/5 of the
             spectrum of the original data. The length of the output
             is 1/5 the length of the input.
  
  Example:
             W2: Decilp(W1, 5, 0.95)

             Same as above, but the transition band is narrowed.

    Remarks:
             DECILP creates an FIR low pass filter using the Kaiser
             Window method. Use

                      (out, lpf) = decilp(s, n)

             to return both the output data and the coefficients of the
             FIR low pass filter.

             The edge frequencies of the filter are determined as
             follows:

                     Fcut = edge * Fstop
                     Fstop = Fs / (2*n)

             Where Fs is the sample rate of the data (i.e. RATE(s)).

             The decimation factor automatically adjusts the sampling rate (1/deltax) 
             of the resulting series.

    See Also:
             DADiSP/Filters
             Decimate
             Padfilt

#endif


/* decimate with low pass filtering to preserve frequency spectrum */
ITERATE decilp(s, n, edge, attn)
{
        local r, fc, f, out;

        /* default args */
        if (argc < 4)
        {
                attn = 60.0;
                
                if (argc < 3)
                {
                        edge = 0.8125;
                        
                        if (argc < 2)
                        {
                                n = 2;
                        }
                }
        }

        r = rate(s);
        n = int(n);

        if (n < 2)
        {
                if (n == 1)
                {
                        return(s);
                }
                else
                {
                        error("decilp - positive decimation factor required");
                }
        }
        
        if (edge <= 0.0 || edge >= 1.0)
        {
                error(sprintf("decilp - edge is %g, must lie between 0.0 and 1.0", edge));
        }
        
        if (attn <= 0.0)
        {
                error("decilp - positive attenuation required");
        }

        /* cutoff freq is 1/n * nyquist frequency */
        fc = r / (2 * n);

        /* build a Kaiser FIR low pass filter */
//        f = firlp(r, edge * fc, fc, attn);
        f = firlp(r, edge * fc, fc * (2 - edge), attn);

        /* FIR filter and decimate it */
        out = decimate(firflt(s, f, 1), n);

        if (outargc > 1)
        {
                /* return result and filter */
                return(out, f);
        }
        else
        {
                return(out);
        }
}