ObieWebApp

1.1 What is a FRF?

We want to express the violin mathematically as a transfer function (also called the frequency response function) — it converts some input (in this case a force from the hammer) to some output (in this case sound or a pressure fluctuation):

Taking the Fourier Transform, we can find the transfer function for the violin as:

where \(P(f)\) and \(F(f)\) are the Fourier transforms of pressure and force and have real and imaginary parts. The Fourier Transform is defined as:

We can simplify — get the denominator to be real — by multiplying numerator and denominator by \(F^*\):

1.2 Least Squares Fit — H1 and H2

In reality, we take 5–10 hits per position so that we can get a "best fit" reading. We know that:

So the error of the best fit for measured values \(F_i\) and \(P_i\) is:

Realizing that \(H(f)\) is a complex number \((a+ib)\) and that \(F, P\) are complex as well:

Summing over all readings and finding the minimum by differentiating with respect to \(a\) and \(b\) and setting equal to zero:

Multiplying the second equation by \(i\) and then adding:

Adding them together:

This gives the H1 estimator (minimises output noise):

Alternatively, we can minimise the effect of noise on the input by using as the error:

Expanding the squared error:

Differentiating and setting to zero:

Multiplying the second by \(i\) and adding:

Adding and rearranging:

This gives the H2 estimator (minimises input noise):

If the data is clean and the noise is truly random, H1 and H2 should be very close and their ratio — the coherence — should be 1. George used the geometric mean as the working estimate:

1.3 Coherence

The coherence \(C\) measures how consistently the system behaves across repeated measurements:

The numerator is the cross-spectrum squared — how much \(F\) and \(P\) move together. The denominator is the product of the individual power spectra — how much each signal varies on its own. It is the frequency-domain analogue of \(R^2\) from linear regression, and answers: "how much of the output can be explained by the input?"

• \(C \approx 1\): output is well explained by input — clean measurement, linear system.
• \(C \approx 0\): output has little to do with input at that frequency — noise dominated, or the system simply doesn't respond there.
• \(0 < C < 1\): partial explanation — could be noise, nonlinearity, or multiple inputs contributing to the output.

We can look at the average distribution of energy by breaking the FRF into frequency bands and averaging the magnitudes over each band:

We also calculate the spectral centroid of the band — the first moment — which tells you where the energy is concentrated within the band:

In Python (Python/processing/bands.py):

f_lo  = float(band['start'])
f_hi  = float(band['end'])
mask     = (freq >= f_lo) & (freq <= f_hi)
H_band   = H_mag[mask]          # linear magnitude, NOT in dB
f_band   = freq[mask]
avg_db   = 20.0 * math.log10(max(float(np.mean(H_band)), 1e-12))
centroid = float(np.sum(f_band * H_band) / np.sum(H_band))

We can estimate what an instrument will sound like by convolving its FRF with a recording — ideally captured with a force sensor at the bridge rather than a microphone. The convolution of two functions \(f(t)\) and \(g(t)\) is:

What it does

A violin has resonances — certain frequencies ring longer and louder than others. The FRF is a fingerprint of exactly how the violin shapes sound. Convolution applies that fingerprint to any audio signal. Feed in a Tchaikovsky recording; get out what it would sound like coloured by that violin's resonance structure.

FFT acceleration

Computing the integral directly is slow. The convolution theorem says we can instead multiply the two signals in the frequency domain:

In Python (Python/processing/convolution.py):

y = np.fft.irfft(
    np.fft.rfft(wav, n_fft) * np.fft.rfft(ir.astype(np.float64), n_fft),
    n_fft,
)[:N].real

The output is the input signal with the violin's character stamped onto it — its resonant peaks amplified, its dead spots suppressed, its natural decay times imposed on every note.

IR length

The impulse response (IR) length determines frequency resolution. ObieApp uses an adaptive length: the next power of 2 above \(\text{sample\_rate} / \Delta f\), capped at 65536:

