pyDOF: a Python library for the design of discrete forward and inverse filters

Abstract: In this work, we present pyDOF, a Python-based software library which provides a domain-specific framework for the design of symmetric, physical-space, forward as well as inverse discrete filters. pyDOF is based on a constrained optimisation framework developed in our previous work [1, 2]. This framework allows the user to impose a wide range of constraints on the discrete filter transfer-function such as monotonicity, positivity, value-fixing, gradient-smoothing etc. amongst many others. pyDOF additionally includes an adaptive filter stencil selection option, and a van Cittert-based inverse-filter design with a user-controlled reconstruction order. The filter coefficients are computed automatically, and saved to a plain text file which can be readily parsed by any programming language. pyDOF can be used to design a wide range of low-pass, high-pass, multi band-pass/band-stop etc. discrete filters. In addition, due to its generality and abstraction, pyDOF can be used to design specific filters for user-defined target filter transfer functions. Although developed primarily for application to computational fluid dynamics simulations, pyDOF can be used to design discrete filters for a wide range of signal processing applications.

• The paper introduces pyDOF, an open-source Python library that uses constrained optimisation to design symmetric physical-space forward and inverse filters with user-controlled van Cittert reconstruction order.
• Benchmark results show pyDOF achieves competitive mean-squared error while strictly limiting pass-band overshoot to 0.05%, enforcing unity response at zero wavenumber, and attenuating the grid cut-off.
• The library’s adaptive designs recover resolved signal content in Gaussian-filter, turbulent-flame, and image tests, but higher inverse orders require larger stencils and greater computational cost.

pyDOF is an open-source Python library for designing symmetric, physical-space discrete filters, covering both forward filtering and inverse (deconvolution) filter design within a single constrained-optimisation framework. The library builds on a framework previously developed by the authors for large eddy simulation (LES) applications (2606.26830), and its distinguishing feature relative to standard signal-processing tools is support for inverse filter design with explicit user control over the reconstruction order.

Motivation and positioning

The primary application domain is computational fluid dynamics, where discrete filters underpin turbulence and combustion models such as dynamic Smagorinsky, scale-similarity, and mixed models, as well as deconvolution-based approaches that reconstruct unfiltered fields from filtered ones. In LES, filtering operations are applied to 3D fields on meshes with millions of points, so filters must be accurate yet compact. Physical-space convolution is preferred over Fourier-based approaches for ease of implementation, parallelisation, and boundary-condition handling.

The authors argue that conventional coefficient derivation via truncated Taylor-series expansions produces filters that perform poorly at high wavenumbers near the grid cut-off $kh=\pi$, are not consistent (constant signals may be damped or amplified), and may have negative transfer functions—rendering them unsuitable for deconvolution modelling since reconstruction becomes unstable when $\hat{G}_d(kh)\le 0$ for some resolved wavenumber. Classic signal-processing methods address forward filter design well: Kaiser windows offer efficient design with controllable transition width; Parks-McClellan yields equi-ripple optimal filters; least-squares methods minimise global error but without ripple or transition-width control. However, none of these provide independent pass-band/stop-band control, consistency constraints, positivity guarantees, or—in the case of most published work—inverse filter design within a unified framework. pyDOF fills this gap using constrained optimisation with van Cittert regularisation for the inverse problem.

Mathematical formulation

For a 1D mesh with spacing $h$, the forward filter is a symmetric stencil convolution $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$ with real transfer function $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$. The inverse filter operates analogously on the filtered field with coefficients $b_l$ and half-stencil $M_I$.

The inverse target derives from van Cittert iteration. For $b=1$, after $N$ iterations the reconstructed transfer function is

$$\hat{Q}_d^N(k_rh)=1-\left(1-\hat{G}_d(k_rh)\right)^{N+1}.$$

Two properties follow directly: any wavenumber with $\hat{G}_d(kh)\le 0$0 can never be recovered, whereas if $\hat{G}_d(kh)\le 0$1 then $\hat{G}_d(kh)\le 0$2 as $\hat{G}_d(kh)\le 0$3. This makes positivity of the forward transfer function a prerequisite for stable deconvolution, which pyDOF enforces via constraints. The user selects $\hat{G}_d(kh)\le 0$4 explicitly, thereby controlling how much of the resolved wavenumber range is reconstructed—a parameter most iterative schemes leave implicit through runtime cost.

Optimisation minimises either mean-squared error ($\hat{G}_d(kh)\le 0$5) or mean absolute error ($\hat{G}_d(kh)\le 0$6) between $\hat{G}_d(kh)\le 0$7 and a target $\hat{G}_d(kh)\le 0$8, optionally augmented by a gradient-smoothing penalty weighted by $\hat{G}_d(kh)\le 0$9. For inverse filters, the objective matches $h$0 against $h$1 (or its continuous-domain analogue $h$2). Built-in targets include Gaussian, Helmholtz, implicit (Padé-type), Butterworth, and arbitrary sharp multi-band ideal filters defined by band edges; users may supply custom targets through class inheritance.

Constraints are specified by keywords: value fixing at zero or one (fixOne, fixZero), lower/upper bounds over wavenumber intervals (luBounds), positivity, bounding within $h$3, monotonicity of the transfer function gradient, and an adaptive switch instructing the optimiser to grow the stencil size until a target error metric is met. The value-fixing constraints deserve emphasis: they guarantee $h$4 (consistency for constant signals) and $h$5 (complete attenuation at grid cut-off), properties that Kaiser-window, Parks-McClellan, and unconstrained least-squares designs cannot enforce simultaneously.

