Time-varying Notch Filter in C

Question

I have a digital notch filter (band-stop) implemented in C. It works very well. The issue is that varying the filter coefficients in order to track a noise signal that varies in frequency results in a transient response that can be more noisy than the incoming signal. We refer to this as "shot" noise and it appears this is a well-known problem with time-varying notch filters. I have looked at a couple of papers that appear to damp the changes to the coefficients in order to avoid this effect, but the maths is beyond me and I am really just looking for a C example that does something like this.

EDIT: this problem appears worse in our case because we have cascading notches - at least 24 all being updated at the same time - it seems that "shot" noise gets amplified in this case if you get the wrong combination of inputs.

EDIT2: An additional problem we are seeing is that the filter cascade becomes unstable at certain points - so we get much larger output than is possible on the input. This only happens when we are linking several notch filters together. This is all in floating point. Is there a way of avoiding this instability?

Answer 1

A few ideas:

If you are implementing this as an IIR filter, use Direct Form I. This minimizes the discontinuity of the state variables (which are simply input and output).

A biquad notch filter can be implemented as the sum of unity and a second order allpass, i.e.

$$H(z) = \frac{1}{2}\left(1+ \frac{a_2+a_1z^{-1} + z^{-2}}{1+a_1z^{-1} + a_2 z^{-2}} \right) $$

In this structure only two coefficients need to be updated.

Don't update too fast. Big jumps in frequency will create big jumps in the frequency response. To slow things down you can put a first order lowpass on the notch frequency and adjust the corner frequency to dial in the trade off between responsiveness and artifacts.

implemented in C.

I would first develop, verify and test the algorithm in something that can easily be instrumented like Matlab, Octave or Python. Once you have a stable fully unit tested reference, you can port to C .

Answer 2

After quite a lot of analysis we concluded that the main culprit causing the instability is the rate of change of the center frequency driving the coefficients. We found that a simple slew of:

$$ f'= \begin{cases} f(1-s) & \text{if $f_\lambda < f(1-s)$}\\ f\frac{1}{(1-s)} & \text{if $f_\lambda > f\frac{1}{(1-s)}$}\\ f_\lambda & \text{otherwise} \end{cases} $$

where $f_\lambda$ is the requested frequency, $f$ is the current frequency and $f'$ is the new frequency, worked quite well for reasonable update rates.

Answer 3

In the most simple cases in discrete time, it is enough to assume your noise is white, for designing a linear filter.

When the noise is not white, such as what you describe, an additional term could be useful.

Linear parameter estimation, such as the ARMAX type models, are based on a specific noise probabilistic distribution:

$$ A(z)y(t)=B(z)u(t)+C(z)e(t) $$

Without more information about these spots, it is impossible to say anything else. But with data, you can estimate $C(z)$ and eliminate those noise effects straightforwardly without trouble.

Most discrete time filter design packages usually integrate the 'Moving Average' component, which you can use specifically in this case for this effect. But if that is not the case, or it results too cumbersome to obtain from your package, perhaps you must guide the design manually to address this, or even, to calculate manually a estimate of the noise response $F(z)$, and including that estimate explicitly as a first stage compensator $C(z)=1/F(z)$ in a proper smooth way for this compensator to work, to make the noise $C(z)e(t)$ enough white for the subsequent filter to operate well.