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);
}