Problem 438663 · hard · Level 04 Non-Linear Data Structures

The Guitar Body's Resonances, Read From Its Poles

poles · resonance · pole angle · decay time · complex conjugates · frequency

A luthier taps an acoustic guitar body and records the sound at fs samples per second. A model fitted to the recording is a recursive filter whose denominator a (powers of z^-1) describes how the body rings by itself. Its poles tell the luthier the body's resonances without any Fourier transform.

A pair of conjugate poles p = r * exp(+-1j * angle) (with 0 < angle < pi for the one above the real axis) contributes r ** k * cos(angle * k + phase) to the impulse response: a sine that turns by angle radians per sample, shrinking by the factor r every sample. So:

  • its frequency is angle * fs / (2 * pi) hertz (angle = cmath.phase(p)), in the same way as the turning arrows of Level 3;
  • its decay time, the time for the ringing to fall to 1 / e (about 37 %) of its size, follows from r ** n = exp(-1): n = -1 / log(r) samples, which is -1 / (fs * log(r)) seconds (natural logarithm, r = abs(p)). Poles close to the circle ring for a long time.

Real poles (on the real axis) give plain decays (or, for negative ones, sample-by-sample sign flips), not resonances.

Use the setup's roots_of(a) to find the poles (the roots of a read as powers of z, highest first). Write resonances(a, fs) that returns one pair (frequency_hz, decay_s) for every pole with imaginary part greater than 1e-6 (one from each conjugate pair), sorted by frequency. To see the tap ring, paste your run_transfer from the seismometer problem and press Run with plot(run_transfer([1], a, [1] + [0] * 400)).

Examples

Input:  a = [1, -0.9, 0.81], fs = 100
Output: [(16.666666666666668, 0.09491221581029916)]
Explanation: the poles are 0.9 * exp(+-1j * pi / 3): pi / 3 rad per sample is a sixth of a cycle,
so 100 / 6 = 16.7 Hz; -1 / (100 * log(0.9)) = 0.095 s.

Input:  a = [1, -1.4, 1.9, -1.301, 0.8064, -0.2592], fs = 100
Output: [(16.666666666666664, 0.09491221581029893), (25.0, 0.044814201177245536)]
Explanation: this is (1 - 0.9 z^-1 + 0.81 z^-2)(1 + 0.64 z^-2)(1 - 0.5 z^-1): the pair above,
a pair 0.8 * exp(+-1j * pi / 2) at a quarter of the sample rate, and a real pole at 0.5.

Input:  a = [1, -0.5], fs = 100
Output: []

Constraints

  • a[0] != 0, len(a) <= 11; every pole has abs < 1 (the body stops ringing), and no pole is repeated
  • the resonances' frequencies are distinct, and every complex pole has an imaginary part of at least 1e-3; roots_of returns real poles with an imaginary part far below 1e-6
  • answers are compared with a tolerance of 1e-6; a list may replace a tuple

Goals

  • Read a resonant frequency from the angle of a complex pole: f = angle * fs / (2 pi)
  • Read how long the resonance rings from the pole's distance to the origin
  • Ignore real poles, which decay without ringing
Starting Python…