def _adaptive_ir_length(freqs, sample_rate):
    df = freqs[1] - freqs[0]
    natural = int(sample_rate / df)
    next_pow2 = 1 << (natural - 1).bit_length()
    return min(next_pow2, 65536)

Minimum-phase reconstruction

For real-only FRF files (AvR format, no phase information), the phase can be estimated by treating the system as minimum-phase — for such systems the log-magnitude and phase are a Hilbert transform pair, so knowing one determines the other. In practice ObieApp stores all FRFs as complex (AvC or TRF), so minimum-phase reconstruction is not required.

All multi-byte fields are little-endian. Fields prefixed f in the MtxVec TSingleArray descriptor are 32-bit single-precision floats; all others are as specified. Sources: george.pas (Dew Research / DSP Master) and the LabVIEW writer code.

4.1 TRF / TRV — Transfer Function Format

Binary format used by Rational Acoustics / MtxVec. The header is 110 bytes, followed immediately by the data payload.

TFRFHeader (110 bytes, packed — no padding)

4.2 AvR / AvC — Averaged Spectrum Format

AvR files store a real-valued averaged spectrum; AvC files store a complex-valued averaged spectrum. The file layout is identical except for the data section, which is \(N \times 8\) bytes for AvR and \(N \times 16\) bytes for AvC. \(N\) is given by fLength in the TSingleArray descriptor.

This section derives the linearised acoustic wave equation from first principles and applies it to the mechanics of violin plates and the sound hole.

5.1 Wave Equation Derivation

Assume that the pressure is \(p\), the speed of a fluid particle is \(u\) with a mean of zero, and the density everywhere is \(\rho_0 + \rho\). The 1D continuity equation gives:

Expanding and dropping the nonlinear term \(\partial(\rho u)/\partial x\):

Differentiating again with respect to \(t\):

The momentum equation (1D, no gravity or friction) linearised for small fluctuations:

Differentiating with respect to \(x\) and dropping the middle (nonlinear) term:

Combining the continuity and momentum equations:

Adiabatic closure

For an adiabatic process \(p = K\rho^\gamma\), so:

Linearising the second time derivative:

Substituting into the combined equation gives the acoustic wave equation:

where the speed of sound \(c = \sqrt{\gamma RT}\).

5.2 Acoustic Impedance

Specific acoustic impedance \(z\) is the ratio of pressure to acoustic particle velocity: \(p = zu\), analogous to voltage driving a current.

For a plate with mass \(m\), stiffness \(k\), and damping \(\beta\), area \(A\), the equation of motion is:

Solving for impedance:

Assuming a harmonic solution \(u = \dot{x} = Ce^{i\omega t}\):

where the lumped parameters are:

This is the direct mechanical analogue of an electrical RLC circuit: resistance \(r\), inductance (mass) \(M\), and compliance (capacitance) \(C\).

For the sound hole, modelled with frictional resistance \(R_h\), and the top plate with restoring force:

5.3 Sound Intensity

Sound intensity is defined as the sound power per unit surface area. For a point source of power \(P\):

Sound Power Level and Sound Intensity Level are referenced to \(P_0 = 10^{-12}\,\mathrm{W}\) and \(I_0 = 10^{-12}\,\mathrm{W/m^2}\) respectively:

Since \(I = pu = p^2/z\), the Sound Pressure Level is:

referenced to \(p_0 = 20\,\mu\mathrm{Pa}\).

5.4 Quality Factor

The quality factor \(Q\) characterises the sharpness of a resonance peak. For the damped oscillator \(m\ddot{x} + R\dot{x} + kx = 0\):

where \(\Delta f\) is the half-power bandwidth (the −3 dB width of the peak), \(r\) is the exponential decay rate, and \(R\) is the damping coefficient. Higher \(Q\) means a sharper, more persistent resonance.

