Hacker News Reader

Comment by willis936 | Hacker News Reader

This is off topic but I’m assuming difficult-to-gather DSP peers will be in this comments section.

I have run into the limit of double precision floating point when making high order (20+) IIR filters in MATLAB. My solution in the past has been FIR filters and waiting a long time. I recently experimented with alternative forms of IIR representation. I found that converting from transfer function (B,A) to zero-pole-gain (z,p,k) form, then to second order sections (sos,g) form, then using the MATLAB filter function, the issues with double floating point precision go away in my use cases.

Is this a well known trick that I just stumbled upon?

View on HN

Yes, second-order/biquad sections are less sensitive to coefficient quantization than high-order filters. I think they benefit from interleaving zero+pole stages as well.

You can even ask MATLAB to generate the (z,p,k) form directly rather than (b,a) form when using its filter-design functions, saving a step (just ask for a three-argument result instead of a two-argument result).

(I am not an expert but I just so happened to have run across this feature and explanation last week, after vaguely remembering from college that biquad sections are "better" for some reason.)

This is a well known trick. See “Filter Realization” section here: https://en.m.wikipedia.org/wiki/Digital_filter

High-order polynomials are very sensitive to small deltas in the coefficients of their high-order terms. Representing the same function as the product of low-order factors is much more stable numerically.

I'm no DSP expert, but I did run into a problem in GNU Radio which I think applies because it happens to floating point numbers in general.

GNU radio has a component called "DC Block" which has a length of N samples... it averages those samples, and subtracts that average from the input... giving an output.

If you exceed the dynamic range of floating point by mixing inputs of lets say 1e12 and 1e-12... the filter can develop a permanent output bias.

The reason is performance... the samples are stored in a ring buffer, and they subtract the value when it leaves the buffer, and add the new value as it goes it. This is N times faster than adding up all the values, and subtracting. However... this bug can happen.

I suspect you have similar behavior occurring due to the smaller coefficients at the tails of the function.

That's my 2cents, I hope it made sense, and helps explain what you've encountered.

Yeah it’s definitely a FP dynamic range issue. Moving to extended point is impractical for a great many reasons. Transfer function form IIRs put high demands on FP dynamic range that other forms do not. This has real world parallels too, since the most basic way to realize an analog filter is with the transfer function method. I don’t know of anyone who makes a 2 F capacitor with tolerances that are within 2 pF.

Matlab has a helper function for this exactly, tf2sos. It's also covered in most applied DSP texts.

Reply on news.ycombinator.com