11institutetext: Computer Vision Laboratory, Dept. EE, Linköping University,
Sweden 22institutetext: Saab, Sweden
22email: {emil.brissman, per-erik.forssen<EMAIL_ADDRESS>
# Camera Calibration without Camera Access -
A Robust Validation Technique for Extended PnP Methods

Emil Brissman 1122 0000-0002-0418-9694 Per-Erik Forssén 11 0000-0002-5698-5983
Johan Edstedt 11 0000-0002-1019-8634

###### Abstract

A challenge in image based metrology and forensics is intrinsic camera
calibration when the used camera is unavailable. The unavailability raises two
questions. The first question is how to find the _projection model_ that
describes the camera, and the second is to detect incorrect models. In this
work, we use off-the-shelf extended PnP-methods to find the model from 2D-3D
correspondences, and propose a method for model validation. The most common
strategy for evaluating a projection model is comparing different models’
residual variances—however, this naive strategy cannot distinguish whether the
projection model is potentially underfitted or overfitted. To this end, we
model the residual errors for each correspondence, individually scale all
residuals using a predicted variance and test if the new residuals are drawn
from a standard normal distribution. We demonstrate the effectiveness of our
proposed validation in experiments on synthetic data, simulating 2D detection
and Lidar measurements. Additionally, we provide experiments using data from
an actual scene and compare non-camera access and camera access calibrations.
Last, we use our method to validate annotations in MegaDepth.

Figure 1: Left: Standardised residuals for a correct model with one distortion
parameter using our robust scale estimate. I.e. the residuals are not affected
by the model error. Middle: Standardised residuals from images under more
distortion, for an incorrect model using a non-robust scale estimate. Right:
Standardised residuals for an incorrect model, using a robust scale estimate.

##  1 Introduction

Intrinsic camera calibration is a fundamental computer vision problem. It
involves finding the parameters that allow the conversion of pixel coordinates
to bearing angles [12]. It is possible to use the camera for metrology using a
calibration. In the single-view case, metrology means measuring the lengths
and angles of objects depicted in an image. As an extension, it is the
underpinning of single view 3D reconstruction [4]. Metrology has many
applications, including non-contact measurements, sensor fusion, and forensic
analysis.

Traditionally, intrinsic calibration is a semi-automatic process, which
involves imaging of calibration objects [30, 28]. Such calibration allows
controlled accuracy; however, access to the camera is required. In forensic
analysis, the camera is only sometimes available, depending on the received
material. Therefore, we aim to facilitate measurements in an image when the
camera is unavailable. Using a calibration profile from a camera of the same
model often works well, but the accuracy is unknown in this approach and
should thus be avoided in forensics.

In the Perspective-n-Point (PnP) problem, the goal is to estimate the camera
pose given a set of 2D-3D point correspondences. Early methods assume a
calibrated camera, and only estimate translation and rotation parameters [7].
More recent variants of PnP also estimate the intrinsic camera parameters [14,
21]. These _extended PnP methods_ (xPnP) do not require the camera to be
available, in contrast to calibration pattern methods [30, 28]. However, they
introduce new challenges such as 2D-3D matching and validation.

In this work, we attend to the validation of camera calibration for forensic
metrology applications [2]. Usually, a model is assumed to be validated if it,
on average, has low residuals. However, this approach will not provide any
measure of uncertainty in the image plane. Moreover, deriving the uncertainty
is challenging because the amount of distortion scales non-linearly with the
distance to the camera centre. Thus, we treat noise modelling as a robust
regression problem and predict a residual scaling for each 2D-3D
correspondence. When the model is correct, we assume the scaled residuals to
follow a standard normal distribution (Figure 1). Next, to verify this
assumption, we use a hypothesis test. Simulated data, an indoor scene and
MegaDepth [18], with annotated cameras depicting different scenes, demonstrate
our proposed validation.

Contributions Our contributions are as follows: (i) We propose a method for
testing residuals based on variance predictions and standardisation. (ii) We
suggest using xPnP methods for unavailable cameras as input to our method,
given 2D-3D correspondences. (iii) An empirical estimate of the variance
scales residuals poorly. Instead, we propose a predictive noise model to scale
individual residuals over the 2D detector and projected 3D noise. (iv) We
analyse the effectiveness of our method in quantitative and qualitative
experiments and demonstrate its ability to significantly predict incorrect
models, also when the mean of the residuals is low.

###  1.1 Background Motivation

At the Swedish National Forensic Center (NFC), the task is to collect
information linked to crimes without the possibility of misinterpretations
when used in the Swedish court system. At the time of writing, NFC uses the
Zhang method [30] for metrology. However, this commonly accepted practice only
applies to images where the camera is available. Lidar scanning is a standard
technique in forensic investigations, and in many legal cases, the depicted
location is revisited for scanning. Using this working methodology, Olsson
[22] first investigated the validation of xPnP methods that forms the basis of
this work. In [22], model correctness is assessed by checking if the empirical
mean of the re-projected sample distances is within the two centre quartiles.
However, this decision will prefer incorrect models since outliers will expand
the decision range.

###  1.2 Ethical Consideration

This work does not concern any police investigations or legal cases. Instead,
the method we propose analyses residuals using synthetic data, available
benchmark data, and a snapshot of a fictional crime scene provided by a police
agency. These are all free from apparent ethical dilemmas.

##  2 Related Work

Semi-automatic Calibration Camera calibration is a broad subject found in many
areas of industry and research. However, the most common camera calibration
practice is to use a printed pattern on a planar surface. This strategy was
proposed by Zhang [30], who suggested using a checkerboard pattern with
equidistant squares of black and white. The inner corners of the pattern form
unambiguous features that are easy to find. Detecting several of these
features, also called saddle points, between different views allows camera
parameters to be estimated. Each detected saddle point in each picture is
assigned to its corresponding point on the checkerboard. This set of
correspondences represents a series of homographies, determining the intrinsic
and extrinsic parameters for one or more cameras. In the case of Tsai [28],
camera calibration depends only on one view of a co-planar checkerboard
pattern. The Zhang method [30], instead depends on at least three pictures of
a planar checkerboard pattern. More recently, deep methods like Li _et al._
[17] take a single image as input and jointly learn to predict distortion
coefficients and optical flow from images with lens-type annotations. However,
this problem only concerns visual quality and provides no model accuracy
assessment.

Perspective-n-Point The PnP problem [7] refers to finding a rigid
transformation from 2D-3D point matches. That is, to estimate the rotation
matrix ${\bf R}$ and the translation vector ${\bf t}$ describing the camera
pose in the coordinate system of the 3D points, assuming that the intrinsic
parameters never changed during the sampling of a scene. The minimal but
ambiguous case, P3P [9, 19, 23] is not considered in this work.

Lepetit _et al._ [16] (EPnP), reduced the computational complexity to $O(n)$
operations. Although the convergence is fast, the solution depends on
initialization and global convergence is not guaranteed. Later works
recognized the need to include intrinsic parameters to generalize application
tasks [21][14]. Nakano [21] extends the PnP problem by including intrinsic
parameters and dividing the parameter estimation into different stages. Radial
distortion and equal focal length horizontally and vertically are assumed, as
well as fixating the principal point to the image center. Larsson _et al._
[14] instead require a minimal correspondence set and add a local optimization
step [15]. In our work, we propose a method to validate xPnP methods, for
application in forensic analysis, by _Goodness-of-Fit_ (GoF) testing between
distributions. That is, we do not improve the methods [21] and [14].

Empirical Performance Evaluation

Works by Wang _et al._ [29] and Thai _et al._ [27] are related to our work.
Wang _et al._ [29] propose a method to test hypotheses about the effect of
conceptual changes in deep classification models. That is, if the difference
is the probable reason behind the increase in accuracy. Thai _et al._ [27]
propose to identify cameras by raw pixel intensities. Similar to our approach,
two quantities parametrize the intensity variance—analog gain, controlled by
the camera ISO setting, and electronic noise caused during sensor readout. For
an unknown camera the parameters are first estimated and secondly tested
against known parameter values (null hypothesis).

Goodness of Fit The goodness of fit testing is one of the fundamental tasks in
statistics. In this work, we focus on normality testing due to normal
distributions being a good model for uncertainty in projective geometry [13,
8]. Still, our approach could easily be generalized to GoF tests for arbitrary
distributions.

There are a large variety of proposed statistics for normality testing, of
which the Kolmogorov-Smirnov (KS) [20], D’Agostino-Pearson (DAP) [6], and
Shapiro-Wilk (SW) [26] tests are well known. These tests all seek to maximize
the power,  _i.e._ , minimizing the risk of the null hypothesis being
accepted, given that an alternative hypothesis is correct. We discuss those
further in Section 3.5 and test all three in Section 4.

Single View Metrology Metrology is the study of measurement. In the context of
computer vision, single view metrology [5] involves estimating, _e.g._ ,
angles and lengths, from a single image. In all metrology, an accurate measure
of uncertainty is crucial, and in particular in the forensic setting. Previous
work in single view metrology has focused on the undistorted (but often
uncalibrated) case [5], with model uncertainty assumed to be normally
distributed [5, 13, 8, 3]. At inference time, these uncertainties can be
propagated by first order propagation or by Monte Carlo simulation. However,
in those works, both the uncertainty estimation and propagation requires _a
priori_ knowledge of the noise levels and estimation method, and implicitly
assume the estimated model is approximately correct. In contrast to those
methods, our approach

  1. 1.

Is estimation agnostic, _i.e._ , we can treat the estimators as black boxes.

  2. 2.

Generalises to arbitrary projection models.

  3. 3.

Does not implicitly assume that the estimated model is approximately correct.

In particular, perturbation theory, as used in previous work, does not provide
a reliable measure of the trustworthiness of the estimated model, it simply
provides an approximate measure of the estimation sensitivity to the input. In
contrast, our method directly measures trustworthiness by testing the
hypothesis of the matches being generated from the estimated model.

##  3 Method

We propose a method that compares observed and expected noise levels. The
method takes residual values as input, given a calibration computed from an
xPnP method and 2D-3D correspondences. We decompose the residual error for
each correspondence as three additive terms: (i) 2D detector noise, (ii) 3D
detector noise projected into the image, and (iii) model noise. The expected
model noise is zero for the correct model, not affecting the residual
distribution in any direction. We describe this in Section 3.3. To handle
unexpected model noise, Section 3.4 details robust regression over (i) and
(ii) to obtain a scale value for each 2D point. Finally, we assume the scaled
residuals are drawn from a standard normal distribution and test this using a
GoF test. We motivate our preferred choice of test in Section 3.5. We consider
all points to influence the validation decision and believe this to improve
applications in forensic analysis. We begin with an example to get a good
intuition of our approach.

###  3.1 Motivating Example

Consider a correct data model $y=x$ and an (incorrect) hypothesis
$h_{\text{bad}}:y=x+0.5x^{5}$. Under the assumption that $y$ is observed with
some Gaussian noise, the residuals $r$ of the true model will be distributed
as $\mathcal{N}(0,\sigma_{y}^{2})$. In contrast, the residuals of the
incorrect hypothesis are typically _significantly_ different from the expected
distribution (as shown in Figure 1). Thus, if $\sigma_{y}$ is known, a simple
hypothesis test is whether $\frac{r}{\sigma_{y}}\sim\mathcal{N}(0,1)$.
However, in real world scenarios $\sigma_{y}$ is typically not known and needs
to be estimated. Since incorrect hypotheses typically contain outliers, it is
important with a robust estimate of the noise level. We show these steps in
Figure 1. It is clear that $h_{\text{bad}}$ produces a tailed residual
distribution that does not follow the expected Gaussian curve. Hence we can
use the KS test [20] to validate the produced models.

Underfittning and Overfitting It is common for a complex model to be optimized
to fit the data $y$ perfectly. We can describe overfitting and underfitting as
a constant multiplication of $\sigma_{y}^{2}$, yielding residuals distributed
as $\mathcal{N}(0,a\sigma_{y}^{2})$. When $a<1$ the model is overfitted, and
when $a>1$ it is incorrect (underfittning). The following sections describe
how we can apply this intuition to validate a camera calibration.

###  3.2 Camera Calibration

Calibration fundamentally depends on correspondences of point coordinates. An
arbitrary camera, $c$, observes a set of $K$ 3D points $\\{{\bf
X}_{k}\\}_{k=1}^{K}$, and a set of corresponding image points $\\{{\bf
x}_{k}^{c}\\}_{k=1}^{K}$. Point sets and correspondences are known $\forall
k$, and for each camera. In this work, we consider xPnP based camera
calibration using the methods proposed by [21] and [14]. Both extrinsic
(rotation and translation) and intrinsic (focal length and distortion)
parameters in (3) are computed to enable measurement of length and angles in
the camera image.

Distortion Depending on the optical system of a camera, small or large
displacements of image coordinates can be introduced, called image distortion.
Unlike the focal length, which scales the image uniformly, distortion is
characterised as scaling the image differently depending on the distance to a
distortion centre. The farther the pixels are from the centre of distortion,
the more they are distorted. We let the same point represent the distortion
and optical centra, which is assumed to be fixed and in the centre of the
image.

${\bf y^{\prime}}=g({\bf y},\boldsymbol{\theta})$ (1)

We model the distortion as in (1), and let
$\boldsymbol{\theta}=\left[\theta_{1},\theta_{2},\theta_{3}\right]$ specify
the non-linear distortion terms. When $g$ uses one distortion term, it will be
denoted as D(1,0) and as D(3,0) when all three terms are used, according to
[14].

Correspondences The calibration uses Lidar measurements, which map physical
features with high precision by emitting narrow laser beams that are reflected
back. Even if the image, whose camera we want to calibrate, and the lidar map
are recorded at separate times, there should be enough overlapping features
left for calibration. That is, consistent physical properties. Such
properties, which are more likely to be consistent, are, for example, those
found on buildings, vegetation, paintings, furniture, etc. In practice,
correspondences can be of varying quality, making robust estimation a critical
importance, when computing an xPnP solution [14]. Therefore, we use only the
residuals from correspondences marked as inliers by the model estimator for
model validation.

###  3.3 Residual error model

Regardless of whether the corresponding coordinates are found by an interest
point detector, or whether they are manually annotated by a human, they will
suffer from detection noise. This means that a location estimate $\tilde{\bf
x}$ has a residual $\boldsymbol{\epsilon}_{detector}$, compared to the ideal
point location $\hat{\bf x}$. This residual is typically modelled as a 2D
normal distribution:

$\boldsymbol{\epsilon}_{detector}=\tilde{\bf x}-\hat{\bf
x}\sim\mathcal{N}({\bf 0},\sigma_{d}^{2}{\bf I})\,.$ (2)

For a successful calibration, the residual between a detected point, and the
projection of the corresponding 3D point using the estimated parameters,
should also satisfy (2). In other words:

$\hat{\bf x}=\textup{proj}_{\boldsymbol{\Theta},{\bf P}}(\hat{\bf X})={\bf
K}g(\pi({\bf R}\hat{\bf X}+{\bf t}),\boldsymbol{\theta})\,.$ (3)

Here $\hat{\bf X}$ is the ideal 3D point, and $\pi$ is the pinhole projection.
The intrinsic calibration, $\boldsymbol{\Theta}=({\bf
K},\boldsymbol{\theta})$, and the _extrinsic calibration_ , ${\bf P}=({\bf
R},{\bf t})$ (the camera pose) are of course estimates in practice. We
summarize the error caused by the estimation in an additive _modelling noise_
term $\boldsymbol{\epsilon}_{model}$. We intend to explain the residuals by
detection noise in the image and in the Lidar, and test whether the
explanation holds using a test on the residual data,  _e.g._ by the DAP test
[6], or by the KS test [20], testing the GoF.

For 2D-3D matches, the 3D points are also affected by noise
$\boldsymbol{\epsilon}_{3D}$. Thus, (3) should be replaced by:

$\hat{\bf
x}-\boldsymbol{\epsilon}_{lidar}-\boldsymbol{\epsilon}_{model}=\textup{proj}_{{\boldsymbol{\Theta}},{\bf
P}}(\hat{\bf X}-\boldsymbol{\epsilon}_{3D})\,,$ (4)

where $\boldsymbol{\epsilon}_{3D}$ is the detection noise in 3D, and
$\boldsymbol{\epsilon}_{lidar}$ its projection. By combining (4) with (2) we
obtain the following residual model:

$\boldsymbol{\epsilon}=\boldsymbol{\epsilon}_{detector}+\boldsymbol{\epsilon}_{lidar}+\boldsymbol{\epsilon}_{model}\,.$
(5)

We model the detection error as in (2), and describe the lidar error model in
detail below.

####  3.3.1 Lidar error model

For a Lidar sensor, the 3D noise has both angular and depth components.
However, when the camera and 3D-sensor are close to being co-axial, and point
in roughly the same direction (i.e. ${\bf t}$ is small, and ${\bf
R}\approx{\bf I}$ in (3)), the depth error becomes irrelevant, and the
projection in the image $\boldsymbol{\epsilon}_{lidar}$ is dominated by ${\bf
K}$, which is affine. This means that the shape of
$\boldsymbol{\epsilon}_{lidar}$ is a simple, but location dependent scaling.

We thus model the projection of the Lidar error
$\boldsymbol{\epsilon}_{lidar}$ as:

$\boldsymbol{\epsilon}_{lidar}\sim\mathcal{N}(0,\sigma_{l}^{2}\textup{diag}(a_{x}^{2},a_{y}^{2}))\,,$
(6)

where $\sigma^{2}_{l}$ is a noise variance, and $a_{x},a_{y}$ are the noise
scalings in horizontal and vertical directions. These depend on the location
in the image. To estimate $a_{x},a_{y}$, we can project the current 3D point
and its neighbours in pan and tilt directions to obtain:

$\displaystyle a_{x,k}$ $\displaystyle=\|\textup{proj}({\bf
X}_{k})-\textup{proj}({\bf X}_{k}^{P})\|$ (7) $\displaystyle a_{y,k}$
$\displaystyle=\|\textup{proj}({\bf X}_{k})-\textup{proj}({\bf
X}_{k}^{T})\|\,,$ (8)

where ${\bf X}_{k}$ is the current 3D point, and ${\bf X}_{k}^{P}$, and ${\bf
X}_{k}^{T}$ are its neighbours in pan and tilt directions.

###  3.4 Noise Estimation

The parameters of the detector errors in the model (5) can be fitted to the
observed residuals using robust linear regression. However, estimating the
variance of $\boldsymbol{\epsilon}_{model}$ is neglected since its observed
values are those that will remain in order to test whether the model is
incorrect. I.e. we assume:

$E\left\\{\epsilon^{2}\right\\}=\sigma_{d}^{2}+\sigma_{l}^{2}.$ (9)

By using a common $\sigma_{d}$ for x, and y image residuals, and the aspect
ratio model in (6) we obtain:

$E\left\\{\begin{pmatrix}\epsilon_{x}^{2}\\\
\epsilon_{y}^{2}\end{pmatrix}\right\\}=\begin{pmatrix}1&a_{x}^{2}\\\
1&a_{y}^{2}\end{pmatrix}\begin{pmatrix}\sigma_{d}^{2}\\\
\sigma_{l}^{2}\end{pmatrix}.$ (10)

We fit these to the observed $K$ residuals, for each camera, $c$, separately,
to obtain the regressor parameters $(\sigma_{d},\sigma_{l})$.

$\begin{pmatrix}\epsilon_{x,1}^{2}&\epsilon_{y,1}^{2}&\ldots&\epsilon_{x,K}^{2}&\epsilon_{y,K}^{2}\end{pmatrix}^{\intercal}=\begin{pmatrix}1&1&\ldots&1&1\\\
a_{x,1}^{2}&a_{y,1}^{2}&\ldots&a_{x,K}^{2}&a_{y,K}^{2}\end{pmatrix}^{\intercal}\begin{pmatrix}\sigma_{d}^{2}\\\
\sigma_{l}^{2}\end{pmatrix}$ (11)

In practice, we do not use linear regression by solving the normal equations
to (11), but use robust regression using IRLS [11] with initial weights
$1/p(\boldsymbol{\epsilon}_{k}|\sigma=5.0)$. We can now obtain standardised
residuals:

$\tilde{\boldsymbol{\epsilon}}_{k}=\begin{pmatrix}\epsilon_{x,k}/\sigma_{k,x}\\\
\epsilon_{y,k}/\sigma_{k,y}\end{pmatrix}=\begin{pmatrix}\epsilon_{x,k}/\sqrt{\sigma_{d}^{2}+a_{x,k}^{2}\sigma^{2}_{l}}\\\
\epsilon_{y,k}/\sqrt{\sigma_{d}^{2}+a_{y,k}^{2}\sigma^{2}_{l}}\end{pmatrix}.$
(12)

###  3.5 Hypothesis Testing

When the modelling error is low, the standardised residuals in (12) should
pass a statistical test, such as, _e.g._ , the KS test. We can thus use the
test to check whether the calibration worked for a particular set of 2D-3D
correspondences. More formally, we test the null hypothesis $\mathcal{H}_{0}$
: The standardised residuals (12) are distributed explicitly according to a
standard normal, against $\mathcal{H}_{1}$: at least one value does not match
that distribution. Related to this classical approach is that the data we are
testing is random, so the test decision is random too, which means there is
still a small probability of an incorrect decision. Nevertheless, tests are
useful to detect low model errors and thus further validate the calibration.

Evidence The approach involves comparing the samples (residuals) with a
statistical model under $\mathcal{H}_{0}$, where a test statistic measures the
discrepancy between the data and the model. To this end, we use the KS test
[20] and compute a _p-value_ , measuring the error size of rejecting
$\mathcal{H}_{0}$. Commonly, when the _p-value_ is below $5\%$,
$\mathcal{H}_{0}$ can be rejected in favour of $\mathcal{H}_{1}$. That is, the
error probability is sufficiently low. However, this probability does not
directly infer confidence for the data distributed as a standard normal.

Other tests also calculate a _p-value_ to test the normality of data. For
example, the DAP test [6] sums the discrepancies from a skewness test and a
kurtosis test into a single _p-value_. Skewness is the asymmetry about the
mean, and kurtosis is the measure of the ”tailedness”. Although parametric
tests are preferable to non-parametric ones, and the SW test is one of the
more powerful [10], we believe their null hypotheses to be non-directional,
where a broader chance of normality is possible, leading to unstable
decisions.

##  4 Experiments

We first evaluate the proposed method for testing a calibration using 2D-3D
correspondences on synthetic data simulating detection and Lidar errors. Next,
we provide results on a real scenario using Lidar measurements and compare
this with a semi-automatic calibration. Last, we analyse a large-scale dataset
using our method.

Figure 2: Using residuals computed from $56$ sets of simulated 2D-2D
correspondences, we show the results for the incorrect ${\tt D(1,0)}$ model
(left) and the correct ${\tt D(3,0)}$ model (right). We sort the
correspondence sets in ascending order, using the area of the 2D points’
convex hull (_coverage_). The second, third and fourth rows show the outcome
of the KS [20], DAP [6] and SW [26] tests at level $5\%$. When the model is
correct, the standard deviation is low (first row), and our predicted variance
follows the corresponding empirical value. The KS test rejects images under an
incorrect model, while accepting images under the correct model. In contrast,
[6] is too strict, while [26] is too permissive.

Model | $f_{x}$ | $f_{y}$ | $c_{x}$ | $c_{y}$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | std(axis223m) | std(axisp3364)
---|---|---|---|---|---|---|---|---|---
$\textsc{m}_{1}$ | ✓ |  |  |  | ✓ |  |  | 2.99 | 2.91
$\textsc{m}_{2}$ | ✓ |  | ✓ | ✓ | ✓ |  |  | 2.96 | 2.73
$\textsc{m}_{3}$ | ✓ |  |  |  | ✓ | ✓ | ✓ | 3.03 | 2.95
$\textsc{m}_{4}$ | ✓ |  | ✓ | ✓ | ✓ | ✓ | ✓ | 3.19 | 2.52
$\textsc{m}_{5}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | 0.81 | 0.19

Table 1: Projection models used with the Zhang method [30]. Models
$\textsc{m}_{1}$ and $\textsc{m}_{3}$ use the same number of parameters as the
PnP methods D(1,0) [14], D(3,0) [14] and D(3,0) [21]. The standard deviation
on a set of test images, determines how accurately the two cameras have been
calibrated, _axis223m_ and _axisp3364_. $\textsc{m}_{5}$ is the most accurate
model. Models $\textsc{m}_{1}$-$\textsc{m}_{4}$ (approximately) share the same
standard deviation, although models $\textsc{m}_{4}$ and $\textsc{m}_{5}$ only
differ with one parameter. Our method instead decides $\textsc{m}_{1}$ as
incorrect but $\textsc{m}_{3}$ as a plausible model, still usable for
metrology.

###  4.1 Synthetic data

We implemented the simulator in OpenCV [1] and will provide code upon
publication. We aim to render a fictitious checkerboard pattern with
equidistant squares of black and white into a camera $c$ with a small angular
rotation maintaining the image centre as its viewpoint at a distance
${\boldsymbol{t}}$. The pattern contains equally many saddle points (inner
checkerboard corners) vertically as horizontally ($15\times 15$). The
simulator iterates three main tasks to render each image, which is presented
next. (i) The projection model has fixed intrinsic parameters according to
D(3,0) [14]. That is, $f=800$ and
$\boldsymbol{\theta}=\left[-0.0684,0.0100,0.0006\right]$. We let the
distortion centre coincide with the image centre. Rotation parameters are
randomly sampled in the range $\pm 15^{\circ}$ relative to the z-axis. The
translation can also be random, but in our generated synthetic dataset, we
move the pattern closer and closer to the camera. (ii) Next, we smooth the
image (to avoid aliasing), add image noise and interpolate it to size
$1600\times 1600$. A saddle point detector [1] locates the 2D position of
these features with sub-pixel precision. We observed the position error to lie
within a small range of $0.03$ pixels. This corresponds to $\sigma_{d}$ in
(9).(iii) To simulate the Lidar, we add noise on the corresponding 3D points
in all images. We transform the noise such that it lives on the sphere with
origin ${\boldsymbol{t}}^{c}$ and radius $||{\boldsymbol{t}}^{c}||$. On the
sphere, the noise magnitude is dependent on ${\boldsymbol{t}}^{c}$, and in the
vertical direction, the noise is always $90\%$ lower compared to the
horizontal direction. This reweighing aims to simulate the resolution in the
Lidar array, and it replaces (7) and (8), which are used on real datasets.
This gives us the resolution aspect, which is used as _explanatory variables_
in the estimate, (11).

Figure 3: Incorrect models can have small residual errors, but _Godness-of-
Fit_ testing exposes them. Left: A distorted image from which the model is
estimated. Middle: Undistorted image using true model (black) and wrong model
(grey) overlaid. Even though the estimated model is incorrect, when detector
errors are small the residual errors are often small, hence simply checking
the standard deviation of the residuals is insufficient. Right: Outliers
indicate model failure. Given known outlier-free correspondences, deviations
from the expected noise distribution expose incorrect models.

Results In Figure 2, we show the output of our method on models D(1,0) [14]
and D(3,0) [14] for a set of synthetic 2D-3D correspondences, generated under
the D(3,0) model. These are sorted according to the increasing spatial spread
of correspondences. The closer the points are to the image edge, the more they
are affected by the distortion (Figure 3). The first row in Figure 2 shows the
empirical and estimated standard deviations over all points per image, and
rows 2-4 whether the scaled residuals are sufficient evidence to reject the
null hypothesis for three different tests. Given in the second row of Figure 2
is the output of the KS test for both D(1,0) [14] and D(3,0) [14]. To the
left, in the same row, our method rejects the images with correspondences more
uniformly spread over the entire image at level $5\%$ using D(1,0) [14].
Making the same test using a projection model with more parameters fits the
data more accurately, as indicated by both the low standard deviation and the
test. We also test the scaled residuals in the third and fourth rows using
[6], and [26], respectively. While these tests are parametric, they test for
any normal distribution. Compared to the proposed KS test, this leads to false
positives and negatives, see rows 3-4 of Figure 2. In the second column, the
stronger D(3,0) projection model is tested. In most images, $\mathcal{H}_{0}$
can not be rejected as expected due to the simulated data conforming to D(3,0)
[14]. However, the output does not reveal the tests’ differences in this case.
Thus, [20] generally leads to more accurate decisions.

PnP Method | DT | $\sqrt{\textrm{Var}\left[\boldsymbol{\epsilon}\right]}$ | KS[20] | DAP[6] | SW[26]
---|---|---|---|---|---
EAPRD [14] | D(1,0) | 2.23 | ✓ | ✓ |
EAPRD [14] | D(3,0) | 0.98 |  | ✓ |
PNPRF [21] | D(3,0) | 1.05 |  | ✓ |

Table 2: The results for the _axis223m_ camera. For each PnP method, with
distortion type DT, we report the empirical standard deviation and whether
$\mathcal{H}_{0}$ can be rejected at level $5\%$ using the KS, DAP or SW test.

PnP Method | DT | $\sqrt{\textrm{Var}\left[\boldsymbol{\epsilon}\right]}$ | KS[20] | DAP[6] | SW[26]
---|---|---|---|---|---
EAPRD [14] | D(1,0) | 1.71 | ✓ | ✓ |
EAPRD [14] | D(3,0) | 1.47 |  |  |
PNPRF [21] | D(3,0) | 1.08 |  |  |

Table 3: The results for the _axisp3364_ camera. For each PnP method, with
distortion type DT, we report the empirical standard deviation and whether
$\mathcal{H}_{0}$ can be rejected at level $5\%$ using the KS, DAP or SW test.

###  4.2 Lidar Measurements

Next, we compare our method using images from two real cameras, _axis223m_ and
_axisp3364_ , respectively, and a 3D point cloud from a _Leica RTC360_
scanner, with semi-automatic calibration. The second camera offers lower-
quality images than the first, which is visible in Figure 4. There is also no
verified annotation for the cameras; in practice, there is none, and the
camera can be inaccessible. The cameras instead depict a scene such that their
optical axes have a relatively small angle to the 3D point cloud coordinate
system’s z-axis, similar to the simulations in Section 4.1.

Semi-automatic For semi-automatic calibration, we collect images for both
cameras using a checkerboard pattern. The pattern has $6\times 6$ saddle
points. We split the set of images and estimate the model parameters on the
first set using the Zhang method [30] implemented in [1]. Then, we calibrate
using different projection models where $\textsc{m}_{5}$, in Table 1, achieves
the lowest standard deviation on the second set of images (test). A factor of
almost 4 differs between the accuracy of models $\textsc{m}_{4}$ and
$\textsc{m}_{5}$. In Table 1, the number of model parameters differs by one.

Without Camera Access We show in Tables 2 and 3 that D(1,0) [14] is incorrect
compared to D(3,0) [14] and D(3,0) [21] using the KS test on 24 manually
annotated correspondences. Each of the annotated 3D points is visible in both
cameras and projected to consistent features,  _e.g._ corners. Similar to
simulation, we observe that the parametric tests contradict each other and are
thus infeasible for our application. While models $\textsc{m}_{1}$ to
$\textsc{m}_{4}$, in Table 1, obtain higher standard deviations using [30],
models $\textsc{m}_{1}$ and $\textsc{m}_{3}$ are equvivalent to the models
used from [14] and [21], and thus there is possibility that $\textsc{m}_{5}$
is overfitted.

Finally, we found that the computed distortion parameters of D(3,0) [21] were
all zero, shown in the rightmost column of Figure 4. To our knowledge, [21]
divides the xPnP problem into subproblems. In the subproblem that solves
distortion, we can’t find a condition on $\boldsymbol{\theta}$ preventing the
_normal equations_ from giving the trivial solution. Thus, our method can not
make the correct decision to either reject $\mathcal{H}_{0}$ or not based on
residuals from D(3,0) [21].

Figure 4: The first column shows two images of two cameras, _axis223m_ and
_axisp3364_ , respectively. In the second and third columns, the undistortion
looks to be visually removed for both D(1,0) [14], and D(3,0) [14]. Our method
correctly detects D(1,0) [14] as incorrect for both cameras (Tables 2 and 3).
The undistorted images in the fourth column are visually similar to their
original, but this is not detected. For more details, see Section 4.2.

###  4.3 Structure-from-Motion

In this experiment, we use our proposed method on annotations computed from a
Structure-from-Motion (SfM) pipeline to get a broader insight into its
effectiveness. To this end, we use 1000 images from each scene of MegaDepth
[18]. This dataset contains many scenes with 2D-3D correspondences and camera
intrinsic and extrinsic parameters given. The SfM pipeline, COLMAP [24][25],
estimates the annotation parameters of the widely used benchmark for state-of-
the-art comparison. In the dataset, the assumed projection model, a _simple
radial_ , models a single focal length, one distortion parameter, the
distortion centre, rotation and translation. The histogram to the left in
Figure 5 shows that residuals are overall low. However, in 70 out of 100
images, our method rejects the null hypothesis at level $5\%$. The two images
in the middle and to the right, in Figure 5, show when $\mathcal{H}_{0}$ can
not be rejected at level $5\%$ and when $\mathcal{H}_{0}$ is rejected in
favour of $\mathcal{H}_{1}$. We can thus assume mostly overfitted projection
models in [18].

Figure 5: Left: Density plot of residuals from 1000 images in all scenes, on
which the annotation in MegaDepth depends. It is unlikely residuals will be
high for images in [18] measuring a good performance. However, our proposed
method tests each image and rejects the null hypothesis, $\mathcal{H}_{0}$, on
70 out of 100 images. Middle: Example of when $\mathcal{H}_{0}$ can not be
rejected, and the _simple radial_ projection model is suitable. Right: Example
of when our method rejects $\mathcal{H}_{0}$. As can be seen,  _e.g._ on the
flagpole to the right, the images are distorted.

##  5 Conclusion

We suggested that metrology applications in forensic analysis use xPnP methods
and use our proposed method to validate the calibration without camera access.
The method formulation processes a single image, estimating a robust scaling
of each correspondence and tests if the scaled set of residuals is drawn from
a standard normal distribution. We demonstrate via qualitative and
quantitative experiments that the KS test is most suitable and provide further
insight from an extensive collection of annotated cameras.

Although we are sufficiently confident that the test can determine models as
incorrect with a small margin of error, the challenge remains to infer
confidence in the image measurements. A test is not a classification, and the
_p-value_ does not imply measurement confidence. However, when rejection of
the null hypothesis is not possible at the acceptable error level, our error
model explicitly provides the expected measurement errors over the image.
Depending on the number of correspondences, we can get local estimates of
expected measurement error from our assumptions of normally distributed
residuals. Therefore, our method is a useful tool for xPnP camera calibration.

### Acknowledgement

_This work was partially supported by the Wallenberg AI, Autonomous Systems,
and Software Program (WASP) funded by the Knut and Alice Wallenberg
Foundation; and the computations were enabled by the Berzelius resource
provided by the Knut and Alice Wallenberg Foundation at the National
Supercomputer Centre; and a point cloud of a realistic scene was provided by
the Swedish National Forensic Centre (NFC)._

## References

  * [1] Bradski, G.: The OpenCV Library. Dr. Dobb’s Journal of Software Tools (2000)
  * [2] Bramble, S., Compton, D., Klasén, L.: Forensic image analysis. In: Proceedings of the 13th INTERPOL Forensic Science Symposium (2001)
  * [3] Brandner, M.: Bayesian uncertainty evaluation in vision-based metrology. In: Gallegos-Funes, F. (ed.) Vision Sensors and Edge Detection, chap. 5. IntechOpen, Rijeka (2010). https://doi.org/10.5772/10135, https://doi.org/10.5772/10135
  * [4] Criminisi, A.: Single-view metrology: Algorithms and applications. In: Pattern Recognition, 24th DAGM Symposium. Zurich, Switzerland (January 2002)
  * [5] Criminisi, A., Reid, I., Zisserman, A.: Single view metrology. International Journal of Computer Vision 40(2), 123–148 (2000)
  * [6] D’Agostino, R., Pearson, E.S.: Tests for departure from normality. Biometrika 60, 613–622 (1973)
  * [7] Fischler, M.A., Bolles, R.C.: Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. Commun. ACM 24(6), 381–395 (jun 1981)
  * [8] Förstner, W.: Uncertainty and projective geometry. In: Handbook of Geometric Computing, pp. 493–534. Springer (2005)
  * [9] Gao, X.S., Hou, X.R., Tang, J., Cheng, H.F.: Complete solution classification for the perspective-three-point problem. IEEE Transactions on Pattern Analysis and Machine Intelligence ( Volume: 25, Issue: 8, Aug. 2003) (2003)
  * [10] Ghasemi, A., Zahediasl, S.: Normality tests for statistical analysis: a guide for non-statisticians. Int. J. Endocrinol. Metab. 10(2), 486–489 (April 2012)
  * [11] Green, P.J.: Iteratively reweighted least squares for maximum likelihood estimation, and some robust and resistant alternatives. Journal of the Royal Statistical Society. Series B (Methodological) 46(2), 149–192 (1984)
  * [12] Hartley, R., Zisserman, A.: Multiple View Geometry in Computer Vision, 2nd ed. Cambridge University Press (2003)
  * [13] Heuel, S.: Uncertain projective geometry: statistical reasoning for polyhedral object reconstruction, vol. 3008. Springer (2004)
  * [14] Larsson, V., et al.: Revisiting radial distortion absolute pose. International Conference on Computer Vision (ICCV) (2019)
  * [15] Lebeda, K., Matas, J., Chum, O.: Fixing the locally optimized ransac. British Machine Vision Conference (BMVC) (2012)
  * [16] Lepetit, V., Moreno-Noguer, F., Fua, P.: Epnp: An accurate o(n) solution to the pnp problem. International Journal Of Computer Vision 81, 155–166 (2009)
  * [17] Li, X., Zhang, B., Sander, P.V., Liao, J.: Blind geometric distortion correction on images through deep learning. Conference on Computer Vision and Pattern Recognition (CVPR) (2019)
  * [18] Li, Z., Snavely, N.: Megadepth: Learning single-view depth prediction from internet photos (2018)
  * [19] Lu, X.X.: A review of solutions for perspective-n-point problem in camera pose estimation. Journal of Physics: Conf. Ser. 1087 052009 (2018)
  * [20] Massey, F.J.: The kolmogorov-smirnov test for goodness of fit. Journal of the American Statistical Association 46(253), 68–78 (1951)
  * [21] Nakano, G.: A versatile approach for solving pnp, pnpf and pnpfr problems. European Conference on Computer Vision (ECCV) (2016)
  * [22] Olsson, E.: Lens Distortion Correction Without Camera Access. Master’s thesis, Linköping University, Sweden (2022)
  * [23] Persson, M., Nordberg, K.: Lambda twist: An accurate fast robust perspective three point (p3p) solver. In: Proceedings of the European Conference on Computer Vision (ECCV) (September 2018)
  * [24] Schönberger, J.L., Frahm, J.M.: Structure-from-motion revisited. In: Conference on Computer Vision and Pattern Recognition (CVPR) (2016)
  * [25] Schönberger, J.L., Zheng, E., Pollefeys, M., Frahm, J.M.: Pixelwise view selection for unstructured multi-view stereo. In: European Conference on Computer Vision (ECCV) (2016)
  * [26] Shapiro, S.S., Wilk, M.B.: An analysis of variance test for normality (complete samples)†. Biometrika 52(3-4), 591–611 (1965). https://doi.org/10.1093/biomet/52.3-4.591
  * [27] Thai, T.H., Cogranne, R., Retraint, F.: Camera model identification based on the heteroscedastic noise model. IEEE Transactions on Image Processing 23(1), 250–263 (2014)
  * [28] Tsai, R.Y.: A versatile camera calibration technique for high-accuracy 3d machine vision metrology using off-the-shelf tv cameras and lenses. IEEE (1987)
  * [29] Wang, Q., Alexander, W., Pegg, J., Qu, H., Chen, M.: Hypoml: Visual analysis for hypothesis-based evaluation of machine learning models. IEEE Transactions on Visualization and Computer Graphics 27(2), 1417–1426 (2021)
  * [30] Zhang, Z.: A flexible new technique for camera calibration. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(11), 1330–1334 (2000)


# Blind estimation of audio effects using an auto-encoder approach and
differentiable signal processing

Côme Peladeau and Geoffroy Peeters

(LTCI, Télécom-Paris, Institut Polytechnique de Paris, France)

###### Abstract

Blind Estimation of Audio Effects (BE-AFX) aims at estimating the Audio
Effects applied to an original, unprocessed audio sample solely based on the
processed audio sample. To train such a system traditional approaches optimize
a loss between ground truth and estimated AFX parameters. This involves
knowing the exact implementation of the AFXs used for the process. In this
work, we propose an alternative solution that eliminates the requirement for
knowing this implementation. Instead, we introduce an auto-encoder approach,
which optimizes an audio quality metric. We explore, suggest, and compare
various implementations of commonly used mastering AFXs, using differential
signal processing or neural approximations. Our findings demonstrate that our
auto-encoder approach yields superior estimates of the audio quality produced
by a chain of AFXs, compared to the traditional parameter-based approach, even
if the latter provides a more accurate parameter estimation.

##  1 Introduction

Audio Effects (AFXs) play an essential role in the music production. They are
used during mixing to sculpt sounds for artistic purposes or context
requirements (such as when a sound needs to be mixed with others). They are
used during mastering, the final stage of production, to improve the clarity
of a given mix, adapt it for a given media (such as Vynil or streaming), or
harmonize it with other tracks. For these reasons, its automatization has been
the subject of several softwares111such as Izotope Ozone 11 or FabFilter Pro Q
3 which allows learning the mastering of a given track to apply it to new dry
mixes. However, those only focus on the equalization. In this work, we study
the generalization to other common mastering AFXs.

Blind Estimation of Audio Effects (BE-AFX) aims at estimating the Audio
Effects (AFXs) applied to an original, unprocessed (dry) audio sample
$\mathbf{x}$ solely based on the observation of the processed (wet) audio
sample $\mathbf{y}$. This estimation takes the form of the AFXs and their
parameters $\mathbf{p}$.

###  1.1 Related works.

For long BE-AFX techniques has been based on explicit rules and assumptions.
For example, Ávila et al. [1] proposed to estimate the curve of a memoryless
non-linear distortion by assuming that the unprocessed signal has the
statistics of a Gaussian white noise. However, nowadays, most BE-AFX
approaches rely on training neural networks. Indeed, following SincNet [2] and
DDSP [3], modeling audio processes as differential processes has allowed
developing differentiable AFXs as specialized neural networks layers with
interpretable parameters [4, 5, 6]. Because they are differentiable, they can
be integrated transparently in a neural network. Since then, differentiable
AFXs have been used for many tasks: automatic mixing and mastering [7, 8],
production style transfer [6] or estimation of audio effects [9].

In the case of BE-AFX, neural networks are usually trained to minimize a loss
function that aims at reconstructing the AFX ground truth parameters [10, 11].
However, as we will highlight in this work, a parameter distance does not
translate well to a perceptual distance between audio effects. This is why we
will propose the use of an audio distance here.

Estimation of audio effects with an audio loss function and differentiable
audio effects has already been investigated. For example, Colonel et al. [12]
used it for non-linear distortion using a differentiable Wiener-Hammerstein
model and also in [9] for a complete mixing setting. However, in both cases,
their approaches require paired $\mathbf{x}$ and $\mathbf{y}$ data for the
estimation, i.e. they did not perform the blind estimation. In this work, we
perform blind estimation, i.e. we aim at estimating the AFXs applied to
$\mathbf{x}$ using only the knowledge of $\mathbf{y}$.

Figure 1: Proposed auto-encoder approach for Blind Estimation of Audio Effects
(BE-AFX).

###  1.2 Proposal.

To solve the BE-AFX problem, we propose an auto-encoder approach which is
illustrated in Figure 1. In the left part, we construct synthetic processed
mixes by applying a set of (synthesis) audio effects $\\{e^{s}\\}$ with known
parameters $\mathbf{p}$ to an unprocessed mix $\mathbf{x}$. The results are
our ground-truth processed mixes $\mathbf{y}$. Using only $\mathbf{y}$, an
analysis network $f^{a}$ then estimates the set of parameters
$\hat{\mathbf{p}}$ to be used to process $\mathbf{x}$ with (analysis) audio
effects $\\{e^{a}\\}$ to produce an estimated audio sample $\hat{\mathbf{y}}$.
The analysis network $f^{a}$ is trained to minimize an audio loss function
between $\hat{\mathbf{y}}$ and $\mathbf{y}$ so that
$\hat{\mathbf{y}}\approx\mathbf{y}$. It therefore closely matches the
formulation of an auto-encoder.

Doing so, when $\\{e^{a}\\}$=$\\{e^{s}\\}$, $f^{a}$ implicitly learns to
replicate the parameters $\mathbf{p}$ given only $\mathbf{y}$. When
$\\{e^{a}\\}\neq\\{e^{s}\\}$, $f^{a}$ learns parameters to be used for
$\\{e^{a}\\}$ such that the effect of the analysis chain sounds similar to the
synthesis chain.

###  1.3 Paper organization.

To be able to estimate $\mathbf{p}$ using this auto-encoder, we should be able
to “differentiate” the audio effects $\\{e^{a}\\}$ (i.e. compute the
derivative of their outputs $\hat{\mathbf{y}}$ w.r.t. their inputs
$\hat{\mathbf{p}}$). We therefore discuss various implementations of the audio
effects in part 2.1. To be able to estimate $\mathbf{p}$ we should define an
architecture for the analysis network $f^{a}$. We therefore discuss various
possible architectures in part 2.2. During evaluation (part 3), we first
decide for each type of effect what is the best implementation $e^{a}$ (among
those of 2.1) and architecture of $f^{a}$ (among those of 2.2) to reconstruct
$\mathbf{y}$. We then compare our proposed approach (audio reconstruction
$\hat{\mathbf{y}}\approx\mathbf{y}$) to the previously proposed approach
(parameter reconstruction $\hat{\mathbf{q}}\approx\mathbf{q}$). Finally, we
evaluate the joint estimation of the whole AFX chain defined as the succession
of an equalizer, a compressor, and a clipper. We conclude in part 4 and
propose future directions.

To ensure replicability, we provide the code of this study.222https://
github.com/ peladeaucome/ ICASSP-2024-BEAFX-using-DDSP.git

##  2 Proposal

###  2.1 Audio effects

In this work, we only consider the 3 following AFXs commonly used for
mastering: equalizer, dynamic range compressor, and clipper. We distinguish
their implementation for synthesis $\\{e^{s}\\}$ and for analysis
$\\{e^{a}\\}$.

$\\{e^{s}\\}$ is the implementation of the effects that have been used to
create the observed master $\mathbf{y}$. In a real scenario, this
implementation is unknown.

$\\{e^{a}\\}$ is the implementation we use in our model to (a) predict the
parameters $\hat{\mathbf{p}}$ of the effects or (b) replicate the resulting
process of the mastering. To perform (a) (comparing the estimated
$\hat{\mathbf{p}}$ to the ground-truth $\mathbf{p}$) the implementation of the
effect in $\\{e^{s}\\}$ and $\\{e^{a}\\}$ should be similar. To perform (a)
and (b), the implementation in $\\{e^{a}\\}$ should be differentiable (in
order to be able to estimate the parameters) or, if not, we will have to use
neural networks to approximate the effects..

We now detail the implementation of the three AFXs for synthesis and analysis
and list them in Table 1.

The audio is normalized to 0 dBFS before each effect.

Table 1: Considered audio effects and their implementations. Effect | Implementation | Synthesis | Analysis
---|---|---|---
Equalizer | Parametric | $\surd$ | $\surd$
| Graphic |  | $\surd$
Compressor | DSP | $\surd$ |
| Simplified DSP |  | $\surd$
| NP |  | $\surd$
|  Hybdrid NP |  | $\surd$
Clipper | Parametric | $\surd$ | $\surd$
| Taylor |  | $\surd$
| Chebyshev |  | $\surd$

####  2.1.1 Equalizer

For synthesis, we use a 5-band parametric equalizer: 1 low-shelf, 3 peak, 1
high-shelf.

For analysis, we use either this parametric equalizer or a 10 band graphic
equalizer. Each band of a graphic equalizer has a bandwidth of 2 octaves [13].

Each parametric band has 3 parameters: center frequency, gain, and quality
factor. Each band of the graphic equalizer has only 1 parameter: the gain.

Frequencies of parametric bands are logarithmic parameters, gains (in dB) and
quality factors are linear. Differentiable filters are implemented in the
frequency domain as we find the time-aliasing error small enough for training
a neural network.

####  2.1.2 Dynamic range compressor

For synthesis, we use the DSP compressor proposed in [14]. It has 5
parameters: threshold, ratio, attack time, release time, and knee width.

For analysis, we either use a simplified DSP compressor or a Neural Proxy
(NP). The simplified DSP compressor is the DSP compressor of [14] but with the
attack and release time linked, has proposed by [6] to reduce the computation
time. The NP compressor [15] is trained to approximate a DSP compressor333 To
train it, we first process a set of $\mathbf{x}$ with the DSP compressor using
known parameters $\mathbf{p}$. The output provides ground-truths $\mathbf{y}$.
We then train the NP compressor conditioned with the same parameters
$\mathbf{p}$ such that its output $\hat{\mathbf{y}}\approx\mathbf{y}$.. It
uses a TCN conditioned with FiLM [16] layers. In [15] the NP directly outputs
$\hat{\mathbf{y}}$. In our case, we use the same TCN architecture but propose
to replace its output activation by a sigmoid such that it provides the
compressor gain factor $\mathbf{g}$ to be applied over time $n$:
$\hat{y}[n]=g[n]\cdot x[n]$. 444 We found by experiment that this modification
allows to largely reduce the number of parameters (number of TCN channels)
with equivalent performances. With 8 channels, our causal model with a
receptive field of 3 s duration has a test mean average error (MAE) of 0.0060
while the TCN from [15] has a test MAE of 0.050.

Once trained, the NP compressor, being differentiable, can be inserted in the
analysis chain $\\{e^{a}\\}$ to train the analysis network $f^{a}$ and obtain
its compressor parameters $\hat{\mathbf{p}}$. We can of course use
$\hat{\mathbf{p}}$ to process $\mathbf{x}$ with the NP compressor but also use
$\hat{\mathbf{p}}$ to process $\mathbf{x}$ with the DSP compressor. We name
the latter “hybrid NP compressor”. Since it is not differentiable, we only use
it during validation and testing (not during training). It has already been
used in [6].

The compressors’ ratio and time parameters are logarithmic while their
threshold and knee (both in dB) are linear.

####  2.1.3 Clipper

We propose 3 implementations of the clipper: parametric (both for synthesis
and analysis), Taylor, and Chebyshev (only for analysis).

The parametric clipper is defined by the function $f$ defined in (1) with
hardness parameter $h$ which blends between tanh, cubic and hard clipping:

$f(x,h)=\begin{cases}(1-h)\tanh(x)+hf_{\text{cubic}}(x),&h\in[0;1],\\\
(2-h)f_{\text{cubic}}(x)+(h-1)f_{\text{hard}}(x),&h\in[1;2].\end{cases}$ (1)

with :

$\displaystyle f_{\text{hard}}(x)=\max(-1,\min(1,x))$ (2) $\displaystyle
f_{\text{cubic}}(x)=\begin{cases}x+4x^{3}/27,&x\in[-\frac{2}{3};\frac{2}{3}],\\\
\operatorname{sign}(x),&|x|>\frac{2}{3}.\end{cases}$ (3)

The effect used for synthesis is constructed with the following parameters:
gain $g$ (in dB), offset $o$, and hardness $h$.

$y[n]=\left(f(g\cdot x[n]+o,h)-f(o,h)\right)/g$ (4)

We also use two other memory-less models. Both have been proposed for memory-
less distortion identification. The Taylor clipper is inspired by Taylor
series expansions [1]:

$y[n]=\sum_{h=0}^{H-1}g_{h}x[n]^{h}$ (5)

The Chebyshev clipper is inspired by Chebyshev’s polynomials as used for non-
linear audio effect identification [17]:

$y[n]=\sum_{h=0}^{H-1}g_{h}T_{h}(y[n]).$ (6)

with $g_{h}\in[-1;1]$ being the effect’s parameters and :

$\displaystyle T_{n}(x)=2xT_{n-1}(x)-T_{n-2}(x),\qquad\forall n\geq 2,$ (7)
$\displaystyle T_{0}(x)=1,\qquad T_{1}(x)=x.$

In both cases, the parameters to be estimated are the $\\{g_{h}\\}$ and we set
$H=24$. All clipper parameters are linear.

####  2.1.4 Parameter ranges

All AFX parameters, $p_{c}/\hat{p}_{c}$, are derived from “normalized” AFX
parameters, $q_{c}/\hat{q}_{c}$ $\in[0,1]$. Linear parameters (see above) are
derived from $q_{c}/\hat{q}_{c}$ using an affine transformation:

$p_{c}=(p_{c,\text{max}}-p_{c,\text{min}})q_{c}+p_{c,\text{min}}.$ (8)

Logarithmic parameters (see above) are derived from $q_{c}/\hat{q}_{c}$ using
an exponential transformation:

$p_{c}=e^{\left(\log\left(p_{c,\text{max}}\right)-\log\left(p_{c,\text{min}}\right)\right)q_{c}}p_{c,\text{min}}.$
(9)

In all cases, $p_{c,\text{min}}/p_{c,\text{max}}$ are the extreme values of
the range of parameter $c$.

###  2.2 Analysis network $f^{a}$

The analysis network is divided into two parts:

  1. 1.

An encoder; which outputs a time invariant embedding. We compare 3
implementations of the encoder described below.

  2. 2.

A MLP with 4 layers of size 2048, 1024, 512, and $C$ where $C$ is the total
number of parameters $\hat{\mathbf{p}}$ of the effect chain $\\{e^{a}\\}$.
Each hidden layer is followed by a Batchnorm1d and a PReLU. The output layer,
which estimates normalized parameters $\hat{\mathbf{q}}$, is followed by a
sigmoid.

For the encoder, we compared three architectures proposed in the literature:

MEE

the Music Effects Encoder proposed by [18]. It consists of a cascade of
residual 1D convolutional layers,

TE

a Timbre Encoder inspired by [19]. It consists of a single 2d convolution
layer with multiple sizes of filters. The conv2d is applied on the CQT [20] of
the signal,

TFE

a Time+Frequency Encoder inspired by [21]. It consists of two 2d convolutional
nets, one focusing on highlighting temporal motifs, the second frequency
motifs. The network’s input is the CQT of the signal.

MEE has 88M parameters, TE 2.8M, and TFE 3.4M.

##  3 Evaluation

###  3.1 Dataset

For training, validation and testing we used the mix files of MUSDB18 [22].
From those, we extracted randomly clips of 10 s duration. Those are then
converted to mono and peak-normalized to 0 dBFS. For each, we randomly pick
the normalized parameters $\mathbf{q}\sim\mathcal{U}(0,1)$, convert them to
$\mathbf{p}$ and apply them to the clip. We used the training, validation and
testing splits proposed by MUSDB18 [22].

###  3.2 Training

Each model $f^{a}$ is trained using the ADAM algorithm with a learning rate of
$10^{-4}$ and a batch size of 16 during a maximum of 400 epochs, where a
single epoch is defined as 430 training examples (5 from each song of the
training subset). The learning rate is scheduled to decrease by a factor of 10
when the best validation score has not been improved for 30 epochs. Training
stops after 150 epochs without improvement. To ensure the reliability of
computed scores the validation is run 5 times and the test 10 times.

In the following, we compare 2 approaches for training $f^{a}$:

Audio reconstruction $\hat{\mathbf{y}}\approx\mathbf{y}$

(our proposal): we minimize the $\ell^{1}$ norm between the log-magnitude Mel-
spectrograms of $\mathbf{y}$ and $\hat{\mathbf{y}}$ [23]:

$\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}=\left\|\log(|\text{Mel}(\hat{\mathbf{y}}|)-\log(|\text{Mel}(\mathbf{y}|)\right\|_{1}.$
(10)

Parameter reconstruction $\hat{\mathbf{q}}\approx\mathbf{q}$

(previously proposed approach): we minimize the $\ell^{2}$ norm between
$\hat{\mathbf{q}}$ and $\mathbf{q}$:

$\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}{}=\frac{1}{C}\sum_{c=0}^{C-1}(\hat{q}_{c}-q_{c})^{2}.$
(11)

###  3.3 Performance metrics.

For evaluation, we compute the following metrics:

  * •

$\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$: the MSE between $\hat{\mathbf{y}}$
and $\mathbf{y}$

  * •

$\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$as defined above

  * •

$\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$: the MSE between $\hat{\mathbf{q}}$
and $\mathbf{q}$

To make audio loss and metrics independent of sound level, we normalize
$\mathbf{y}$ and $\hat{\mathbf{y}}$ by their respective RMS values.

###  3.4 Results

####  3.4.1 Single effect estimation.

For each type of effect we first decide what is the best implementation (among
those of 2.1) and architecture (among those of 2.2) to reconstruct
$\mathbf{y}$. For each, we also indicate: \- Random $\hat{\mathbf{q}}$: the
results obtained with a random choice of $\hat{\mathbf{q}}$ (rather than the
estimated one), \- $\mathcal{L}(\mathbf{x},\mathbf{y})$: the value of the loss
when comparing the input $\mathbf{x}$ to the output $\mathbf{y}$.

For the equalizer (Table 4), Parametric and Graphic provide similar results.
Since $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ indicates the
difference between spectras, it is more suited to measure the performances of
an EQ than $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$. We therefore focus on
$\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ and conclude that the
best (0.32) configuration is to use the Parametric implementation for
$\\{e^{a}\\}$ and TFE for $f^{a}$.

For the compressor (Table 4), the best (0.011, 0.076) configuration is to use
the Hybrid NP for $\\{e^{a}\\}$ and use MEE for $f^{a}$. As a reminder, the
Hybrid NP compressor uses the NP compressor to estimate $\hat{\mathbf{p}}$ but
the DSP compressor (the same used for synthesis) to get $\hat{\mathbf{y}}$.
This works better than using the NP compressor (0.014, 0.098) or the
simplified DSP compressor (0.041, 0.16) This is due to the fact that the
latter links the attack and release time parameters which might be too
restrictive. The fact that the Hybrid NP works better indicates that our proxy
performs well enough for the task of estimated parameters usable for the DSP
Compressor.

For the clipping (Table 4), the best configuration is to use the Parametric
clipper (the same used for synthesis) for $\\{e^{a}\\}$ and use MEE for
$f^{a}$.

It should be noted that in all three cases, the best results are obtained when
$\\{e^{a}\\}$=$\\{e^{s}\\}$.

Table 2: Results by implementation of $\\{e^{a}\\}$ and encoder $f^{a}$ for
equalisation. The best results are indicated in bold.

Equalizer | Parametric | Graphic
---|---|---
Metrics | $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ | $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ | $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ | $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$
MEE | 0.35 | 0.41 | 0.38 | 0.41
TE | 0.35 | 0.46 | 0.39 | 0.48
TFE | 0.28 | 0.32 | 0.25 | 0.34
Random $\hat{\mathbf{q}}$ | 0.86 | 0.89 | 1.1 | 1.1
$\mathcal{L}(\mathbf{x},\mathbf{y})$ | 0.48 | 0.64 | 0.48 | 0.63

Table 3: Same for dynamic range compression.

Comp. | NP |  Hybrid NP | Simpl. DSP
---|---|---|---
Metrics |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$
MEE |  0.014 |  0.098 |  0.011 |  0.076 |  0.041 |  0.16
TE |  0.020 |  0.12 |  0.019 |  0.11 |  0.042 |  0.17
TFE |  0.015 |  0.11 |  0.012 |  0.086 |  0.042 |  0.17
Random $\hat{\mathbf{q}}$ |  0.038 |  0.16 |  0.036 |  0.15 |  0.13 |  0.66
$\mathcal{L}(\mathbf{x},\mathbf{y})$ |  0.041 |  0.16 |  0.041 |  0.16 |  0.040 |  0.16

Table 4: Same for clipping

Clipper | Parametric | Taylor | Chebyshev
---|---|---|---
Metrics |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$
MEE |  $\mathbf{0.0045}$ |  $\mathbf{0.072}$ |  3.4 |  0.15 |  0.65 |  0.17
TE |  0.011 |  0.090 |  2.9 |  0.21 |  3.2 |  0.15
TFE |  $0.0067$ |  0.075 |  2.1 |  0.16 |  1.34 |  0.12
Random $\hat{\mathbf{q}}$ |  0.078 |  0.29 |  1.9 |  0.39 |  2.0 |  1.1
$\mathcal{L}(\mathbf{x},\mathbf{y})$ |  0.078 |  0.34 |  0.079 |  0.34 |  0.079 |  0.34

Table 5: Comparison between training using $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ and $\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$ for single effects and whole FX chain. Loss | $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ | $\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$
---|---|---
Metrics |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$ |  $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$ |  $\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$
Eq. |  0.28 |  0.32 |  0.089 |  0.27 |  0.40 |  0.072
Comp. |  0.011 |  0.076 |  0.11 |  0.0081 |  0.069 |  0.069
Clip. |  0.0045 |  0.072 |  0.044 |  $\mathbf{0.0041}$ |  0.064 |  0.028
Chain |  0.31 |  0.40 |  0.10 |  0.33 |  0.49 |  0.072

####  3.4.2 Training method comparison.

We now compare our proposed training method (based on audio reconstruction
$\hat{\mathbf{y}}\approx\mathbf{y}$) to the previously proposed one (based on
parameter reconstruction $\hat{\mathbf{q}}\approx\mathbf{q}$). Results are
indicated in Table 5. For each single effect, we use the best configuration
found above: $f^{a}$=TFE for equalizer, MEE for compression and clipping.

For equalization, in terms of audio quality
($\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$), the network
trained to minimize $\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$
outperforms (0.32) the one that minimizes
$\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$ (0.40). But in terms of parameter
estimation ($\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$) minimizing directly
$\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$ leads to better results (0.072).

For compression and clipping, training by minimizing the parameter distance
($\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$) leads to better results both in
terms of audio quality
($\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$=0.069, 0.064) and
parameter estimation ($\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$=0.069,
0.028).

####  3.4.3 Effects chain estimation

We finally evaluate the estimation of the whole chain of effect
(equalizer+compressor+clipper). In this case, we use $f^{a}$=TFE for all. As
far as concerned audio quality we see that our approach (minimizing
$\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$) leads to the best
results: $\text{MSE}_{\hat{\mathbf{y}},\mathbf{y}}$=0.31 and
$\mathcal{L}^{\text{Mel}}_{\hat{\mathbf{y}},\mathbf{y}}$=0.40. However, as can
be predicted, minimizing directly $\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$
leads to better estimation of the parameters:
$\text{MSE}_{\hat{\mathbf{q}},\mathbf{q}}$=0.072. These contrasting outcomes
underscore that achieving accurate parameter estimation
($\hat{\mathbf{p}}\approx\mathbf{p}$) does not guarantee high audio quality
($\hat{\mathbf{y}}\approx\mathbf{y}$). While our approach may not yield the
best parameter estimation, it does yield the best audio transform estimation.

##  4 Conclusion

In this work, we proposed an auto-encoder approach for Blind Estimation of
Audio Effects. Given only wet (processed) audio signals, we train a neural
network to estimate AFX parameters such that when used for effects applied to
a dry (unprocessed) signal it approximates the wet signal. This allows
training a network using real dry/wet data pairs without knowing the exact
effect implementation. We show that our audio-based method better replicates
the audio quality of the mastering process than the previous parameter-based
method.

Future works will focus on including other important mastering effects in the
chain and testing their estimation on real mastered music productions.

## References

  * [1] F. R. Ávila and L. W. P. Biscainho, “ML estimation of memoryless nonlinear distortions in audio signals,” in ICASSP 2014, pp. 4493–4497, IEEE.
  * [2] M. Ravanelli and Y. Bengio, “Speaker recognition from raw waveform with SincNet,” in SLT 2018, pp. 1021–1028.
  * [3] J. Engel, L. H. Hantrakul, C. Gu, and A. Roberts, “DDSP: Differentiable digital signal processing,” in ICLR 2020 Conference Track Proceedings, ICLR.
  * [4] S. Nercessian, “Neural parametric equalizer matching using differentiable biquads,” in Proceedings of the 23rd Int. Conf. on Digital Audio Effects (DAFx2020), pp. 265–272.
  * [5] S. Lee, H.-S. Choi, and K. Lee, “Differentiable artificial reverberation,” vol. 30, pp. 2541–2556.
  * [6] C. J. Steinmetz, N. J. Bryan, and J. D. Reiss, “Style transfer of audio effects with differentiable signal processing,” vol. 70, no. 9, pp. 708–721.
  * [7] C. J. Steinmetz, J. Pons, S. Pascual, and J. Serrà, “Automatic multitrack mixing with a differentiable mixing console of neural audio effects,” in ICASSP 2021, pp. 71–75, IEEE.
  * [8] M. A. Martínez Ramírez, O. Wang, P. Smaragdis, and N. J. Bryan, “Differentiable signal processing with black-box audio effects,” in ICASSP 2021, pp. 66–70, IEEE.
  * [9] J. T. Colonel and J. D. Reiss, “Reverse engineering of a recording mix with differentiable digital signal processing,” vol. 150, no. 1, pp. 608–619.
  * [10] R. Hinrichs, K. Gerkens, and J. Ostermann, “Classification of Guitar Effects and Extraction of Their Parameter Settings from Instrument Mixes Using Convolutional Neural Networks,” in Artificial Intelligence in Music, Sound, Art and Design (T. Martins, N. Rodríguez-Fernández, and S. M. Rebelo, eds.), vol. 13221, pp. 101–116, Springer International Publishing.
  * [11] S. Lee, J. Park, S. Paik, and K. Lee, “Blind estimation of audio processing graph,” in ICASSP 2023, pp. 1–5, IEEE.
  * [12] J. T. Colonel, M. Comunità, and J. D. Reiss, “Reverse Engineering Memoryless Distortion Effects with Differentiable Waveshapers,” in Audio Engineering Society Convention 153, p. 10, AES.
  * [13] R. Bristow-Johnson, “RBJ audio-EQ-cookbook.”
  * [14] D. Giannoulis, M. Massberg, and J. D. Reiss, “Digital dynamic range compressor design— a tutorial and analysis,” vol. 60, no. 6.
  * [15] C. J. Steinmetz and J. D. Reiss, “Efficient neural networks for real-time modeling of analog dynamic range compression,”
  * [16] E. Perez, F. Strub, H. De Vries, V. Dumoulin, and A. Courville, “FiLM: Visual Reasoning with a General Conditioning Layer,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 32, AAAI.
  * [17] A. Novak, L. Simon, P. Lotton, and J. Gilbert, “Chebyshev Model and Synchronized Swept Sine Method in Nonlinear Audio Effect Modeling,” in Proceedings of the 13th Int. Conf. on Digital Audio Effects.
  * [18] J. Koo, S. Paik, and K. Lee, “End-to-end music remastering system using self-supervised and adversarial training,” in ICASSP 2022, pp. 4608–4612, IEEE.
  * [19] J. Pons, O. Slizovskaia, R. Gong, E. Gómez, and X. Serra, “Timbre analysis of music audio signals with convolutional neural networks,” in 2017 25th European Signal Processing Conference (EUSIPCO), pp. 2744–2748.
  * [20] K. W. Cheuk, H. Anderson, K. Agres, and D. Herremans, “nnAudio: An on-the-fly GPU audio to spectrogram conversion toolbox using 1D convolutional neural networks,” vol. 8, pp. 161981–162003.
  * [21] J. Pons, “Deep neural networks for music and audio tagging.”
  * [22] Z. Rafii, A. Liutkus, F.-R. Stöter, S. I. Mimilakis, and R. Bittner, “MUSDB18 - a corpus for music separation.”
  * [23] C. J. Steinmetz and J. D. Reiss, “Auraloss: Audio-focused loss functions in PyTorch,” in DMRN+15.

  *[BE-AFX]: Blind Estimation of Audio Effects
  *[AFX]: Audio Effect
  *[AFXs]: Audio Effect
  *[NP]: Neural Proxy
  *[TCN]: Temporal Convolution Network
  *[FiLM]: Feature-Wise Linear Modulation