Problem 590957 · medium · Level 05 Advanced Algorithms & Graphs

A Minute of Vibration from the Engine Room

numpy.fft · Hann window · amplitude spectrum · peak picking · vibration analysis

A ship's engineer clips a vibration logger to the main engine's mounting and records a whole minute or more at fs samples per second: well over a hundred thousand samples. Each rotating part (the crankshaft, the turbocharger, a pump) shakes the hull at its own frequency, and a new peak in the spectrum is an early sign of a fault. The engineer wants the list of peaks.

A recording this long needs the FFT, and here you may use numpy: numpy.fft.rfft(v) returns the bins X[0], ..., X[N // 2] of the DFT X[k] = sum over n of v[n] * exp(-2j * pi * k * n / N) in a fraction of a second, for any length N. The steps:

  1. Multiply the recording by the Hann window w[n] = 0.5 - 0.5 * cos(2 * pi * n / N) for n = 0 .. N-1 (this reduces the leakage of a wave that does not fit a whole number of times into the recording). Use exactly this formula: numpy.hanning uses N - 1 in place of N and gives slightly different numbers.
  2. Take X = numpy.fft.rfft(v * w) and the amplitude of each bin, A[k] = 2 * |X[k]| / sum(w). With this scaling a wave a * cos(...) that falls exactly on a bin reads A[k] = a.
  3. A bin k with 1 <= k <= N // 2 - 1 is a peak when A[k] > A[k-1], A[k] >= A[k+1] and A[k] >= threshold.

Write vibration_peaks(x, fs, threshold) that returns the list of (frequency, amplitude) pairs, one per peak, in order of frequency: the frequency is k * fs / N hertz and the amplitude A[k].

The tests use the helper engine_log(seconds, fs, seed), which returns a made-up recording of seconds * fs samples: a few waves of random frequency, an offset and some noise. You can call it with Run, for example x = engine_log(60, 2000, 1), and press plot(freqs, amps) with your arrays to see the spectrum.

Examples

Input:  x = [1, 0, -1, 0] * 4, fs = 4, threshold = 0.1
Output: [(1.0, 1.0)]
Explanation: a cosine of amplitude 1 at 1 Hz (a quarter of the sampling rate, bin 4 of 16).

Input:  x = [2 + math.cos(2 * math.pi * 3 * n / 32) + 0.5 * math.cos(2 * math.pi * 10 * n / 32 + 1) for n in range(32)],
        fs = 32, threshold = 0.1
Output: [(3.0, 1.0), (10.0, 0.5)]
Explanation: two waves exactly on bins 3 and 10; the offset 2 lives in bin 0, which is never a peak.

Constraints

  • 16 <= N <= 300000; answers are compared with a tolerance of 1e-6
  • no amplitude in the tests is within a few per cent of threshold

Goals

  • Use numpy.fft.rfft on a recording far too long for a hand-written transform
  • Window the recording and scale the bins so a peak reads as the wave's amplitude
  • Report the local peaks of the spectrum that stand above a threshold
Starting Python…