Implementation

The library exposes a small object-oriented API built around five core classes. Wavenumber defines the sampled range $h$6; Filter subclasses encapsulate continuous target transfer functions; DiscreteFilter holds the coefficients and stencil configuration for forward or inverse filters; ReconstructedTransferFunction computes $h$7 or $h$8 from a discrete filter (optionally combined with a physical-space Filter); and Optimiser assembles and solves the constrained problem, updating coefficients in place. Coefficients are exported as plain text parseable from any language. The design is deliberately decoupled from any specific optimiser, so alternative solvers can be substituted while retaining the constraint and target abstractions.

Validation

Low-pass filters. For an ideal low-pass target with cut-off $h$9 and stop-band attenuation of 60 dB, stencil sizes of 73 and 145 were derived from the standard formula $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$0 for transition widths $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$1 and $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$2. pyDOF was benchmarked against Kaiser-window, Parks-McClellan, and least-squares designs, with a uniform ripple bound $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$3 imposed over both bands plus the two value-fixing constraints. Least-squares achieves the lowest mean-squared error but exhibits severe Gibbs overshoot (~9% in the pass-band). pyDOF's overshoot is exactly bounded at 0.05%—a direct consequence of the imposed constraint—and its MSQE is competitive:

Method	MSQE ($\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$4=73)	Max pass-band overshoot
Kaiser window	4.91E-03	1.20E-01
Parks-McClellan	4.77E-03	6.32E-02
Least-squares	2.80E-03	9.00E+00
pyDOF	4.80E-03	5.00E-02

Trends are similar at $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$5 and for multi-pass-band filters with transition bands centred at $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$6 and $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$7. A physical-space test on a 512-point periodic signal containing pass-band components and stop-band components including the grid cut-off confirms correct attenuation behaviour. Computationally, Kaiser windows remain cheapest, followed by Parks-McClellan, then the optimisation approaches—an expected trade-off given that filter design is a one-time offline cost.

Forward and inverse Gaussian filters. With $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$8 and adaptive stencil selection (forward error tolerance $\bar{u}_i=\sum_{l=-M_F}^{M_F} g_l u_{i+l}$9, inverse $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$0), the adaptive procedure selected $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$1 for the forward filter, while the inverse stencil grows with reconstruction order: $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$2, 14, and 26 for $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$3, 500, and 5000 respectively. Forward MSQE is $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$4; inverse MSQE remains below $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$5 across all orders, with maximum overshoot in $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$6 rising from 0.12% to 0.75% as $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$7 increases. A physical-space test shows recovery of the original signal down to the largest resolved wavenumber component. Note the monotone trend: higher reconstruction order demands larger inverse stencils and slightly larger composite error, so the reconstruction order trades accuracy against convolution cost.

2D turbulent flame and images. Applying dimensional splitting to 1D filters designed for an extreme implicit target ($$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$8, near the boundary of the admissible range $$\hat{G}_d(k_rh)=g_0+2\sum_{l=1}^{M_F}g_l\cos(k_rh\,l)$$9) with reconstruction order limited to $b_l$0 for cost reasons, the adaptive optimiser produced $b_l$1 and $b_l$2. Filtering a DNS density field of a turbulent premixed V-flame removes substantial fine-scale structure; reconstruction at this modest order recovers much of it. The same filter pair applied to 2D images recovers the originals almost completely at $b_l$3, illustrating applicability beyond CFD.

Limitations and open questions

Several caveats bear directly on the results. First, the framework is restricted to symmetric, explicit, real-coefficient filters obtained via dimensional splitting; asymmetric or truly multidimensional stencils are outside its scope. Second, the validation against classical methods covers ideal low-pass and multi-pass targets only; no comparison is reported for smooth targets such as Gaussian or Butterworth against alternative inverse-design approaches, so the claimed advantage rests largely on the unique capability set rather than uniformly superior error metrics—indeed, unconstrained least-squares attains lower MSQE in every benchmarked case. Third, the inverse design inherits van Cittert's structural limitation: wavenumbers where the forward transfer function vanishes are unrecoverable regardless of $b_l$4, making the positivity constraint essential rather than optional. Fourth, computational cost of the optimisation exceeds classical methods, and the paper does not report optimiser runtimes quantitatively. Finally, the 2D demonstrations use a fixed, modest reconstruction order chosen for cost; the interaction between stencil growth, reconstruction order, and 3D convolution expense on production-scale LES meshes is demonstrated only by reference to prior work rather than established here.

Conclusion

pyDOF provides a unified, constraint-driven framework for designing symmetric discrete forward and inverse filters, with adaptive stencil selection, explicit reconstruction-order control, and extensible target transfer functions. Benchmark results show error metrics comparable to Kaiser-window and Parks-McClellan designs while enforcing consistency and complete grid cut-off attenuation that those methods cannot impose, and the inverse-filter capability—validated on turbulent-flame and image data—is the library's principal differentiator. The main open questions concern extension beyond symmetric 1D-split filters and quantitative cost characterisation relative to classical design methods.