Bayesian Deconvolution

Besides FFT filtering, the Bayesian or maximum-likelihood deconvolution introduced by Kennett et al. (1978) offers a robust, non-negative route to recover the logarithmic time-constant spectrum \(R(\zeta)\).

Idea in a nutshell

  • Treat the convolution

    (29)\[h(z)=\left(R\ast w_z\right)(z)\]

    as the generation of observations \(h\) by a source distribution \(R\) through a blurring kernel \(w_z\).

  • Apply Bayes’ theorem iteratively: update the current guess of \(R\) by comparing the synthetic response produced through \(w_z\) with the actually measured \(h\).

Discretisation

Bayesian deconvolution requires an evenly spaced logarithmic grid (\(z\)-axis). The continuous integral becomes a sum

(30)\[h[z] \;=\; \sum_{\zeta} w_z[z-\zeta]\;R[\zeta].\]

With square-bracket indices the kernel is a Toeplitz matrix (\(W_{kj}=w_z[k-j]\)):

(31)\[h_k \;=\; W_{kj}\,R_j.\]

Applying Bayes’ theorem

For vectors Bayes’ rule reads

(32)\[P(R_i\!\mid h_k) \;=\; \frac{ P(h_k\!\mid R_i)\,P(R_i) }{ \sum_j P(h_k\!\mid R_j)\,P(R_j) }.\]

Identify

  • likelihood \(P(h_k\!\mid R_j)=W_{kj}\),

  • prior on data \(P(h_k)=h_k\).

Choosing a flat prior on \(R\) (initialised with \(h_k\)) and enforcing probability conservation gives the classic Richardson–Lucy (RL) iteration specialised to NID:

(33)\[R_i^{(n+1)} \;=\; R_i^{(n)} \sum_{k} \frac{h_k\,W_{ki}} {\sum_{j}W_{kj}\,R_j^{(n)}}.\]

Key properties

  • Energy (area) preserving Because \(\int w_z\,\mathrm{d}z=1\), the cumulative thermal resistance after deconvolution equals the original steady-state value.

  • Non-negativity All terms are positive; negative artefacts cannot occur.

  • No explicit regulariser needed – early stopping of the iterations itself controls over-fitting to noise.

Practical guidelines

  • Initial guess Set \(R^{(0)}=h\) or any smooth positive approximation. Convergence is monotonic but slow.

  • Even spacing Resample the measured impulse response to a uniform \(\Delta z\); the LOWESS/SURE filter (see Optimal regression filtering) can supply both smoothing and resampling.

  • Stopping criterion Iterate until the synthetic step response

    (34)\[h^{(n)} = W\,R^{(n)}\]

    matches the measured one within the estimated noise level (e.g. \(\chi^2\) test) or until a fixed iteration count (1000 – 5000) is reached.

  • Zero handling If the denominator \(\sum_j W_{kj} R_j^{(n)}\) vanishes for wide spectra, skip the affected \(k\) or add a tiny \(\varepsilon\) to avoid division by zero.

  • Underflow Late iterations may drive negligible bins below floating- point precision. Mask them out to save computation without affecting the result.