Problem 586664 · hard · Level 05 Advanced Algorithms & Graphs

Tuning the Step Counter's Filter in Decibels

noise reduction · signal-to-noise ratio · decibels · FIR filter design · numpy.convolve · group delay

A phone's step counter reads the accelerometer fs times per second. Walking makes a smooth wave (the step frequency, around 2 Hz, and a few multiples of it), and the sensor adds random noise across all frequencies. A low-pass filter removes noise, but a cutoff that is too low also flattens the walking wave itself. On a test bench the true wave clean is known exactly alongside the sensor's reading noisy, so the designer can measure which cutoff in the list cutoffs works best.

For each cutoff fc:

  1. Design the Hamming windowed-sinc low-pass with taps weights (odd): with M = taps - 1 and t = k - M / 2 for k = 0 .. taps - 1, ideal[k] = 2 * fc / fs if t == 0, else sin(2 * pi * fc * t / fs) / (pi * t); multiply by w[k] = 0.54 - 0.46 * cos(2 * pi * k / M); divide by the sum so the weights add up to 1.
  2. Filter and align: compute the full convolution y[n] = sum over k of h[k] * noisy[n - k] (zero outside the recording) and keep y[D], ..., y[D + N - 1], where D = (taps - 1) / 2 is the filter's delay and N = len(noisy). The aligned output lines up with clean.
  3. Measure only on the samples i = D, ..., N - 1 - D, where the filter saw real data on both sides. The noise power before is the mean of (noisy[i] - clean[i]) ** 2, after it is the mean of (aligned[i] - clean[i]) ** 2 (it includes the damage done to the wave). The noise reduction is 10 * log10(before / after) decibels: positive means the filtered signal is closer to the truth.

Write best_cutoff(clean, noisy, fs, taps, cutoffs) that returns the pair (best, reductions): the list of noise reductions in decibels, one per cutoff in order, and the cutoff with the largest one (the first such, on a tie).

The recordings are long (a minute at 200 samples per second, filtered with 101 weights, for eight cutoffs: about 10 million multiplications), so this problem lets you use numpy: numpy.convolve(a, b) returns the full convolution. The tests use the helper walk_bench(seconds, fs, seed), which returns a made-up pair (clean, noisy) of seconds * fs samples; you can call it with Run. Press Run with plot(cutoffs, reductions, kind="scatter") to see the trade-off.

Examples

Input:  clean = [1, 1, 1, 1, 1], noisy = [1, 3, -1, 3, 1], fs = 10, taps = 3, cutoffs = [2.5]
Output: (2.5, [1.4504481522878987])
Explanation: the weights are [0.0462, 0.9076, 0.0462] (the outer sinc values 0.3183 times the
window 0.08, the middle 0.5, then normalised). Measured on samples 1 to 3, the squared error
falls from mean([4, 4, 4]) = 4 to mean([2.97, 2.66, 2.97]) = 2.864: 10 * log10(4 / 2.864) = 1.45 dB.

Input:  (clean, noisy) = walk_bench(60, 200, 1), fs = 200, taps = 101,
        cutoffs = [1, 2, 4, 6, 8, 12, 20, 40]
Output: (6, [-0.17226417825631987, 1.3354285254445928, 7.848537373210894, 12.22982467227254,
             11.19780321653176, 9.154329549259847, 6.919428527347156, 3.9721369809866047])
Explanation: below about 6 Hz the filter cuts into the walking wave (at 1 Hz it is worse than
no filter at all); above it, more noise gets through.

Constraints

  • taps odd, 3 <= taps <= 301, taps <= N <= 10**5; 0 < fc < fs / 2 for every cutoff; 1 <= len(cutoffs) <= 20
  • noisy differs from clean somewhere in the measured range; answers are compared with a tolerance of 1e-6

Goals

  • Measure a filter's noise reduction in decibels against a known clean signal
  • Line up the filtered signal by the group delay before comparing it
  • Find the cutoff that trades removed noise against a distorted signal
Starting Python…