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


#if @HELP_ZHPFILT

    ZHPFILT

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

    Syntax:  ZHPFILT(s, fc, mask)

                 s - A series, the input series.

                fc - Optional. A real, the high 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 high pass filtered output.

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

             W1 contains a 20 Hertz sinewave riding on top of a 1 Hertz
             sinusoidal drift.

             W2 filters out the drift an returns the 20 Hertz sinusoidal
             component.

    Remarks:
             ZHPFILT performs high pass filtering by zeroing out the
             low 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 ZLPFILT for a low pass filter.

  See Also:
             DADiSP/Filters
             Zlpfilt
#endif


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

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

                        fc = rate(s) / 10;
                }

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

        len = length(s);

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

        if (mask)
        {
                /* unity frequency domain mask */
                f = ones(len, 1, bw);
        }

        /* no filtering */
        if (fc < 0)
        {
                return(mask ? f : s);
        }

        /* index of samples to zero, note that fc is filtered out */
        q = int(fc / bw) + 1;

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

        if (not(mask))
        {
                /* frequency domain */
                f = fft(s);
        }

        /* zero out DC range */
        f[1..q] = 0;

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

        if (mask)
        {
                /* return mask */
                return(f);
        }
        else
        {
                /* return time domain series */
                f = ifft(f);

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

                setxoffset(f, xoffset(s));

                return(f);
        }
}