Modal Analysis (Web/tools/modeshape) animates how a plate physically deforms at a chosen frequency, using nothing more than each stencil node's measured complex FRF at that frequency. Because every node's hammer-impact FRF is referenced to the same fixed input, comparing their complex values at a single frequency reveals both relative amplitude and relative phase across the surface — which is exactly what a mode shape is.

6.1 The Complex FRF as a Phasor

At frequency \(f\), node \(i\)'s measured FRF is a complex number \(H_i(f) = \mathrm{Re}_i + i\,\mathrm{Im}_i\). Treat this as a phasor: under sinusoidal excitation at that frequency, the steady-state physical displacement of node \(i\) over time is the real part of \(H_i e^{i\theta}\), where \(\theta=\omega t\) is the drive phase:

Sweeping \(\theta\) from \(0\) to \(2\pi\) animates one full vibration cycle. Points where \(H_i \approx 0\) barely move (nodal lines); points with large \(|H_i|\) move the most (antinodes); points whose \(H_i\) have opposite sign move in opposite directions at the same instant, which is how the animation shows a plate's characteristic node lines.

6.2 Building a Continuous Surface

Stencil nodes sit on a rectangular row/column grid at known physical (x, y) positions. For the currently selected frequency, Modal Analysis:

1. Linearly interpolates each node's own FRF between its two nearest measured frequency bins to get \(H_i(f)\) at the exact requested frequency.
2. Arranges the per-node Re and Im values into two 2-D grids indexed by row and column.
3. Fits a tensor-product natural cubic spline to each grid — first across columns per row, then across rows per fine column — producing smooth, \(C^2\)-continuous \(\mathrm{Re}(x,y)\) and \(\mathrm{Im}(x,y)\) surfaces. This is equivalent to scipy.interpolate.RectBivariateSpline, reimplemented in JS so it can run synchronously inside the animation loop.
4. Normalises both surfaces by the peak magnitude \(\max\sqrt{\mathrm{Re}^2+\mathrm{Im}^2}\) over the fine grid, so the animated peak amplitude is always 1 regardless of instrument or frequency — mode shapes stay visually comparable across runs.

Each animation frame then evaluates the phasor formula per fine-grid point:

// modeshape.js — one animation frame
disp[y][x] = amp * (Re_fine[y][x]*cos(θ) - Im_fine[y][x]*sin(θ))

Circle Fit (Web/tools/circlefit, backed by Python/processing/circlefit.py) extracts a natural frequency, damping, and modal constant for each resonance directly from a measured accelerance FRF, using the classical Kennedy–Pancu circle-fitting method (Ewins, Modal Testing: Theory and Practice).

7.1 Receptance vs. Accelerance

For a single-degree-of-freedom (SDOF) resonance with hysteretic (structural) damping, it is receptance — displacement per unit force — that traces an exact circle in the Nyquist (Re vs Im) plane:

where \(A_r\) is the complex modal constant (residue), \(\omega_r\) the natural frequency (rad/s), and \(\eta_r\) the hysteretic loss factor. Acquire's hammer-impact measurements are accelerance, \(A(\omega) = -\omega^2 X(\omega)\), so Circle Fit first converts:

7.2 Algebraic Circle Fit

Within a user-picked candidate band, the receptance points \((x_i,y_i)=(\mathrm{Re}_i,\mathrm{Im}_i)\) are fit to a circle \(x^2+y^2+Dx+Ey+F=0\) by ordinary least squares (Kåsa method) directly on the design matrix \([x,\,y,\,1]\):

giving the centre and radius:

Fit quality is reported as the RMS radial residual as a fraction of the radius (resid_rms_frac) — points that don't actually lie on a circle (no real resonance in the band, or a badly chosen band) show up as a large value here.

7.3 Natural Frequency — Fastest Angular Sweep

