View Raw SPL
/*****************************************************************************
*                                                                            *
*   FPADFILT.SPL     Copyright (C) 2009 DSP Development Corporation          *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    FIR filtering with optional padding using the FFT           *
*                                                                            *
*   Revisions:   22 Dec 2009  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#define ISEVEN(x) ((x)%2==0)

#if @HELP_FPADFILT

    FPADFILT

    Purpose: FIR filtering with optional endpoint padding using the FFT

    Syntax:  FPADFILT(s, h, pad, padlen)

                   s - input series

                   h - FIR filter coefficients (i.e. impulse response)

                 pad - optional integer, start and end point padding method,
                       0:none, 1:flip about end points, 2: flip about 0.0,
                       defaults to 0, no padding

              padlen - optional integer, length of segment to pad with,
                       defaults to length(h) / 2

    Returns: A series


    Example:
             W1: Integ(gnorm(1000,.001))*10+1
             W2: Kwlpass(1000.0, 50.0, 60.0, 70.0)
             W3: Fpadfilt(w1, w2, 0);
             W4: Fpadfilt(w1, w2, 1); overp(w3, lred); overp(w1, lgreen)

             W1 consists of simulated data with a low frequency trend.
             W2 contains an FIR lowpass filter with a cutoff frequency
             of 50 Hertz.  The data is filtered without padding in W3.

             Before filtering, the beginning and end points of the data
             are padded by flipping the end segments in time and
             amplitude in W4.  The overplot of the original and
             convolution data provides a comparison of the effects of
             the padding on the ramp-up and ramp-down transient
             implicit in the filtering process.  Standard filtering
             assumes the data is zero prior to the start and after the
             end of the data.

             W5: Fpadfilt(w1, w2, 2)

             Same as W4, except the start and end points are padded by
             flipping the start and end point segments about 0.0.

   Remarks:
             FPADFILT uses frequency domain convolution to perform the
             filtering.

             FPADFILT automatically sets the XOFFSET of the resulting
             data to the correct value.  For non-causal filters of
             even length, the XOFFSET of the output may be larger
             (i.e. more positive) than the XOFFSET of the input data.

             See PADFILT for a time domain implementation.

    See Also:
             Conv
             Endflip
             Fir
             Filteq
             Padfilt
             Zeroflip
#endif


/* FFT based PADFILT */
fpadfilt(s, h, pad, padlen)
{
        local yout;

        /* default pad mode and length */
        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                error("fpadfilt - input series and FIR filter required");
                        }
                        
                        pad = 0;
                }
                
                padlen = ceil(length(h) / 2) + ISEVEN(length(h));
        }

        /* use FFT to pad filter */
        yout = padfilt(s, h, pad, padlen, "fft");

        return(yout);
}