public double[] filter(double[] numerator, double[] denominator, double[] input) {
    assert numerator.length == denominator.length;

    int order = numerator.length - 1;
    double[] output = new double[input.length];

    output[0] = numerator[0] * input[0];
    for (int i = 1; i <= order; ++i) {
      output[i] = 0.0;
      for (int j = 0; j <= i; ++j) {
        output[i] += numerator[j] * input[i - j];
      }
      for (int j = 1; j <= i; ++j) {
        output[i] -= denominator[j] * output[i - j];
      }
    }

    for (int i = order + 1; i < input.length; ++i) {
      output[i] = 0.0;

      for (int j = 0; j <= order; ++j) {
        output[i] += numerator[j] * input[i - j];
      }
      for (int j = 1; j <= order; ++j) {
        output[i] -= denominator[j] * output[i - j];
      }
    }
    return output;
  }