lfilter function numerical stability

🐛 Bug

When using lfilter function in torchaudio/functional.py, I input two groups of biquad coefficient parameters (one group is normalized by a0, and the other is unnormalized). The maximum error between two output waveforms is 6e-4.

To Reproduce

Steps to reproduce the behavior:

I have a group of coefficients calculated from the bass biquad filter:

[b0, b1, b2, a0, a1, a2] = 
[27.26210639724472 -37.770844011587144 14.561748502365816 20.54382514976342 -39.77708440115872 19.27378936027553]

After normalization by a0:

[b0/a0, b1/a0, b2/a0, a0/a0, a1/a0, a2/a0] = 
[1.3270219250069244 -1.8385497216920241 0.7088138842796521 1.0 -1.9362063350513272 0.9381791959272727]

The details of coefficient calculation are shown in Additional Test (list below). Then I compare the output waveforms and the number of elements that has error (absolute difference of output values) larger than 1e-4 is 12603. The maximum error between two waveforms is 6e-4.

filename = 'whitenoise.wav'
noise_filepath = common_utils.get_asset_path(filename)
waveform, sample_rate = torchaudio.load(noise_filepath, normalization=True)
device = waveform.device
dtype = waveform.dtype

output_waveform1 = F.lfilter(
waveform,
torch.tensor([a0, a1, a2], dtype=dtype, device=device),
torch.tensor([b0, b1, b2], dtype=dtype, device=device),
)

output_waveform2 = F.lfilter(
waveform,
torch.tensor([a0/a0, a1/a0, a2/a0], dtype=dtype, device=device),
torch.tensor([b0/a0, b1/a0, b2/a0], dtype=dtype, device=device),
)

atol = 1e-4
diff = abs(output_waveform2 - output_waveform1)
print((diff > atol).sum())
v, i = diff.topk(1)
print(v)

output result:
12603
6e-4

Expected behavior

There is negligible error between two output waveforms.

Environment

• PyTorch Version: 1.6.0a0+78acc9d
• OS: Linux
• How you installed PyTorch: Conda
• Build command you used : As described above
• Python version: 3.7
• CUDA/cuDNN version: CUDA 10.1

Additional context

The maximum error deceases when gain decreases.
The method to calculate the coefficients of bass (biquad) filter:

def bass_biquad_coefficient(
        #waveform: Tensor,
        sample_rate: int,
        gain: float,
        central_freq: float = 100,
        Q: float = 1
) -> Tensor:
    r"""Design a bass tone-control effect.  Similar to SoX implementation.
    Args:
        waveform (Tensor): audio waveform of dimension of `(..., time)`
        sample_rate (int): sampling rate of the waveform, e.g. 44100 (Hz)
        gain (float): desired gain at the boost (or attenuation) in dB.
        central_freq (float, optional): central frequency (in Hz). (Default: ``100``)
        Q (float, optional): https://en.wikipedia.org/wiki/Q_factor (Default: ``0.707``).
    Returns:
        Tensor: Waveform of dimension of `(..., time)`
    References:
        http://sox.sourceforge.net/sox.html
        https://www.w3.org/2011/audio/audio-eq-cookbook.html#APF
    """
    w0 = 2 * math.pi * central_freq / sample_rate
    alpha = math.sin(w0) / 2 / Q
    A = math.exp(gain / 40 * math.log(10))

    temp1 = 2 * math.sqrt(A) * alpha
    temp2 = (A - 1) * math.cos(w0)
    temp3 = (A + 1) * math.cos(w0)

    b0 = A * ((A + 1) - temp2 + temp1)
    b1 = 2 * A * ((A - 1) - temp3)
    b2 = A * ((A + 1) - temp2 - temp1)
    a0 = (A + 1) + temp2 + temp1
    a1 = -2 * ((A - 1) + temp3)
    a2 = (A + 1) + temp2 - temp1
    
    return [b0, b1, b2, a0, a1, a2]
    
central_freq = 1000
s = 1
gain = 40
sample_rate = 44100
A = math.exp(gain / 40 * math.log(10))
Q = 1/math.sqrt(2 + (A + 1/A) * (1/s - 1)) 

[b0, b1, b2, a0, a1, a2] = bass_biquad_coefficient(sample_rate, gain, central_freq, Q)

Start in torchaudio/functional.py and reproduce the issue with the two coefficient groups and whitenoise.wav shown in the report. Compare the two output waveforms using the stated absolute-error checks; done means normalization produces negligible differences without regressing lfilter behavior.