Physically, the phasor sweeps around the Nyquist circle fastest exactly at resonance. The angle of each point relative to the fitted centre, unwrapped to avoid a branch-cut jump, is:

The natural frequency \(\omega_r\) is taken as the point of maximum \(\left|\dfrac{d\theta}{d\omega}\right|\) (central difference across interior points), with a 3-point parabolic interpolation around the discrete peak for sub-bin precision.

7.4 Damping — Half-Angle Relation

For every pair of points \(a\) (below resonance) and \(b\) (above resonance), the half-angles \(\theta_a,\theta_b\) subtended at the circle centre between each point and the resonance point give an independent damping estimate:

Pairs where either half-angle is below a 5° minimum are discarded (\(\tan(\theta/2)\to0\) is numerically unstable right at resonance), and the final \(\eta_r\) is the median across all surviving pairs — robust to the occasional noisy outlier pair in real hammer data.

7.5 Modal Constant and Sign

The residue magnitude follows directly from the circle diameter:

Its sign is read off the circle centre's vertical position rather than assumed: because the receptance denominator \(\omega_r^2-\omega^2+i\eta_r\omega_r^2\) traces a horizontal line at constant positive imaginary part as \(\omega\) varies, and inverting a line not through the origin (\(z\to1/z\)) always produces a circle on the opposite side of the real axis — a positive real \(A_r\) always puts the fitted circle's centre at \(y_0<0\), and a negative \(A_r\) at \(y_0>0\).

7.6 Residual Compensation and Multi-Pass Fitting

A real structure has more than one mode, and modes outside the current band still contribute a slowly-varying, near-real background to the data inside it. Each already-identified mode \(s\)'s approximate contribution at \(\omega\) is:

which is subtracted from the real part of the band's data before each circle fit. Because every mode's residual depends on every other mode's fit, Circle Fit iterates — up to 5 passes — refitting every band using the latest fits of all the others, stopping early once the largest relative frequency change across modes falls below a tolerance (default \(10^{-4}\)).

7.7 Per-Node Residue Extraction

Once a mode's frequency \(\omega_r\) and damping \(\eta_r\) are fixed from the reference (averaged) FRF, its contribution to any individual node's receptance is linear in that node's own unknown complex residue \(A_r\). This lets every node's residue — and hence a damping-corrected, per-mode mode shape — be solved directly from the whole spectrum in one complex least-squares system, rather than band by band:

where \(C\) and \(M/\omega^2\) absorb the near-constant and low-frequency-mass-like contributions of modes outside the fitted set.

Spectrum Monitor's averaged spectrum is a simpler, live-only cousin of the H1 estimator in §1 — there's no reference "input" channel to divide out, just one microphone channel's own magnitude spectrum, smoothed over the last few snippets so a steady tone or resonance holds still instead of jittering frame to frame.

8.1 Windowed FFT per Snippet

Every snippet (default 0.3 s, user-adjustable) is Hann-windowed and transformed independently:

Dividing by half the window's sum before squaring keeps the resulting power roughly amplitude-correct regardless of window shape or FFT length — the same normalisation convention used for a single-sided amplitude spectrum. This is implemented directly in Web/tools/spectrum/main.py (_process_snippet) using numpy.fft.rfft — a plain FFT, not a reuse of any canonical Python/processing/ module, since no existing module computes a single channel's own magnitude spectrum this way.

8.2 Power Averaging, Not dB Averaging

The displayed spectrum block-averages the last \(N\) snippets' power, then converts to dB once at the end:

Changing the channel resets this running average immediately (mixing two channels' power together mid-average would be meaningless); changing \(N\) or the snippet length does not force a reset — the averaging window just grows, shrinks, or (if the FFT size itself changes) silently restarts once the accumulated spectra no longer match the new bin count.

The Spectrogram tab is unrelated math — it reuses Python/processing/spectrogram.py's compute_spectrogram STFT unchanged, the same function Convolve uses for its input/output spectrogram panes.