A seismometer's datasheet specifies its built-in filter only by its transfer function H(z) = B(z) / A(z), as coefficient lists b and a in powers of z^-1, and the supplier did not bother to make a[0] equal to 1. To process recorded ground motion x you must run the filter yourself.
Y(z) * A(z) = X(z) * B(z), read with z^-i as a delay of i samples, is the difference equation
a[0] * y[k] + a[1] * y[k-1] + ... + a[N] * y[k-N] = b[0] * x[k] + b[1] * x[k-1] + ... + b[M] * x[k-M]
Solve it for the newest output:
y[k] = ( b[0] * x[k] + ... + b[M] * x[k-M] - a[1] * y[k-1] - ... - a[N] * y[k-N] ) / a[0]
The filter starts at rest: every x and y before sample 0 is 0.
Write run_transfer(b, a, x) that returns the outputs y[0], ..., y[len(x) - 1] as a list. Press Run with plot(run_transfer(b, a, [1] + [0] * 30), kind="stem") to see the impulse response.
Examples
Input: b = [2, 2], a = [4, -2], x = [1, 0, 0, 0]
Output: [0.5, 0.75, 0.375, 0.1875]
Explanation: y[0] = (2 * 1) / 4 = 0.5; y[1] = (2 * 0 + 2 * 1 + 2 * 0.5) / 4 = 0.75;
y[2] = (0 + 0 + 2 * 0.75) / 4 = 0.375. (The -a[1] = +2 feeds the output back.)
Input: b = [0.1], a = [1, -0.9], x = [1, 1, 1]
Output: [0.1, 0.19, 0.271]
Constraints
a[0] != 0;0 <= len(x) <= 10**4- answers are compared with a tolerance of
1e-6
Goals
- Turn H(z) = B(z) / A(z) back into a difference equation you can run
- Divide by a[0] when the datasheet does not normalise it to 1
- Run a recursive filter of any order from rest