Generalized IIR Savitzky-Golay Filters: Derivation, Parameterization, Realization and an Adaptive Configuration Hugh L. Kennedy
Abstract—Finite-impulse-response (FIR) Savitzky-Golay (SG) filters, with a finite memory, have been used for decades as datasmoothers for signals that are locally well-approximated by polynomials in the time domain. Analysis of their frequency response also reveals that low group-delay configurations have highly desirable properties, such as good magnitude flatness and phase linearity at low frequencies, which makes them suitable as general-purpose low-pass filters. Their less well-known infinite-impulse-response (IIR) counterparts, with a fading memory, share these properties and it is shown it this paper that they are indeed superior in many respects. They are extended and improved here by generalizing the error weighting function used when deriving their coefficients via discounted least-squares analysis in the time domain. The resulting filters are readily configured using just a few parameters that are linked to the regression-based design procedure. As a useful byproduct, the filters generate the associated Laguerre spectrum, which is used to recursively estimate the variance of the reconstructed signal. An adaptive filter-bank arrangement with a variance-based switching mechanism is also presented, which has a lowgroup delay and a greater bandwidth when levels of ambient noise and interference are low; a longer group delay is applied for greater high-frequency attenuation when required.
Keywords— Digital filters, Estimation theory, IIR filters, Low pass filters, Orthogonal transforms, Polynomials, Recursive filters, Zero assignment.
I. INTRODUCTION
S
avitzky-Golay (SG) filters, with a finite impulse response (FIR), are traditionally used as smoothing filters in dataprocessing applications where local polynomial models are known to be a good representation of the underlying process [1],[2]; however, as awareness of their excellent low-frequency phase and magnitude response grows
[3],[4],[5],[6], they are increasingly being used as general-purpose low-pass filters where it is necessary to track “any reasonably smooth process” [10]. The universality of this problem (feedback control, sensor systems, industrial and consumer electronics etc.) and the suitability/flexibility of these filters, has inspired recent research on their properties in more generic digital-signal- processing (DSP) contexts [7],[8]. SG filters with an infinite impulse response (IIR), or fading-memory polynomial (FMP) filters, have been used for many
H. L. Kennedy is with the School of Engineering, University of South Australia, Adelaide, SA 5095, Australia (e-mail:
[email protected]).
years in data processing applications, (e.g. target tracking) [9],[10],[11]; however their potential utility as low-order lowpass filters in broader application areas has been largely overlooked. One of the purposes of this paper is to compare and contrast the FIR- and IIR-SG filters (section II), to establish a common and general framework for their design (sections III & IV) and to elucidate the benefits of IIR designs (sections IV,V & VII). However, revisiting the foundations of FMP filter design [9], reveals that the full potential of IIR-SG filters has yet to be realized and exploited. Some novel enhancements are described here. Conventional FIR-SG filters are derived using least-squares regression with a uniform error weight and the discrete Legendre polynomials (or Gram polynomials, or Discrete Chebyshev polynomials) [9],[10],[12],[13],[14],[15]; ‘regular’ IIRSG filters are derived using discounted least-squares regression with an exponentially decaying error weight and the discrete Laguerre polynomials [9],[10],[11],[16],[17],[18]; whereas the generalized IIR-SG filters described here are obtained when the exponential weighting function is extended to include a shape parameter 𝜅 and used with the discrete associated Laguerre polynomials [19]. Low-order filters are appropriate in real-time systems that require simple components with a low-complexity. Closedform expressions for the coefficients of low-order IIR-SG filters, as a function of the pole location (𝑝) and the nominal group delay (𝑞), for some specific combinations of model degree orders (𝐵 = 1 or 2) and weighting function forms (𝜅 = 0 or 1), have previously been tabulated by the author and applied to problems in image processing [18] and [19] and feedback control [16],[17]. This paper however, provides a more detailed description of a general process that may be used to derive the coefficients of higher-order filters for arbitrary parameter combinations (𝐵 ≥ 0 and 𝜅 ≥ 0). The generalized IIR-SG filters have repeated poles on the real axis between 0 and 1 (at 𝑝). In this paper it is shown that additional filter poles may be used either to improve low-frequency magnitude flatness and phase linearity (i.e. in the ‘passband’) or to improve mid- to high-frequency attenuation (i.e. in the ‘stopband’). Deriving the filter coefficients from within a regression framework, ensures that the zeroes are appropriately placed in the complex 𝑧-plane for the desired effect. IIR-SG filters are efficient and flexible; although unfortunately, compromises and trade-offs between conflicting requirements such as: frequency selectivity, temporal selectivity/responsiveness, group delay and circuit complexity, are invariably required in any given design problem. This dilemma motivated the derivation of the design guidelines discussed in sections IV and VII and the adaptive filter configuration presented in section VI; these sections rely on the use of the variance reduction factor (VRF) to support offline and online filter analysis. II. OVERVIEW OF FIR- AND IIR-SG FILTERING The analysis of sensor observations using polynomials has played an important role in the development of modern science for more than two centuries [5],[10],[20]. Savitzky and Golay showed that polynomial smoothing (i.e. low-pass filtering) may be implemented using very few logical/numeric operations/instructions using linear-phase filters with an FIR and applied to (off-line data-processing) problems in physical chemistry [1]. Interpolative and extrapolative extensions quickly followed [2], which permit low-frequency phase/delay to be manipulated [12],[16],[17], to customize performance in systems that require causal filters (e.g. for on-line signal-processing). It was then shown that efficient recursive realizations of all the aforementioned FIR-SG filters exist, for the integration of very low-complexity circuits [12]. While the basic linear-phase FIR form is commonly attributed to Savitzky and Golay, the extensive list of references in Morrison’s early monograph on prediction and smoothing suggests that the polynomial filters were already well understood during
the early 1960s, within Bell laboratories at least, and being considered for aerospace/defense applications [9]. Regardless, interest in these filters as simple components in digital systems continues to grow, and their desirable properties, in both the time and frequency domains, are increasingly attracting attention [7],[8],[14],[15]. Unfortunately, one of the less well-known properties of the FIR-SG filters is that the recursive realization is susceptible to long-term rounding-error accumulation, in finite-precision machines, due to the utilization of pole-zero cancellation on the unit circle at 𝑧 = 1. If a pole on the unit circle results in a marginally-stable filter, yielding an impulse response with an envelope that neither grows nor decays over time; then a cancelled pole on the unit circle results in a ‘provisionally’-stable filter, i.e. an impulse response that only vanishes if various conditions are satisfied so that the pole and zero remain coincident. While this issue is summarized succinctly, although somewhat abstractly, in the complex 𝑧-domain; it does have practical ramifications that do need to be considered in real systems. Recursive realizations of the discrete Fourier transform (DFT) present similar problems in broad-band frequency analysis [21]. Fortunately, there is an alternative form of polynomial signal analysis with a recursive realization that is guaranteed to be stable, with all poles inside the unit circle, and an IIR. This filter is referred to here as an IIR-SG and is equivalent to the FMP filters considered by Morrison [9],[10]. The FIR-SG filter and the IIR-SG filter are both designed using linear regression analysis; however the former approach performs a least-squares fit using uniformly weighted samples over a finite time window; whereas the latter approach performs discounted least-squares regression using an exponentially decaying weight over a time window of infinite extent. The use of a finite memory, with the abrupt truncation of old data, in the FIR-SG filter means that the desired frequency response of the ideal filter is convolved with the undesirable frequency response of the rectangular time window (i.e. the Dirichlet kernel), which has a high side-lobe level. In addition to an accurate and efficient recursive structure, the use of a fading memory, with the gradual de-emphasis of old data, in the IIR-SG filter results in a frequency response with some very attractive properties, such as a monotonic magnitude decay (i.e. no side-lobes) outside the passband and excellent high-frequency noise attenuation. Furthermore, it has been shown that a null may be placed at the Nyquist frequency when the nominal low-frequency group-delay is appropriately set, for a zero at 𝑧 = −1 [19]. It is shown here that the use of this optimal delay also minimizes the expected estimation error for a given IIR-SG filter. For an FIR filter, the ideal delay in most cases, is of course at the midpoint of the finite time window, for a linear phase filter; this result should not be surprising, because at this point, there is an equal number of samples before and after the estimate. Deriving filter coefficients using linear regression, promotes a joint time-frequency perspective, which allows filter properties in both domains to be customized using just a handful of adjustable parameters. This feature is particularly helpful in control system design [16],[17], where the transient response (e.g. to step inputs) of a digital controller (i.e. the filter) and the closed-loop system are of primary importance, and where it is customary to fine-tune controllers online after system integration and deployment. IIR-SG or FMP filters, have been used for many years as target tracking filters, in radar data-processing applications, as an alternative to the Kalman filter [9],[10]. In this capacity, an emphasis is typically placed on time-domain performance metrics such as random errors (due to sensor noise) and systematic errors (due to model mismatch) in target state estimates. FIR-SG filters are also used in these roles; however, they are increasingly being used in DSP applications, which has led to greater appreciation of their desirable frequency-domain characteristics, such as their magnitude flatness and phase linearity in the near-dc region. This paper develops the early work described in [9] which has recently been updated in [10]. In both references the emphasis is on (radar) data processing rather than signal processing. This paper
complements these works by focusing on frequency domain design and analysis (via the 𝑧-domain) and some novel extensions are proposed. In the next section – indeed, in the remainder of this paper – it is shown that, relative to the FIR-SG filter, the IIR-SG filter has some highly desirable properties, and that the range of potential benefits and applications is further broadened when the error-weighting function is generalized using the shape parameter 𝜅. III. DERIVATION A. Signal Model and Least-Squares Estimation A continuous-time noise-corrupted input is sampled at time instants 𝑡 = 𝑛𝑇, where 𝑇 is the sampling period and 𝑛 is the sample index. The following model is used to represent the discretized signal, over a specified ‘time-scale’, in the vicinity of 𝑛: 𝑥(𝑛 − 𝑚) = ∑𝐵𝑘̅=0 𝛼𝑘̅ (𝑛)𝜑𝑘̅ (𝑚) = ∑𝐵𝑘=0 𝛽𝑘 (𝑛)𝜓𝑘 (𝑚) 𝑦(𝑛) = 𝑥(𝑛) + 𝜀
(1a) (1b) (1c)
where: 𝑦(𝑛) is the 𝑛th ‘noise-corrupted’ measurement, or sample; 𝑥(𝑛) is the corresponding ‘noise-free’ signal; 𝜀 is a Gaussian-distributed noise term, with 𝜀~𝒩(0, 𝜎𝑦2 ); 𝜑𝑘̅ is the 𝑘̅th discrete component basis function; 𝜓𝑘 is the 𝑘th discrete orthonormal basis function; 𝛼 are the component model coefficients; 𝛽 are the orthonormal model coefficients; 𝐵 is the model degree; and 𝑚 is a delay index (𝑚 ≥ 0 for causality). The model coefficient vector at the time of the 𝑛th sample, using all available measurements (previous and present), is estimated via weighted least-squares regression, which minimizes the weighted sum-of-squared residuals (SSRs) SSR(𝑛) = ∑∞ ̂(𝑛 − 𝑚)]2 𝑚=0 𝑤(𝑚)[𝑦(𝑛 − 𝑚) − 𝑥
(2)
where 𝑤(𝑚) is an error-weighting function. The weighting function need not be normalized and may be arbitrarily selected [13]; however, forms with simple 𝒵-transforms, that lead to recursive structures, are favored. It will later be shown that the choice is as important as basis function selection. Use of a finite-memory error-weighting function, e.g. a uniform weight: 1, 0 ≤ 𝑚 < 𝑀 𝑤(𝑚) = { 0, 𝑀 ≤ 𝑚 < ∞
(3a)
leads to FIR filters; whereas use of a fading-memory error-weighting function: 𝑤(𝑚) = 𝑚𝜅 𝑒 𝜎𝑚 = 𝑚𝜅 𝑝𝑚
(3b)
with 𝑝 = 𝑒 𝜎 and 𝜎 < 0 for a meaningful physical interpretation and a guaranteed-stable filter realization, leads to IIR filters. Generalized window functions of this form have been used to improve spectra (i.e. narrower main-lobe and lower sidelobes) produced by recursive short-time Fourier analyzers [21],[22],[23],[24]. ̅
Unless otherwise stated, it is assumed here that 𝜑𝑘̅ (𝑚) = 𝑚𝑘 , for a polynomial model, and that 𝜓𝑘 is a linear combination of 𝑘 + 1 components. It is also assumed here, that in practice, the use of polynomial filter models is not restricted to polynomial input signals. A polynomial model of a low degree is also an appropriate basis set for low-frequency bandlimited signals if the time-scale or memory of the regression analysis, as determined by 𝑤(𝑚), is short, relative to the ‘timescale’ of signal variation, as determined by the period of the signal’s highest-frequency component. B. Orthonormal Basis-Set Construction Design of FIR or IIR filters, to implement least-squares signal analysis, is greatly simplified if the basis set is orthonormal with respect to the error-weighting function, i.e. ∑∞ 𝑚=0 𝜓𝑘1 (𝑚)𝑤(𝑚)𝜓𝑘2 (𝑚) = 𝛿𝑘1 𝑘2
(4)
where 𝛿 is the Kronecker delta function. For instance, the discrete Legendre polynomials form an orthonormal set when the uniform weight specified in (3a) is applied over a finite discretized interval; whereas the discrete associated Laguerre polynomials form an orthonormal set when the weight in (3b) is applied over an infinite discretized interval. Note that the discrete Laguerre polynomials may be used in the latter case when an exponential weight (𝜅 = 0) is used, which leads to what are referred to here as the ‘regular’ IIR-SG filters [9],[10]. The orthonormal basis set 𝝍(𝑚) = [𝜓0 (𝑚), … , 𝜓𝑘 (𝑚), … , 𝜓𝐵 (𝑚)]T (indexed using 𝑘), is formed from the component basis functions 𝝋(𝑚) = [𝜑0 (𝑚), … , 𝜑𝑘̅ (𝑚), … , 𝜑𝐵 (𝑚)]T (indexed using 𝑘̅), where the T superscript (no italics) is a transpose operation, using the transformation 𝝍(𝑚) = 𝓐𝝋(𝑚). The ‘orthonormalizing’ matrix 𝓐, is a lower-triangular, 𝐵 + 1 by 𝐵 + 1 matrix, which may be found using the Gram–Schmidt procedure or (more conveniently) using a Cholesky decomposition of the 𝐵 + 1 by 𝐵 + 1 summation matrix 𝓑, with elements ℬ𝑘̅1𝑘̅2 = ∑∞ ̅ 1 (𝑚)𝑤(𝑚)𝜑𝑘 ̅ 2 (𝑚) 𝑚=0 𝜑𝑘
(5a)
which may be evaluated in the complex 𝑧-domain using ℬ𝑘̅1𝑘̅2 = 𝒵{𝜑𝑘̅1 (𝑚)𝑤(𝑚)𝜑𝑘̅2 (𝑚)}|
𝑧=1
.
(5b)
If the Cholesky decomposition of 𝓑 yields 𝓑 = 𝑳𝑳T , where 𝑳 is a (lower) triangular matrix, then 𝓐 = 𝑳−1 , such that 𝑰 = 𝓐𝓑𝓐T , where 𝑰 is the identity matrix. For a polynomial model with narrow-band low-pass characteristics:
T
̅
𝝋(𝑚) = [1, … , 𝑚𝑘 , … , 𝑚𝐵 ] ;
(6a)
alternatively, for a sinusoidal model with wide-band low-pass characteristics,
̅
̅
𝝋(𝑚) = [𝜑 −𝐵 (𝑚), … , 𝜑 −𝑘 (𝑚), … ,1, … , 𝜑 +𝑘 (𝑚), … , 𝜑 +𝐵 (𝑚)]
T
(6b)
where 𝜑(𝑚) = 𝑒 𝑗2𝜋𝑚⁄𝑀 for a bandwidth of 𝐵 ⁄𝑀 cycles per sample [21]. Note that in the sinusoidal case, the 𝓐 and 𝓑 matrices are 2𝐵 + 1 by 2𝐵 + 1. Polynomial and sinusoidal models lead to (FIR & IIR) SG filters [1],[9], and frequency sampling filters [21],[22],[23],[24], respectively; finite and fading memories yield FIR and IIR variants, respectively. Thus in the fading-memory polynomial case, ̅
̅
ℬ𝑘̅1𝑘̅2 = 𝒵{𝑚𝑘1+𝑘2 +𝜅 𝑝𝑚 }|𝑧=1 = 𝒬(𝑘̅1 + 𝑘̅2 + 𝜅; 𝑝).
(7)
The resulting coefficient 𝒬, may be simply evaluated, for any power of 𝑚, using Eulerian numbers and polylogarithms via the closed-form expressions provided in [25]. C. Parameter/Spectrum Estimation (Analysis) For a defined signal model and a specified weighting function, the use of an orthonormal basis set allows the model parameters to be estimated at the time of the 𝑛th sample, using all measurements up to and including 𝑦(𝑛), via the analysis equation 𝛽̂𝑘 (𝑛) = ∑∞ 𝑚=0 𝜓𝑘 (𝑚)𝑤(𝑚)𝑦(𝑛 − 𝑚).
(8)
̂ (𝑛) = [𝛽̂0 (𝑚), … , 𝛽̂𝑘 (𝑚), … , 𝛽̂𝐵 (𝑚)]T is the discrete: Fourier spectrum for sinusoidal components The coefficient vector 𝜷 and a finite memory; Legendre spectrum for polynomial components and a finite memory; and the associated Laguerre spectrum for a generalized fading memory. (Note that the “hat” accent is used in this paper to denote estimated quantities.) The summation in (8) is of course truncated after 𝑀 terms if a finite-memory analysis is performed. For fading-memory analysis, stable recursive implementations are reached via the complex 𝑧-domain using ̂ (𝑧) = 𝓐𝛘(𝑧)𝑌(𝑧) 𝔅
(9)
̂ (𝑧) = 𝒵{𝛽̂𝑘 (𝑛)}, 𝑌(𝑧) = 𝒵{𝑦(𝑛)} and 𝛘(𝑧) = 𝒵{𝝋(𝑚)𝑤(𝑚)}. The column vector 𝛘(𝑧), has with 𝒵-transformed quantities 𝔅 ̅
a length of 𝐵 + 1, with elements χ𝑘̅ (𝑧) = 𝒵{𝑚𝜅+𝑘 𝑝𝑚 } which are ratios of polynomials in 𝑧 i.e. 𝐵χ (𝑧)⁄𝐴χ (𝑧) with 𝜅 + 𝑘̅ + 1 repeated (real) poles. The 𝑧-domain analysis operation in (9) may therefore be further simplified by a partial-fraction expansion of each element of 𝛘(𝑧) such that 𝛘(𝑧) = 𝓒𝓕(𝑧) where 𝓒 is a 𝐵 + 1 by 𝐵 + 𝜅 + 1 coefficient matrix. The column vector 𝓕(𝑧), has a length of 𝐵 + 𝜅 + 1, with elements ℐ𝑘̿ (𝑧) = [𝐻𝑝 (𝑧)]
̿ +1 𝑘
̿
= {1⁄(1 − 𝑝𝑧 −1 )}𝑘+1 (note the use of a double-bar accent, to denote
the index of the fundamental analysis unit, see Fig. 8c), thus (9) becomes ̂ (𝑧) = 𝓐𝓒𝓕(𝑧)𝑌(𝑧). 𝔅
(10)
This derivation began with the orthonormal model in (1b), which is amenable to implementation; however, the component model in (1a) may be more amenable to analysis and interpretation. The spectrum in orthonormal coordinates is accordingly transformed back into component coordinates using ̂ (𝑛) ̂ (𝑛) = 𝓐T 𝜷 𝜶
(11a)
with ̂ (𝑛) = [𝛼̂0 (𝑚), … , 𝛼̂𝑘̅ (𝑚), … , 𝛼̂𝐵 (𝑚)]T 𝜶
(11b)
such that ̅
𝑥̂(𝑛 − 𝑚) = ∑𝐵𝑘̅=0 𝛼𝑘̅ (𝑛)𝜑𝑘̅ (𝑚) = ∑𝐵𝑘̅=0 𝛼𝑘̅ (𝑛)𝑚𝑘
(11c)
or ̂ (𝑧) ̂ (𝑧) = 𝓐T 𝔅 𝕬
(11d)
where ̂ (𝑧) = 𝒵{𝜶 ̂ (𝑛)}. 𝕬
(11e)
Substitution of (10) into (11d) yields ̂ (𝑧) = 𝓐T 𝓐𝓒𝓕(𝑧)𝑌(𝑧) 𝕬
(12)
which is simply a bank of 𝐵 + 𝜅 + 1 cascading first-order filters [26], with outputs mixed using the 𝓐T 𝓐𝓒 matrix (see Fig. 8a). The spectrum coefficients resulting from this analysis operation are often useful in their own right [18],[27]; for instance, they may be used for coding, compression, detection, matched filtering, optical-flow determination, transform-domain adaptive filtering, and as will be demonstrated in subsection III.E – residual analysis and variance estimation; however, the standard (non-adaptive) FIR- and IIR-SG filters are reached by performing signal reconstruction or a synthesis operation, for a low-pass filtering effect. D. Signal Reconstruction (Synthesis) With the spectrum constructed, the estimate 𝑥̂(𝑛), of the signal 𝑥(𝑛), derived using noisy measurements 𝑦(𝑛 − 𝑚) for 𝑚 = 0 … ∞ , may be evaluated at time 𝑇[𝑛 − 𝑚], using the estimated model parameters 𝛽̂𝑘 (𝑛) in the synthesis equation, which is simply obtained by substituting 𝛽̂𝑘 (𝑛) for 𝛽𝑘 (𝑛) in (1b), i.e. 𝑥̂(𝑛 − 𝑚) = ∑𝐵𝑘=0 𝛽̂𝑘 (𝑛)𝜓𝑘 (𝑚).
(13a)
For reconstruction at a fixed offset of 𝑞 samples from the current sample at 𝑚 = 0:
𝑥̂𝑞 (𝑛) = 𝑥̂(𝑛 − 𝑞) = ∑𝐵𝑘=0 𝛽̂𝑘 (𝑛)𝜓𝑘 (𝑚)|𝑚=𝑞 =
(13a)
∑𝐵𝑘̅=0 𝛼̂𝑘̅ (𝑛)𝜑𝑘̅ (𝑚)|𝑚=𝑞
(13b)
(𝑛)𝑚𝑘̅
(13c)
= ∑𝐵𝑘̅=0 𝛼̂𝑘̅
|𝑚=𝑞
̅
= ∑𝐵𝑘̅=0 𝛼̂𝑘̅ (𝑛)𝑞 𝑘 .
(13d)
Or in the 𝑧 domain, with the assistance of (12): ̂ (𝑧) 𝑋̂𝑞 (𝑧) = 𝒒T 𝕬
(14a)
where 𝑋̂𝑞 (𝑧) = 𝒵{𝑥̂𝑞 (𝑛)}
(14b)
and ̅
T
𝒒 = [1, … , 𝑞 𝑘 , … , 𝑞 𝐵 ] .
(14c)
(Note that only smoothers are considered here; differentiators are derived by evaluating the 𝐷th derivative of 𝑥̂ at 𝑞.) Substitution of (12) into (14a), and dividing both sides by 𝑌(𝑧), yields the following discrete-time transfer function of the generalized IIR-SG filter 𝐻𝑥̂ (𝑧), as ratios of polynomials in 𝑧 i.e. 𝐵𝑥̂ (𝑧)⁄𝐴𝑥̂ (𝑧), with analysis and synthesis stages combined: 𝐻𝑥̂ (𝑧) = 𝑋̂𝑞 (𝑧)⁄𝑌(𝑧) = 𝒒T 𝓐T 𝓐𝓒𝓕(𝑧)
(15a)
which may be expanded as 𝑀 −1
𝑀 −1
𝑏 𝑎 𝐻𝑥̂ (𝑧) = 𝐵𝑥̂ (𝑧)⁄𝐴𝑥̂ (𝑧) = ∑𝑚=0 𝑏𝑚 𝑧 −1 ⁄∑𝑚=0 𝑎𝑚 𝑧 −1 (15b)
which, after taking inverse 𝒵-transforms, is realized in the time domain using the linear difference equation 𝑀 −1
𝑀 −1
𝑏 𝑎 𝑥̂𝑞 (𝑛) = ∑𝑚=0 𝑏𝑚 𝑦(𝑛 − 𝑚) − ∑𝑚=1 𝑎𝑚 𝑥̂𝑞 (𝑛 − 𝑚)
(15c)
where 𝑀𝑏 = 𝐵 + 𝜅 + 1 and 𝑀𝑎 = 𝐵 + 𝜅 + 2; the filter order is equal to 𝑀𝑎 − 1. Note that 𝐵𝑥̂ (𝑧) is a function of 𝑝 and 𝑞; whereas 𝐴𝑥̂ (𝑧) is a function of 𝑝 only. The frequency response 𝐻𝑥̂ (𝜔), of the filter, as a function of the angular frequency 𝜔 = 2𝜋𝑓 (radians/sample), where 𝑓 is the relative frequency (cycles/sample), is determined using 𝐻𝑥̂ (𝜔) = 𝐻𝑥̂ (𝑧)|𝑧=𝑒 𝑗𝜔 . The SG filters are maximally flat at 𝜔 = 0 thus 𝑑 𝐷 |𝐻𝑥̂ (𝜔)|2 ⁄𝑑𝜔𝐷 |𝜔=0 = 0, for 𝐷 = 0 … 𝐷max , where 𝐷max increases with 𝐵. E. Analysis of Variance Expansion of (2) followed by the substitution of (13a) for 𝑥̂(𝑛 − 𝑚) and simplifying, with the assistance of the orthonormality condition stated in (4), yields
𝐵 2 ̂2 SSR(𝑛) = ∑∞ 𝑚=0 𝑤(𝑚)𝑦 (𝑛 − 𝑚) − ∑𝑘=0 𝛽𝑘 (𝑛).
(16)
Note that the SSR, as it is defined in (4) and (16), involves squared differences of 𝑥̂ and 𝑦 over all currently available samples (i.e. past and present), using the latest parameter estimates derived from all currently available samples. Like the other quantities of interest computed so far, recursive difference equations for the SSR are obtained via the complex 𝑧-domain (see Fig. 8d). Furthermore, let 𝜎𝑦2 = E〈(𝑦 − 𝑥)2 〉
(17a)
and 2
𝜎𝑥2̂ = E 〈(𝑥̂𝑞 − 𝑥𝑞 ) 〉
(17b)
where E〈∙〉 is the expectation operator. While the variance of the sensor-noise 𝜎𝑦2 , remains constant, the variance of the estimate 𝜎𝑥2̂ , reduces to a constant value, which is determined by the 𝑝 & 𝑞 parameters. In the latter case, i.e. for (17b), this constant limit is found using the variance reduction factor (VRF, a constant) as follows [9],[10]: 𝜎𝑥2̂ = VRF𝜎𝑦2 .
(18)
In the absence of the prior knowledge regarding the 𝜎𝑦2 and 𝜎𝑥2̂ parameters, 𝜎̂𝑦2 (𝑛) = 𝑐𝑝 ∙ SSR(𝑛)
(19a)
and 𝜎̂𝑥2̂ (𝑛) = 𝑐𝑝 ∙ VRF ∙ SSR(𝑛)
(19b)
where 𝑐𝑝 is normalizing factor with 𝑐𝑝 = 1⁄𝒬(𝜅; 𝑝), such that 𝑐𝑝 ∑∞ 𝑚=0 𝑤(𝑚) = 1. The statistical properties of these variance estimates are not considered here, as such an analysis typically assumes that the model stated in (1) actually holds; however, in this work (8) and (13a) are deliberately used in cases where this is known/expected to not necessarily be the case (e.g. sinusoidal signals, non-Gaussian noise, interference, time varying parameters, etc.). In these instances, gross errors in (19) are to be expected. Adapting the definition of the VRF provided elsewhere – see (13.6.12) in [9] – to suit the approach adopted here (i.e. generalized weight, smoothing only with no derivatives, utilization of the 𝑧 domain, a preference for matrix/vector notation, different symbols, etc.) yields VRF = 𝝓T (𝑞)𝓐T 𝓐𝓖𝓐T 𝓐𝝓(𝑞).
(20)
The elements of 𝓖, a 𝐵 + 1 by 𝐵 + 1 coefficient matrix, are the infinite summations ̅
̅
2𝑚 𝑘1 +𝑘2 +2𝜅 𝒢𝑘̅1𝑘̅2 = ∑∞ 𝑚 = 𝒬(𝑘̅1 + 𝑘̅2 + 2𝜅; 𝑝2 ). (21) 𝑚=0 𝑝
IV. PARAMETERIZATION A. The Variance Reduction Factor (𝑉𝑅𝐹) The VRF is the white noise gain for an SG filter with 𝐵 set to match the degree of a given polynomial signal, thus it provides an indication of signal estimation accuracy under ideal circumstances when the model assumed in (1) is valid. However, a frequency response plot is a more meaningful representation of an SG filter’s properties and robustness in the general case where the noise is a sum of high-frequency sinusoids and the signal is a sum of low-frequency sinusoids. The VRF thus provides an indication of the ‘frequency selectivity’ of an SG filter – the passband width decreases, and the stopband frequency attenuation increases, as the VRF decreases. The VRF is a useful figure-of-merit of the low-pass filtering ability for a given SG filter, or its ability to track “any reasonably smooth process” [10]; it is computed offline for filter analysis and design purposes and/or used online to compute 𝜎𝑥2̂ from the SSR – see (19b). Closed-from expressions for the VRF of FIR-filters and ‘regular’ IIR-SG filters, i.e. with 𝜅 = 0, are provided in [9]; for the ‘generalized’ IIR-SG filters, i.e. with 𝜅 ≥ 0, (20) should instead be used. The VRF of an FIR-SG filter, for a given 𝐵, is a function of the 𝑀 and 𝑞 parameters; whereas the VRF of a generalized IIR-SG filter, for a given 𝐵 and 𝜅 combination, is a function of the 𝑝 and 𝑞 parameters; FIR- and IIR-SG filters with the same VRF are said to have the same “effective memory length” [10]. Thus if the VRF is specified, it is possible to derive the coefficients of an IIR-SG filter either by setting 𝑞 and solving for 𝑝 (as done in Figs. 1, 2 & 6) or by setting 𝑝 and solving for 𝑞. Alternatively, if the derivative of the VRF with respect to 𝑞 is determined and set to zero, it is possible to derive the ‘optimal’ value of 𝑞 that minimizes the VRF, for a specified value of 𝑝 (as done in Figs. 3-5); for multiple minima, the solution that is closest to 𝑞 = 0 is chosen. These solutions place a filter zero at 𝑧 = −1, for a Nyquist null, which results in infinite attenuation (i.e. zero transmission gain) at 𝑓 = 1⁄2. B. The Pole Position Parameter (𝑝) As 𝑤(𝑚) = 𝑚𝜅 𝑝𝑚 , this parameter ultimately determines the impulse response duration, or the length of its ‘tail’ in the general case for 𝜅 > 0. For a given IIR-SG filter, the VRF decreases as 𝑝 is increased from zero and approaches unity; consequently, the impact of random noise decreases and the frequency selectivity increases. These generally desirable behavioral changes are a consequence of a longer filter memory which means that more data are considered when forming the signal estimate. However, shifting 𝑝 to the right along the real axis also: 1) increases systematic (or bias) errors in cases of persistent model mismatch, 2) increases the phase lag at mid frequencies, which may be destabilizing in closed-loop systems, and 3) increases the time required for the filter to accommodate signal parameters changes or to ‘age out’ transient artefacts due to parametric discontinuities (e.g. ‘maneuvers’) [9],[10]. The determination of 𝑝 in an IIR-SG filter is analogous to the selection of 𝑀 in an FIR-SG filter [8]. This “bias-variance tradeoff” is considered at depth in [8],[9] and [10]; however, the issue of bias invariably requires a specification of alternative forms that the input signal might take. In [8] and [10] it is assumed that bias is due to a polynomial model of insufficient order; however, here it assumed that inputs have a sinusoidal form. Taking a frequency-domain perspective – the aforementioned tradeoff may be restated as passband width vs stopband attenuation. Decisions regarding the pole position are best made at the start of the design process, as the form of the orthogonal basis set (as specified by 𝓐), the structure of the analyzers (as defined by 𝓒 and 𝓕) and the estimation accuracy limits (as computed using 𝓖) all depend on 𝑝. Inspection of (15a) and (20) suggests that it is much easier to vary the zero locations, for a given 𝑝 value, because powers of the 𝑞 parameter only appear as pre/post multipliers which collect elements of the aforementioned terms.
C. The (Nominal) Group-Delay Parameter (𝑞) Inspection of (13) suggests that using 𝑞 < 0 corresponds to polynomial prediction (for a low-frequency phase-lead); whereas using 𝑞 > 0 corresponds to polynomial interpolation for non-integer 𝑞 (for a low-frequency phase-lag) [16],[17]. In all cases, some low-frequency magnitude flatness and phase linearity is achieved; although in general, at mid-to-high frequencies, attenuation turns to amplification for 𝑞 ≪ 0, i.e. long-term prediction amplifies noise, especially when short ‘data records’ are used (with 𝑝 → 0). As a ‘rule of thumb’, for a low-frequency estimate of reasonable quality (as indicated by a low VRF) 𝑞 should be close to the centroid of the weighting function to ensure that the estimate is well ‘supported’ by the data. Thus the 𝑞 parameter determines the nominal group delay, where the “nominal” qualifier is used here to highlight the fact that the SG filters (FIR & IIR) are only approximately linear phase in the low-frequency region and that only the FIR-SG filters with 𝑞 = 12[𝑀 − 1] are exactly linear phase – over the entire frequency domain. If the group delay is not a strict requirement, then the 𝑞 parameter may be set to minimize the VRF (see Fig. 3) for a Nyquist null (see Fig. 4), which results in exceptional high-frequency attenuation (see Fig. 4). D. The Polynomial Model Degree (𝐵) Increasing 𝐵 increases the passband width (see Fig. 5) but it also increases the minimum VRF (see Fig. 3) and the corresponding stopband attenuation (see Fig. 4), because it introduces more degrees of freedom in the regression analysis through which noise or out-of-band interference may be accommodated. The passband is defined as the region over which the gain is unity and the phase shift is linear, to within an acceptable tolerance (e.g. 3 dB and/or 10°). Note that for the SG filters, both the magnitude and phase response are perfect at 𝑓 = 0 then deteriorate at approximately the same rate with increasing frequency, for a given combination of parameters. E. The Weighting-Function Shape-Parameter (𝜅) Use of the generalized weighting function introduces an extra degree of design freedom. If 𝑞 is fixed then 𝑝 may be varied to yield a specified VRF for a given pairing of 𝐵 & 𝜅. Incrementing 𝜅 while reducing 𝑝 to keep the VRF constant, increases the passband width and the stopband attenuation (see Figs. 1, 2 & 6). However if 𝑝 is fixed and 𝑞 is set to minimize the VRF for a given pairing of 𝐵 & 𝜅, then incrementing 𝜅 increases the frequency selectivity of the filter by decreasing the passband width and dramatically increasing the stopband attenuation (see Figs. 3-5). Increasing 𝜅 ‘rebalances’ the error weight to emphasize older samples, which provides better ‘support’ for signal reconstruction at greater delays, and a reduced VRF, by using an increased value of 𝑞 (see Figs. 3 & 4). Examination of Fig. 3 reveals the main benefit of using generalized IIR-SG filters: For a given number of filter poles (𝐵 + 𝜅 + 1), a lower VRF is attainable, i.e. a narrower passband and greater stopband attenuation, if the poles are ‘allocated’ to the weighting function, rather than the polynomial model in the regression formulation. Note that in Fig. 1, the stopband attenuation increases with incremented 𝜅, for a constant VRF; however, Fig. 2 reveals that this is accompanied by an increase in the bandwidth, due to a contraction of the pole radius 𝑝. A similar relationship is apparent in Fig. 6. F. Comparison with Other Filters Generalized IIR-SG filters (𝜅 ≥ 0) proposed in this paper, are compared with ‘regular’ IIR-SG filters (𝜅 = 0) in Figs. 1-6 and FIR-SG filters (with the same effective memory length) in Figs. 1, 2 & 6. In Fig. 6, the SG filters are also compared with FIR filters designed in the frequency domain using the weighted integral squared error (WISE) [28]; and a “maxflat
fractional delay” (MFFD) IIR filter [29]. The WISE filter minimizes the weighted error between the desired response and the actual response (using a passband edge frequency of 𝑓𝑐 = 0.14 cycles/sample and no transition band). The MFFD technique resulted in a filter with
𝐻(𝑧) = [−0.1 + 0.5𝑧−1 ]⁄[1 − 𝑧−1 + 0.5𝑧−2 − 0.1𝑧−3 ]. All filters were tuned to yield a
similar passband width and the same passband group-delay. The two non-SG design methods were chosen because they are also able to produce low-pass filters with very flat pass bands; furthermore, like the SG filters, iterative optimization algorithms are not required to derive the filter coefficients; The shortcomings of the SG filters, and the IIR-SG filters in particular, are evident in the counter-example presented in Fig. 6. The emphasis of this paper is on (very) low-pass filters, and in this role, the IIR-SG filters excel. However in Fig. 6 it is assumed that a wider passband is required, which can only be satisfied using an SG filter with a relatively high-order polynomial model and a moderate effective length (as indicated by a VRF that is not too close to zero) to prevent the contraction of the passband. The resulting numeric complexity of the regular IIR-SG filter is approximately equal to the equivalent FIR-SG filter (i.e. with the same effective memory length) in this case and the IIR complexity only increases further if more poles are added using 𝜅 > 0 to increase the passband width and stopband attenuation. Indeed, other filters such as an FIR-WISE filter or an IIR-MFFD filter, are arguably more appropriate in this case. The passband edge of the FIR-WISE filter may be arbitrarily set at any frequency; furthermore, if it is assumed that passband performance is a priority, then some stopband attenuation can be traded for increased low-frequency magnitude flatness and phase linearity by increasing the relative passband weight. Unfortunately, the parameters of an SG filter do not naturally follow from a cutoff frequency specification; which is why a plot of 3 dB widths for various 𝐵 and 𝑀 combinations is presented in [7] for the FIR-SG filters. The use of multiple poles at a common point on the positive real axis, restricts the ‘influence’ of generalized IIR-SG filters to the low-frequency region. From this crude starting point, the zeros are then automatically placed throughout the complex 𝑧-plane to ‘sculpt’ the response around the unit circle and to realize the ‘intent’ of the time-domain linearregression framework. In contrast, the IIR-MFFD filters are able to have much wider passbands at a lower computational complexity because their poles may be positioned anywhere in the complex 𝑧-plane. However, this additional flexibility is also a potential weakness because it is not known from the outset whether the maximally flat constraint may be satisfied at 𝜔 = 0 for a given combination of numerator and denominator polynomial orders and a specified group delay; and if a solution is feasible, it is not known in advance whether a causal filter will be stable; furthermore, if a stable causal filter is produced, it may have undesirable properties far from 𝜔 = 0. Therefore, many iterations were required to produce an acceptable filter, based mainly on trial and error; however the final result is clearly satisfactory. This highlights one of the most appealing features of SG filters (which are also maximally flat at 𝜔 = 0): the design process is guided by a simple conceptual framework, which has meaning in both the time and frequency domains.
Fig. 1. Comparison of FIR- and generalized IIR-SG filters with 𝐵 = 3 (i.e. cubic model). Impulse response (upper subplot) and frequency response (lower subplots). FIR-SG filter with 𝑀 = 16 and group delay set for linear phase (𝑞 = 7.5). Generalized IIR-SG filters with the same effective memory length and same nominal group delay also shown for different 𝜅 values. Note the monotonic magnitude decay in the stop band for the generalized IIR-SG filters and the increased high-frequency attenuation as the pole multiplicity is increased (via 𝜅).
Fig. 2. Low-frequency magnitude response (upper subplot) and phase response error, showing the deviation from linearity (lower subplot) for the filters in Fig. 1. Generalized IIR-SG filters parameterized by setting the pole location 𝑝, to yield the same VRF as the FIR-SG filter (see legend). The greater highfrequency attenuation of the regular IIR-SG filter (𝜅 = 0) relative to the FIR-SG filter, apparent in Fig. 1 is clearly achieved at the expense of reduced bandwidth although the bandwidth is restored as 𝜅 is incremented.
Fig. 3. Shape of error-weighting function 𝑤(𝑚), for various 𝜅 values (upper subplot) and the VRF as a function of the 𝑞 parameter for various (𝐵, 𝜅) combinations (lower subplot) for the generalized IIR-SG filters with constant 𝑝 = exp(− 1⁄4); locations of VRF minima (closest to zero) are also marked. For a given 𝐵, increasing 𝜅 allows the VRF to be reduced by increasing 𝑞; for a given 𝜅, increasing 𝐵 increases the VRF minima.
Fig. 4. Magnitude response (upper subplot) and phase response (lower subplot) of the filters shown in Fig. 3, for various (𝐵, 𝜅) combinations, with 𝑞 set to minimize the VRF, for a Nyquist null; corresponding lines of constant group delay, for an ideal linear-phase filter, also shown. Note that increasing 𝜅 increases attenuation in the mid-to-high frequencies (i.e. in the stopband); however this is achieved at the expense of increased group delay at low frequencies (i.e. in the passband).
Fig. 5.
Low-frequency magnitude
response (upper subplot) and phase response error (lower subplot) for the filters shown in Fig. 4, for various (𝐵, 𝜅) combinations and optimal 𝑞. Note that improved stopband attenuation brought about by increasing 𝜅 (thus 𝑞) also decreases the width of the low-frequency passband, where magnitude is approximately flat and phase is approximately linear; although this trend is offset by increasing 𝐵.
Fig. 6. Comparison of various non-linear phase low-pass filters with a reduced-delay (𝑞 = 4). SG filters (FIR and IIR) with a model of high degree (𝐵 = 8). Poles of generalized IIR-SG filters with two different 𝜅 values are set for the same effective memory length as the FIR-SG filter. Responses of comparable filters also shown: IIR-MFFD filter and FIR filters designed by minimizing the WISE for two different passband weights (relative to a unity stopband weight). Passband magnitude response (upper subplot), full-band magnitude response (upper subplot inset) and passband phase non-linearity (lower subplot).
V. REALIZATION For low-order polynomial models, the filter response flexibility brought about by the use of generalized IIR-SG filters, as discussed in the previous section (see Figs 1 & 2), is achieved with a lower numeric complexity – the LDE of the (recursive) generalized IIR-SG filters with 𝜅 = 0, 1 & 2 have 6, 7 and 9 non-zero and non-unity coefficients; whereas the (nonrecursive) FIR-SG filter has 12. The complexity differential further increases as the effective length of the FIR- and IIR-SG filters increase (using 𝑀 ≫ 0 and 𝑝 → 1, respectively), for increased frequency selectivity (i.e. lower VRF). As discussed in [12], recursive FIR-SG filter realization are possible; however, when floating-point arithmetic is used, rounding errors due to the cancellation of marginally stable poles with zeros leads to errors when finite-precision arithmetic is used. Due to the maintenance of running sums in recursive FIR filters, numeric errors, once introduced, are never completely removed, unless the filter is reset and reinitialized. Errant behavior is most noticeable: over long time
periods; and/or when processing signals with a large dynamic range or a constant bias; in machines with low numeric precision. This problem is illustrated in Fig. 7. A recursive FIR-SG filter with 𝐵 = 1, 𝑀 = 9, and 𝑞 = 2 was implemented using double-precision arithmetic and the circuit depicted in Fig. 3 of [12] (Note: The author believes that there are some sign errors in [12] – if (5) is correct, then the sign of 𝑐1 in (8) and in Fig. 3 should be reversed). The root-mean-squared error (RMSE) of the estimator was computed over 5000 contiguous segments, where each segment is 105 samples long. The filter was used to process a unity input signal. As time elapses, small differences between the two running sums, which grow ever larger with time, become increasingly difficult to resolve using the available machine precision. The problem is however, largely avoided if a zero-mean input is supplied, as this provides opportunities for the sums to shrink as well as grow. At the end of the final segment for the run shown in Fig. 7, the RMSE of the non-recursive FIR-SG filter and the recursive IIR-SG filter was zero and 1.8 × 10−15 , respectively. A possible realization of a generalized IIR-SG filter, with a pole multiplicity of 5, is given in Fig. 8. This form is appropriate when access to the spectrum is required, as demonstrated in the next section.
Fig. 7. RMSE of the recursive FIR-SG filter as a function of the segment index, showing divergence over time due to a finite numeric precision for a unity input signal.
𝐻𝑥̂ (𝑧) 𝐻ana (𝑧) 𝑦(𝑛)
̂ (𝑛) 𝜷
𝓐
̂ (𝑛) 𝜶
𝓐T
𝑥̂𝑞 (𝑛)
p
𝒒T
p
𝓒
p
p
𝓕(𝑧)
𝐻syn (𝑧)
𝐻pol (𝑧)
(a)
𝑦(𝑛)
𝐻𝜎̂ (𝑧)
𝜎̂𝑥2̂ (𝑛) 𝐻zer (𝑧)
𝓕(𝑧)
𝐻𝑝 (𝑧)
𝐻𝑝 (𝑧)
𝐻𝑝 (𝑧)
+
𝐻𝑝 (𝑧)
𝐻𝑝 (𝑧)
(b)
(c)
̂ (𝑛) 𝜷
𝐻𝜎̂ (𝑧)
(d)
⋅
−1
p
T 𝑦(𝑛)
1 𝑧
𝑝
𝐻𝑝 (𝑧)
+
⋅
𝑐𝑝 VRF
𝜎̂𝑥2̂ (𝑛)
𝑊(𝑧)
Fig. 8. Conceptual circuit schematic for a generalized IIR-SG filter. Thin and thick arrows represent scalar- and vector-valued connections, respectively. Had the spectrum and variance not been required, a much simpler canonical form could instead have been used, with analysis and synthesis operations ̂ (𝑧)⁄𝑌(𝑧) = 𝓐𝓒𝓕(𝑧) and 𝐻syn (𝑧) = combined. (a) The generalized IIR-SG filter is factored into analysis and synthesis sections – where 𝐻ana (𝑧) = 𝔅 ̂ (𝑧) = 𝒒T 𝓐T, respectively – to expose the Laguerre spectrum, which is used here for variance estimation by 𝐻𝜎̂ (𝑧); the component spectrum is 𝑋̂𝑞 (𝑧)⁄𝔅 also optionally provided as an output. (b) & (c) The pole section is a cascading bank of ‘leaky’ integrators. (d) Variance estimation block consisting of parallel vector dot-product and scalar-squaring paths; implements (19b).
VI. ADAPTIVE CONFIGURATION Given the flexibility of the IIR-SG filters, and the ease with which the response may be modified, it is natural to ask whether it is possible to vary the filter parameters online, in response to the characteristics of the input 𝑦(𝑛). Some possible IIR approaches are discussed in [10] and [17], the FIR case is considered in [8]; however, an alternative approach is considered here. The rationale behind the proposed approach is based on the observation that: On the one hand, in closedloop real-time systems, it is generally preferable to have a low filter delay for robust stability margins and system responsiveness; On the other hand, this must be balanced against a desire to strongly attenuate high-frequency noise and interference, which is best achieved with a high delay. In an attempt to satisfy both requirements, an adaptive configuration that employs two complementary generalized IIRSG filters in parallel is proposed. The filters share a common pole section which is cascaded with two customized zero sections (see Fig. 9). After post multiplication by an 𝓐𝓒 matrix of appropriate dimensions to complete the analysis process, the pole section 𝐻pol (𝑧), recursively generates the Laguerre spectrum which is used here for two purposes: firstly, to form the signal estimate; and secondly, to compute the variance of the estimate. Using a simple switching mechanism, a more ‘conservative’ filter with a lower VRF is used when the variance exceeds a specified threshold 𝜆𝜎̂ . A 5-pole IIR filter was considered to yield a reasonable balance between performance and complexity. A pole location of
𝑝 = exp(− 1⁄4) was also deemed to provide a satisfactory impulse response duration and a reasonable starting point for the frequency response. For each filter, the next design issue was to decide whether the poles should be assigned to the polynomial model using 𝐵, or to the error-weighting function using 𝜅. The default filter 𝐻lo (𝑧), assumes a benign noise/interference environment; it uses 𝐵lo = 4 and 𝜅lo = 0 for a relatively wide passband and moderate attenuation in the stopband. Minimization of the VRF yields a low nominal group delay of 𝑞lo = 2.5242 samples. The other, more ‘conservative’, filter 𝐻hi (𝑧), assumes a malign environment; it employs a lesssophisticated polynomial model (𝐵hi = 2) and a more balanced error-weighting function (𝜅hi = 2). Reallocating the poles in this way allows a lower VRF to be attained. Minimization of the VRF yields a high nominal group delay of 𝑞hi = 10.2528 samples. See the Appendix for filter details, Fig. 9 for a schematic of the adaptive configuration, Fig. 10 for the filter responses, and Figs. 11 & 12 for idealized simulation results. 𝑥̂lo (𝑛) 𝑦(𝑛)
lo (𝑧) 𝐻zer
2 (𝑛) 𝜎̂lo
𝐻pol (𝑧)
0 > 𝜆𝜎̂
hi (𝑧) 𝐻zer
𝑥̂hi (𝑛)
𝑥̂𝑞 (𝑛)
1
Fig. 9. Conceptual circuit schematic for the adaptive filter configuration. Two generalized IIR-SG filters with complementary characteristics – 𝐻lo (𝑧) & lo (𝑧) 𝐻hi (𝑧) – with switching logic driven by a variance estimate. Both filters share a common pole section 𝐻pol (𝑧) with customized zero sections, 𝐻zer and hi (𝑧) 𝐻zer in parallel, which are complementary instantiations of the 𝐻zer (𝑧) block in Fig. 8.
Fig. 10. Frequency response of 𝐻lo (𝑧) and 𝐻hi (𝑧), with a low and high nominal group delay in the passband, respectively (𝑞lo = 2.52 & 𝑞hi = 10.25). Both filters have the same number of poles (at 𝑧 = 0.7788) and a zero positioned (at 𝑧 = −1) for a Nyquist null. The former filter (with 𝐵lo = 4, 𝜅lo = 0) is tuned for a wide passband and reasonable stopband attenuation; whereas the latter filter (with 𝐵hi = 2, 𝜅hi = 2) is tuned for enhanced stopband attenuation.
Fig. 11. Illustrative example showing the operation of the adaptive filter in a simulation. Upper subplot: 𝑥(𝑛) is the desired (true) signal – a sinusoid with a period of 80 samples; 𝑦(𝑛) is the input waveform, consisting of signal plus interference – a pulsed envelope modulated by a sinusoid with a period of 4 samples. Lower subplot: using the simple switching mechanism (see Fig. 9), a longer group-delay is only applied when increased high-frequency attenuation is required.
Fig. 12. Operation of the switching mechanism. The filter selection switch (see Fig. 9) is in the ON state, enabling the application of extra high-frequency attenuation (see Fig. 9), when the estimated out-of-band noise power 𝜎̂𝑥2̂ (𝑛) produced by 𝐻lo (𝑧) exceeds a specified threshold. Note the elevated variance output by 𝐻hi (𝑧) prior to the introduction of the interference pulse, due to its narrower passband (see Fig. 10); this filter has a lower expected random error (when processing white noise) but a larger bias error (when processing sinusoids or polynomial models with a degree greater than 𝐵).
VII. DISCUSSION/CONCLUSION A. Why use SG Filters? Use of a polynomial signal model to design digital filters is ideal for data processing in target-tracking systems involving radio-frequency (RF, e.g. radar), electro-optic (EO, e.g. video), infrared (IR) or acoustic (e.g. sonar) sensors; particularly when: the signal-to-noise ratio (SNR) is high, so that measurement-to-track assignment (i.e. data association) is unambiguous, and when the sampling rate is constant, so that filter coefficients do not need to be computed online. When these conditions cannot be satisfied, Kalman filter variants are likely to be more appropriate, particularly when prior models of measurement and kinematic processes are available. In such systems, due to Newton’s laws of motion, sequences of target measurements are likely to be well approximated by a polynomial over long time scales, in azimuth/elevation coordinates (and range for active sensors). A polynomial basis set is also suitable for signal processing in systems when it is necessary to have fine control (and/or make precise measurements) of magnitude and phase in the very-low-frequency region of the spectrum, while strongly attenuating (and/or ignoring) other out-of-band frequencies. SG filters guarantee unity gain and phase linearity at dc and the deviation from ideality, with increasing frequency, decreases as the degree of the polynomial model increases (i.e. they are maximally flat). The SG filters have a ripple-free passband although there is a magnitude ‘hump’ at the band edge which increases as 𝑞 is decreased from the optimal value. Other design techniques, such as those that employ a sinusoidal basis set, are more appropriate when other regions of the spectrum are of interest. B. Why use IIR-SG Filters? IIR-SG filters should be considered when very-low-pass filters with a long impulse response and a short group delay are appropriate. A longer impulse response – as determined by the filter length 𝑀, in the FIR case; or the pole radius 𝑝, in the IIR case – increases the frequency selectivity of a filter. The order of a non-recursive FIR-SG filter is equal to 𝑀 − 1; whereas the order of a IIR-SG filter is equal to 𝐵 + 1. The order of recursive FIR-SG filters is also independent of 𝑀; however, their use is not recommended in machines with finite numeric precision. As the number of multiply and accumulate operations (i.e. the numeric complexity) is approximately proportional to the filter order, IIR-SG filters become more attractive as 𝐵 + 𝜅 decreases, and as the effective length of the filter increases. This observation motivated the emphasis on relatively loworder IIR filters in this paper. Long impulse responses however, reveal the true nature of SG filters in the frequency domain – it is not possible to represent sinusoidal signal components with a short period, using a polynomial model of low degree, over a long timeframe. Thus the effects of increasing 𝑀 (in the finite-memory case) or 𝑝 (in the fading-memory case) should be offset by a commensurate increase in 𝐵 to ensure that the bandwidth (and VRF) of the filter does not vanish. Therefore the advantages of IIR-SG filters discussed above are most likely to be realized when processing low-frequency signals with a narrow bandwidth (or when sampling at a high rate) where SG filters with a low VRF are required. A short group delay (as determined by the 𝑞 parameter) is required for robust stability margins in closed-loop systems and to reduce the processing latency in open-loop systems, particularly in systems with long sampling periods. Using an FIR-SG filter with 𝑞 = 12[𝑀 − 1] is the only way to satisfy a requirement for perfect phase linearity. For short filters, this is certainly not a prohibitive constraint; however as the filter length grows, a point is eventually reached where some loss of phase linearity may be considered to be an acceptable price to pay for reasonable group delay in real-time systems. A short
group-delay decreases magnitude flatness and phase linearity in the passband, decreases stopband attenuation, and increases the overall VRF of a given SG filter – FIR and IIR variants alike. The relative suitability of FIR- and IIR-SG filters is therefore determined by the extent to which their responses are degraded in this low-delay regime. The absence of ripples in the stopband of IIR-SG filters makes them the more attractive alternative when a uniform/consistent response to highfrequency inputs is required. C. Why use Generalized IIR-SG Filters? The ‘regular’ IIR-SG filters with 𝜅 = 0 are ideal for low-latency filters because the exponential weighting function has a maximum at the most recent sample at 𝑚 = 0 and the majority of its weight concentrated near zero. IIR-SG (or FMP) filters have historically been used in closed-loop data-processing (e.g. target tracking) systems, where low, zero or negative delays are desired. This may account for why the generalized form has not so far been reported (to the author’s knowledge). Using 𝜅 > 0 shifts the maximum and mean of the weighting function backwards in time; furthermore as 𝜅 increases, the weight becomes more symmetric around its mean (i.e. decreased skew) and increases less abruptly from 𝑚 = 0. When the production of the estimate is deferred so that polynomial signal is reconstructed near the delayed centroid of a more symmetric weighting function, the frequency selectivity of the filter is greatly improved, i.e. narrower passband width and lower high-frequency gain. Thus significant improvements in low-pass filtering performance are achievable (i.e. the VRF is lowered), when the filter pole multiplicity is increased using 𝜅 rather than 𝐵, provided a small increase in group delay is tolerable. D. Why use Adaptive Generalized IIR-SG Filters? An adaptive filter-bank arrangement was primarily used here as a simple illustrative example to emphasize the three main points raised in this paper: firstly, that generalized IIR-SG filters with complementary responses may be formed by allocating a given number of filter poles to either the degree of the polynomial model for a wider bandwidth or to the shape of the error-weighting function for a lower VRF; secondly that setting the group delay to minimize the VRF results in very good low-pass filters; and thirdly, that factoring the filter into analysis and synthesis stages, allows the Laguerre spectrum to be exposed and used to recursively compute the variance, which permits the quality of the signal estimate to be monitored. Like the Fourier spectrum, the Laguerre spectrum is a very rich source of information which may be used to guide and inform high-level processes, yielding adaptive behavior. Any number of other switching/mixing approaches could have been used; however, the adopted method is sufficient to illustrate the potential of an adaptive approach and the desirable properties of the generalized IIR-SG Filters. APPENDIX: FILTER PARAMETERS Pole stage, common to both filters, 𝐻lo (𝑧) and 𝐻hi (𝑧): 𝑝 = exp(− 1⁄4) = 0.7788 𝒂 = [1 −3.8940 6.0653 −4.7237 1.8394 −0.2865] Zero stage, of low-delay filter, 𝐻lo (𝑧): max 𝐵lo = 4, 𝜅lo = 0, 𝑞lo = 2.5242, VRFlo = 0.1592, 𝐷lo =5
𝒃lo = [0.8142 −1.3017 −0.2775 1.3044 −0.5340] × 10−1 𝓐lo =
4.7032 0 0 0 0 −4.1505 1.1789 0 0 0 × 10−1 3.6628 −2.2284 0.1477 0 0 −3.2324 3.1701 −0.4282 0.0123 0 [ 2.8526 −4.0230 0.8296 −0.0482 0.0008] 1 0 0 0 0 −1 1 0 0 0 𝓒lo = 1 −3 2 0 0 −1 7 −12 6 0 [ 1 −15 50 −60 24] Zero stage, of high-delay filter, 𝐻hi (𝑧): max 𝐵hi = 2, 𝜅hi = 2, 𝑞hi = 10.2528, VRFhi = 0.0592, 𝐷hi =3
𝒃hi = [0 3.5830 −0.7098 −3.3182 0.9746] × 10−3 0.8839 0 0 𝓐hi = [−1.5310 0.1276 0 ] × 10−1 2.1656 −0.3609 0.0113 1 −3 2 0 0 𝓒hi = [−1 7 −12 6 0] 1 −15 50 −60 24
REFERENCES [1] A. Savitzky and M.J.E. Golay, “Smoothing and differentiation of data by simplified least squares procedures,” Anal. Chem. vol. 36, pp. 1627–1639, 1964. [2] P. A. Gorry, “General least-squares smoothing and differentiation by the convolution (Savitzky-Golay) method,” Anal. Chem. vol. 62, no. 6, pp. 570573, 1990. [3] A. D. Thrall, “Some Properties of Least-Squares Digital Filters,” SIAM Journal on Applied Mathematics, vol. 40, no. 2, pp. 169-178, Apr. 1981. [4] P. Steffen, “On digital smoothing filters: A brief review of closed form solutions and two new filter approaches,” Circuits, Syst. Signal Process., vol. 5, no. 2, pp. 187-210, 1986. [5] S. Samadi and A. Nishihara, “The world of flatness,” IEEE Circuits Syst. Mag., vol. 7, no. 3, pp. 38-44, 3rd Qtr 2007. [6] Mohammad Mahdi Jahani Yekta, “A frequency domain proof for the equivalence of the maximally flat FIR fractional delay filter and the Lagrange interpolator”, Digit. Signal Process., vol. 21, no. 1, pp. 13-16, Jan. 2011. [7] R.W. Schafer, “What Is a Savitzky-Golay Filter? [Lecture Notes],” IEEE Signal Process. Mag., vol.28, no.4, pp.111-117, Jul. 2011. [8] S.R. Krishnan and C.S. Seelamantula, “On the Selection of Optimum Savitzky-Golay Filters,” IEEE Trans. Signal Process., vol. 61, no. 2, pp. 380-391, Jan. 2013. [9] N. Morrison, Introduction to Sequential Smoothing and Prediction, McGraw Hill, New York USA, 1969; available online: http://goo.gl/CM9cN. [10] N. Morrison, Tracking Filter Engineering: The Gauss-Newton and Polynomial Filters, Institution of Engineering and Technology, London, 2013. [11] R. Nadjiasngar, M. Inggs, Y. Paichard and N. Morrison, “A new probabilistic data association filter based on composite expanding and fading memory polynomial filters,” in Proc. IEEE Radar Conf., pp. 152-156, May 2011. [12] T. G. Campbell and Y. Neuvo, “Predictive FIR filters with low computational complexity,” IEEE Trans. Circuits Syst., vol. 38, no. 9, pp. 1067-1071, Sep. 1991. [13] D.J. Thornley, “Novel Anisotropic Multidimensional Convolutional Filters for Derivative Estimation and Reconstruction,” in Proc. IEEE Int. Conf. Signal Process. and Commun., pp. 253,256, 24-27 Nov. 2007. [14] Quan Quan and Kai-Yuan Cai, “Time-domain analysis of the Savitzky–Golay filters,” Digital Signal Process., vol. 22, no. 2, pp. 238-245, Mar. 2012. [15] Ç. Candan and H. Inan, “A unified framework for derivation and implementation of Savitzky–Golay filters,” Signal Process., vol. 104, pp. 203-211, Nov. 2014. [16] H. Kennedy, “Recursive Digital Filters With Tunable Lag and Lead Characteristics for Proportional-Differential Control,” IEEE Trans. Control Syst. Technology, vol. 23, no. 6, pp. 2369-2374, Nov. 2015. [17] H. L. Kennedy, “An Adaptive Digital Filter for Noise-Attenuation in Sampled Control Systems,” Int. J. Adapt. Control Signal Process. 2015 (to appear). [18] H. L. Kennedy, “Multidimensional digital smoothing filters for target detection,” Signal Process., vol. 114, pp. 251-264, Sep. 2015. [19] H. L. Kennedy, “Improved IIR Low-Pass Smoothers and Differentiators with Tunable Delay”, in Proc. Int. Conf. Digital Image Computing: Techniques and Applications (DICTA), Adelaide, Australia, 23-25 Nov. 2015 (to appear).
[20] M. Gasca and T. Sauer, “On the history of multivariate polynomial interpolation,” J. Comput. and Applied Math., vol. 122, no. 1–2, pp. 23-35, Oct. 2000. [21] H. L. Kennedy, “Digital-Filter Designs for Recursive Frequency-Analysis,” Journal of Circuits, Systems, and Computers, vol. 25, no. 2, 2016 (to appear). [22] W. Chen, N. Kehtarnavaz and T.W. Spencer, “An efficient recursive algorithm for time-varying Fourier transform,” IEEE Trans. Signal Process., vol. 41, no. 7, pp. 2488-2490, Jul. 1993. [23] S. Tomazic, “On short-time Fourier transform with single-sided exponential window,” Signal Process., vol. 55, no. 2, pp. 141-148, Dec. 1996. [24] M.G. Amin and K.D. Feng, “Short-time Fourier transforms using cascade filter structures,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 42, no. 10, pp. 631-641, Oct. 1995. [25] Weisstein, Eric W. “Z-Transform.” From MathWorld – A Wolfram Web Resource. http://mathworld.wolfram.com/Z-Transform.html [26] M. A. Masnadi-Shirazi and N. Ahmed, “Optimum Laguerre networks for a class of discrete-time systems,” IEEE Trans. Signal Process, vol. 39, no. 9, pp. 2104-2108, Sep. 1991. [27] G. Mandyam and N. Ahmed, “The discrete Laguerre transform: derivation and applications,” IEEE Trans. Signal Process., vol. 44, no. 12, pp. 29252931, Dec. 1996. [28] C. S. Burrus, A. W. Soewito, R. A. Gopinath, “Least squared error FIR filter design with transition bands,” IEEE Trans. Signal Process., vol. 40, no. 6, pp. 1327-1340, Jun. 1992. [29] Xi Zhang, “Maxflat Fractional Delay IIR Filter Design,” IEEE Trans. Signal Process., vol. 57, no. 8, pp. 2950-2956, Aug. 2009.
