Butterworth Filter Prototype

An \(N\)-th order analog Butterworth lowpass filter with arbitrary gain \(H_0\) has transfer function

\[ H_a(s) = \frac{H_0}{\prod_{k=1}^{N}(s - s_k)}, \]

where the poles \(s_k\) are located on a semicircle of radius \(\Omega_c\) in the left half of the \(s\)-plane,

\[ s_k = \Omega_c \, e^{\,j\left(\frac{\pi}{2} + \frac{(2k-1)\pi}{2N}\right)}, \quad k = 1, \ldots, N, \]

ensuring stability.

In this framework, \(H_0 = 1\) is adopted to preserve the spectral shape and establish a direct correspondence with the AR model in the \(z\)-domain.

Matched Z-Transform

The MZT maps each analog pole \(s_k\) to a discrete-time pole \(z_k = e^{s_k T}\), where \(T = 1/f_s\) is the sampling period. Since the Butterworth prototype is an all-pole filter with \(H_0 = 1\), the system has no finite zeros, and the discrete-time transfer function reduces to

\[ H(z) = \frac{1}{\prod_{k=1}^{N}\left(1 - z_k z^{-1}\right)}, \quad z_k = e^{s_k T}, \]

which preserves the all-pole structure required by an AR model. Expanding the denominator polynomial using Viète’s formulas yields

\[ \prod_{k=1}^{N}\left(1 - z_k z^{-1}\right) = 1 + a_1 z^{-1} + a_2 z^{-2} + \cdots + a_N z^{-N}, \]

and the AR coefficients are obtained directly as \(\phi_i = -a_i\), so that the transfer function matches the standard AR form

\[ H(z) = \frac{1}{1 - \phi_1 z^{-1} - \cdots - \phi_N z^{-N}}. \]

Correspondence Between Filter Order and AR Order

The following result, established in Bras-Geraldes, Rocha, and Martins (2026), formalises the equivalence between the number of digital poles and the number of lags in the autoregressive model.

Proposition.

Let \(H(z)\) be the transfer function of an autoregressive model \(y(n)\) of order \(p\). Then \(H(z)\) has precisely \(p = N\) lags, where \(N\) is the number of discrete-time poles \(z_k = e^{s_k T}\), \(k = 1, \ldots, N\). In addition, the AR coefficients satisfy \(\phi_i = -a_i\), where the coefficients \(a_i\) are given by Viète’s formulas, for \(i = 1, \ldots, p\).

This Proposition has a direct practical consequence: the Butterworth filter order \(N\) simultaneously determines the autoregressive model order \(p\). This property provides a direct interpretability link between spectral design parameters and the temporal dynamics of the resulting stochastic process. In the example that follows, the partial autocorrelation function (PACF) is used as a heuristic to fix \(p = 3\), which immediately sets \(N = 3\) for the Butterworth design.

The Butterworth-Induced Autoregressive Model

The principal theoretical result of Bras-Geraldes, Rocha, and Martins (2026) is the following theorem, which formally defines the BTW-AR model and establishes the conditions under which it is well-defined.

Theorem (Butterworth-Induced Autoregressive Model).

Let \(y(n)\) be an autoregressive model of order \(p\), let \(H_a(s)\) be the Butterworth filter used to model its all-pole structure, and let \(\phi_i^{(BW)}\) denote the induced AR coefficients obtained via Viète’s formulas, for \(i = 1, \ldots, N\), where \(N\) is the minimum filter order required to satisfy the stopband specification

\[ N = \left\lceil \frac{1}{2} \frac{\log\!\left(10^{A/10} - 1\right)} {\log\!\left(f_{Nyquist}/f_c\right)} \right\rceil, \]

and \(A\) is the attenuation in decibels at the Nyquist frequency \(f_{Nyquist}\). If a scaling factor \(\alpha\), estimated by minimising the sum of squared residuals, satisfies \(\varepsilon(n) \sim N(0,\sigma^2)\), then the Butterworth-induced autoregressive model is defined by

\[ y(n) = \alpha \sum_{i=1}^{N} \phi_i^{(BW)}\!\left(N,\, A,\, f_{Nyquist}\right)\, y(n-i) + \varepsilon(n). \]

Model Fitting Pipeline

Given a training series \(\{y_t\}\), the fitting procedure is:

1. Stationarity — apply successive differencing (guided by the Augmented Dickey-Fuller test) until stationarity is achieved.
2. Centering — subtract the training mean \(\bar{y}\).
3. Cross-validation — search over a grid of candidate \((A, N)\) values using nested rolling-origin CV to select the combination that minimises out-of-sample RMSE.
4. Scaling — estimate a scalar \(\alpha\) by regressing the observed series onto the one-step-ahead AR predictions using the chosen method (LS, LAD, or Huber).

Setup

library(BTWAR)

set.seed(123)

fs <- 12000   # sampling frequency (Hz)

BTWAR.R

Simulating Data

We simulate an AR(3) process with known coefficients and split it into training (70%) and test (30%) sets using simulate_ar_split().

phi_obs <- c(0.34, -0.22, 0.16)

sim_data <- simulate_ar_split(
  phi_real   = phi_obs,
  n          = 1000,
  burn       = 200,
  prop_train = 0.7
)
#> Simulating AR(3) series ... done.

dados <- list(
  train = as.numeric(sim_data$train),
  test  = as.numeric(sim_data$test)
)

