/*****************************************************************************
*                                                                            *
*   ZLPFILT.SPL Copyright (C) 2008-2012 DSP Development Corporation          *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Low pass filtering by FFT zeroing                           *
*                                                                            *
*   Revisions:   18 Jul 2008  RRR  Creation                                  *
*                 9 Mar 2012  RRR  Optional frequency domain mask return     *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_ZLPFILT

    ZLPFILT

    Purpose: Performs low pass filtering by zeroing the FFT values.

    Syntax:  ZLPFILT(s, fc, mask)

                s - A series, the input series.

               fc - A real, the low pass cutoff frequency. Defaults
                    to rate(s) / 10.

             mask - Optional. An integer, return the frequency domain mask
                    only. Defaults to 0, return processed input.

    Returns:
             A series, the low pass filtered output.

    Example:
             W1: gsin(1000, 1/1000, 20) + gsin(1000, 1/1000, 1)
             W2: zlpfilt(w1, 10)

             W1 contains a 1 Hz sinewave with 20 Hz high frequency
             sinusoidal noise.

             W2 filters out the noise an returns the 1 Hz sinusoidal
             component.

    Remarks:
             ZLPFILT performs low pass filtering by zeroing out the
             high frequency components of the FFT of the input.

             FFT zeroing is a crude but fast filtering operation.
             However, the pass band and transition edge of the
             equivalent filter will exhibit significant ripple. For a
             flatter response, use traditional filter design methods
             such as those offered by DADiSP/Filters.

             If fc < 0, no filtering is performed.

             If mask == 1, the frequency domain binary mask used to zero
             out the stopband frequencies is returned. 

             See ZHPFILT for a high pass filter.

  See Also:
             DADiSP/Filters
             Zhpfilt
#endif


/* FFT zeroing low pass filter */
ITERATE zlpfilt(s, fc, mask)
{
        local len, bw, f, q;

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

                        fc = rate(s) / 10;
                }

                /* default to process input */
                mask = 0;
        }

        len = length(s);

        /* bin frequency width */
        bw  = rate(s) / len;

        /* build low pass frequency mask */
        f = zeros(len, 1, bw);

        if (fc < 0)
        {
                /* no filtering, mask is inverted here */
                return(mask ? not(f) : s);
        }

        /* index of samples to retain */
        q = int(fc / bw) + 1;

        if (q > len / 2)
        {
                q = len;
        }

        /* DC range passband */
        f[1..q] = 1;

        /* nyquist range passband */
        if (q >= 2 && q < len)
        {
                f[(len - q + 2)..len] = 1;
        }

        if (mask)
        {
                /* low pass mask */
                return(f);
        }
        else
        {
                /* frequency domain processing */
                f = ifft(fft(s) * f);

                if (isreal(s))
                {
                        f = real(f);
                }

                setxoffset(f, xoffset(s));

                return(f);
        }
}