{Function: _filter2pole -- 2-Pole Lowpass/Highpass IIR filters
by Alex Matulich, August 2004
Copyright (c) 2004- by Unicorn Research Corporation and Alex Matulich.
All rights reserved.
This function implements either a Butterworth, Critically-damped, or
Bessel filter, as either low-pass or high-pass. All filters share the
same cutoff frequency, and this cutoff frequency is preserved when
cascading multiple filters of the same type.
See http://unicorn.us.com/trading/2polefilters.html for the calculation
steps, as well as plotted curves for frequency response and temporal
response to a step function, for all the filter types.
FUNCTION ARGUMENTS
------------------
p: This is the input series to filter.
length:
The length parameter is the 3 dB cutoff wavelength. Other lowpass
filters like the exponential moving average (known as XAverage, EMA, or
XMA) and the T3 Average have a cutoff wavelength at pi*length. To give
this filter the same cutoff wavelength as the EMA, you must first
multiply your length parameter by pi (3.14159265359) before applying
this filter. In practice, you'll never really need to do this.
type:
The default type is "bessel_lowpass" if you pass an unrecognized
string to this parameter. The recognized filter types are:
"butter_lowpass" = Butterworth lowpass
"butter_highpass" = Butterworth highpass
"critical_lowpass" = Critically-damped lowpass
"critical_highpass" = Critically-damped highpass
"bessel_lowpass" = Bessel lowpass
"bessel_highpass" = Bessel highpass
passes:
The passes parameter sets the number of filters to cascade to get the
final result. IIR filters like these are extremely fast to calculate,
so cascading multiple filters shouldn't adversely affect execution
time. You can, of course, cascade the filters externally yourself, but
doing it via the passes parameter will result in an internal adjustment
that preserves the 3 dB cutoff at the wavelength you set in the
length parameter.
NOTES
-----
Recommendation: For general filtering, use type bessel_lowpass, with
passes=2. A Bessel filter has the most consistent phase delay, and
two passes will result in better noise rejection than the T3 Average,
but without overshoot, unlike the T3. All of the filters implemented
here also have far superior noise rejection compared to an EMA having
the same cutoff wavelength.
The high pass filters, even critically-damped, will have oscillations
at passes >= 2. One can expect 2*passes-1 zero-crossings, actually.
To create a bandpass filter, simply cascade a lowpass and a highpass.
For example, a bandpass that passes all wavelengths between 10 and 20
would need a highpass filter of length 20 (passing frequencies > 1/20)
and a lowpass filter of length 10 (passing frequencies < 1/10). It
would be implemented as
bp = _filter2pole(_filter2pole(price, 20, "critical_highpass", 2),
10, "critical_lowpass", 2);
The bulk of the algorithm consists of initializations which are
executed only once. The filter itself is very short and fast.
Execution time slows down at low values of length because the sampling
rate must be increased, to avoid a ringing or numerical overflow when
the adjusted angular cutoff frequency Wc > tan(pi/4) for the lowpass
filter or when Wc < tan(pi/4) for the highpass filter. A sample rate
that maintains 0 < Wc < tan(pi/8) for the lowpass, or
tan(pi*3/8) < Wc < tan(pi/2) for the highpass, avoids this problem.
}
Inputs:
p(NumericSeries), {input series}
length(NumericSimple), {3 dB cutoff wavelength)}
type(String), {filter type (see comments above)}
passes(NumericSimple); {number of cascaded filters (MAX 5)}
vars:
sig0(0), sig1(0), sig2(0), {oversampled signal history 0,1,2}
a0(0),a1(0),a2(0), b1(0),b2(0), {filter coefficients}
j(0), k(0), m(0), {counters}
filtertype(0), type_butter(1), type_critical(2), type_bessel(3),
lowpass(True), {flag for low pass filter}
cascades(0), {number of passes corrected for bad user input}
Fs(0), {sampling frequency (normally 1.0 for large length)}
len(0), aFc(0), {adjusted filter length and adjusted cutoff freq}
d0(0), d1(0), {intermediate coefficients}
cc(0), {correction to preserve 3 dB cutoff frequency}
Wc(0), {angular cutoff frequency}
K1(0), K2(0); {K1,K2 = intermediate calculation values}
array: bw0[4](0), bw1[4](0), bw2[4](0); {filter history (0,1,2) for each pass}
{========== INITIALIZATION ==========}
if currentbar < passes+2 then begin
if type = "butterworth_highpass" then begin
filtertype = type_butter;
lowpass = false;
end else if type = "butterworth_lowpass" then begin
filtertype = type_butter;
lowpass = true;
end else if type = "critical_highpass" then begin
filtertype = type_critical;
lowpass = false;
end else if type = "critical_lowpass" then begin
filtertype = type_critical;
lowpass = true;
end else if type = "bessel_highpass" then begin
filtertype = type_bessel;
lowpass = false;
end else begin
filtertype = type_bessel;
lowpass = true;
end;
{limit bad user input}
cascades = passes;
if cascades < 1 then cascades = 1;
if cascades > 5 then cascades = 5;
len = length;
if len < 1.0001 then len = 1.0001;
{initialize cutoff correction and intermediate coefficients}
if filtertype = type_butter then begin
cc = Power(Power(2, 1/cascades) - 1, 0.25);
d0 = 1; d1 = SquareRoot(2);
end else if filtertype = type_critical then begin
cc = SquareRoot(Power(2, 0.5/cascades) - 1);
d0 = 1; d1 = 2;
end else begin
cc = SquareRoot(3)*SquareRoot(SquareRoot(power(2,1/passes)-0.75)-0.5);
d0 = 1/3; d1 = 1;
end;
{determine sampling rate Fs and adjust len and aFc accordingly}
if lowpass then begin
cc = 1/cc;
{limit corrected angular freq to 0 < Wc < tan(pi/8)}
Fs = IntPortion((8*cc/len) + 1);
len = len * Fs; {new length to use for resampled data}
aFc = cc / len; {adjusted lowpass cutoff freq}
end else {highpass} begin
{limit corrected angular freq to tan(pi*3/8) < Wc < tan(pi/2)}
Fs = IntPortion((8*cc/len) + 1);
len = len * Fs; {new length to use for resampled data}
aFc = 0.5 - cc/len; {adjusted cutoff to use in lowpass calc}
{normally aFc = Fs/2 - cc/len, but because we are outputting
a value every Fs samples rather than every sample, we use Fs=1
in this formula.}
end;
{calculate filter coefficients}
Wc = tangent(180 * aFc);
K1 = d1 * Wc;
K2 = Wc * Wc;
a0 = K2 / (d0 + K1 + K2);
a1 = 2 * a0;
a2 = a0;
b1 = a1 * (d0/K2 - 1);
b2 = 1 - (a0 + a1 + a2 + b1);
if lowpass = false then begin {adjust coefficients for highpass}
a1 = -a1;
b1 = -b1;
end;
{initialize first filter value}
if lowpass then value1 = p else value1 = 0;
for j = cascades-1 downto 0 begin
bw0[j] = value1; bw1[j] = value1; bw2[j] = value1;
end;
sig1 = p; sig0 = p;
end {END of initialization}
else begin {actual filtering starts here}
{========== N-PASS FILTER ==========}
for m = 1 to Fs begin {Fs=1 usually, except for small length}
j = cascades-1;
sig2 = sig1; {prior previous value of input signal}
sig1 = sig0; {previous value of input signal}
sig0 = p; {current input value}
{this is the 1-pass filter}
bw2[j] = bw1[j]; {prior previous filter value}
bw1[j] = bw0[j]; {previous filter value}
bw0[j] = a0*sig0 + a1*sig1 + a2*sig2 + b1*bw1[j] + b2*bw2[j];
{cascade additional filters if needed}
while j > 0 begin
k = j; j = j - 1;
bw2[j] = bw1[j]; {shift current filter values to prior values}
bw1[j] = bw0[j];
bw0[j] = a0*bw0[k] + a1*bw1[k] + a2*bw2[k] + b1*bw1[j] + b2*bw2[j];
end;
end;
end;
_filter2pole = bw0[0]; {final filter value for this bar}