BTWAR.R

The simulated AR coefficients and their corresponding Z-plane poles are:

cat("True AR coefficients:\n")
#> True AR coefficients:
print(phi_obs)
#> [1]  0.34 -0.22  0.16

cat("\nTrue Z-plane poles:\n")
#> 
#> True Z-plane poles:
print(sim_data$poles)
#> [1]  0.5154289+0.0000000i -0.0877144+0.5502066i -0.0877144-0.5502066i

BTWAR.R

Exploratory Analysis

Following Bras-Geraldes, Rocha, and Martins (2026), the autoregressive order \(p\) is first identified using the PACF as a heuristic diagnostic tool. Because the MZT establishes a correspondence between the AR order and the Butterworth filter order, fixing \(p\) also determines the filter order \(N\).

Once \(N\) is fixed, the nested cross-validation procedure searches only over the stopband attenuation \(A\), reducing the hyperparameter search to a one-dimensional grid.

pacf(dados$train, main = "PACF – Training Series")

PACF of the training series. Significant spikes at lags 1, 2, and 3 suggest AR order 3, which fixes the filter order N = 3.

BTWAR.R

An inspection of the PACF suggests that an AR model with \(p = 3\) lags provides an adequate representation of the training series. This fixes \(N = 3\) prior to cross-validation.

Fitting the BTW-AR Model

p <- 3

res <- btwar_fit(
  y_tr_raw = dados$train,
  y_te_raw = dados$test,
  fs       = fs,
  method   = "ls",
  N_vec    = p:p
)
#> Warning in tseries::adf.test(x): p-value smaller than printed p-value

cat("BTW-AR Performance:\n")
#> BTW-AR Performance:
print(res$performance)
#> $rmse_train
#> [1] 1.013619
#> 
#> $rmse_test
#> [1] 0.9561448

BTWAR.R

Time-Series Predictions

The plots below show the one-step-ahead predictions of the BTW-AR model on the training and test sets. The model tracks the observed signal closely, capturing the dominant low-frequency dynamics of the underlying process. .

plot(res, dataset = "train", fs = fs,
     colour_observed = "goldenrod",
     colour_btwar    = "blue",
     lwd_observed    = 1.2,
     lwd_btwar       = 1.2)

Training set: observed series vs. BTW-AR predictions.

BTWAR.R

plot(res, dataset = "test", fs = fs,
     colour_observed = "goldenrod",
     colour_btwar    = "blue",
     lwd_observed    = 1.2,
     lwd_btwar       = 1.2)

Test set: observed series vs. BTW-AR predictions.

BTWAR.R

Spectral Interpretation: Bode Magnitude Diagram

A distinctive feature of the BTW-AR framework is that the fitted model admits a direct frequency-domain interpretation. The Bode magnitude diagram below shows the frequency response of the Butterworth prototype underlying the selected model.

plot_bode(res, fs = fs)

Butterworth magnitude response for the selected BTW-AR model.

BTWAR.R

The diagram reveals two distinct spectral regions:

• Passband \([0, f_c]\) — the model captures this band with a nearly flat gain, meaning that the low-frequency structure of the signal is represented with minimal distortion.
• Transition band and stopband \([f_c, f_{Nyquist}]\) — components in this region are progressively attenuated, reaching \(A\) dB of suppression at the Nyquist frequency.

The cutoff frequency \(f_c\) is therefore a direct, interpretable summary of the model: it marks the spectral boundary between the signal dynamics that the BTW-AR model captures and the higher-frequency variability that it treats as noise. This frequency-domain transparency is not available in conventional data-driven AR models, where the spectral implications of the estimated coefficients are implicit and not immediately readable.

Frequency Amplitude Spectrum

The amplitude spectrum of the training series provides a complementary view to the Bode magnitude diagram. While the Bode diagram characterizes the frequency response of the Butterworth filter prototype, the amplitude spectrum reveals the actual distribution of energy in the observed signal, allowing a direct visual assessment of how well the cutoff frequency \(f_c\) separates the dominant spectral content from the high-frequency components treated as noise by the model.

plot_freq_amplitude(
  dados$train,
  fs = fs,
  fc = res$parameters$fc_opt
)

Amplitude spectrum of the training series. The cutoff frequency fc separates the dominant signal energy from the high-frequency components.

BTWAR.R

The vertical dashed line marks the cutoff frequency \(f_c\) selected by cross-validation. Components to the left of this boundary are captured by the BTW-AR model; components to the right are attenuated according to the Butterworth roll-off visible in the Bode diagram above.

Z-Plane Pole Structure

The Z-plane pole plot compares the poles of the fitted BTW-AR model with those of the true data-generating process. Poles inside the unit circle correspond to stable models.

df_true <- data.frame(
  Re = Re(sim_data$poles),
  Im = Im(sim_data$poles)
)

plot_zpoles(
  res,
  external_list = list(
    "True Signal Z-Poles" = df_true
  )
)

Z-plane poles: BTW-AR (selected) vs. true signal.

BTWAR.R

The pole locations of the BTW-AR model are entirely determined by the Butterworth design parameters \((N, A)\) and are not estimated from the data. Despite this structural constraint, the selected model places its poles in a region of the Z-plane consistent with the dominant dynamics of the true process, reflecting the effectiveness of the cross-validation procedure in identifying a spectrally appropriate model.