{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}