# IIR Bandpass Filter - Numerical Precision Causing Instability?

15 views (last 30 days)
Tom Andersson on 8 Aug 2017
Commented: Star Strider on 8 Aug 2017
I am using a Chebychev Type 2 Bandpass IIR filter to filter some audio data. I am using the filtfilt() command to compensate for the phase distortion introduced by the IIR filter. My low frequency cut-off of the bandpass is ~1% of the Fs/2 frequency, giving a number of poles very close to the 1+0j point in the z-plane, but none outside of the unit circle and hence isstable(bandpassFilter) returns true. See the pole and zero locations below:
However, when I implement the filter with y = filtfilt(bandpassFilter,y) I get a crazy output:
I suspect that this is due to numerical precision issues resulting from the filter coefficients being stored as floating-point values, causing some of the ~1+0j poles to move outside of the unit circle, or perhaps some kind of ill-conditioning?
Strangely, if I exactly mimic the bandpass filter in the frequency domain by applying a highpass followed by a lowpass Cheby2 IIR with filtfilt(), I get the stable output that I expect. I also get a stable output if instead of using filtfilt() I use filter() followed by flipud() twice. Hence, could the issue lie in the way filtfilt() handles data?
Cheers

Star Strider on 8 Aug 2017
The problem could be with the way you are designing the filter. Some functions are better than others (I prefer the â€˜z,p,kâ€™ rather than the transfer function form for the filter), and all need to be converted to second-order-section form for stability.
Try something like this:
Fs = 44100; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency (Hz)
Wp = [5000 10000]/Fn; % Passband Frequency (Normalised)
Ws = [4990 10010]/Fn; % Stopband Frequency (Normalised)
Rp = 10; % Passband Ripple (dB)
Rs = 75; % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs); % Filter Order
[z,p,k] = cheby2(n,Rs,Ws); % Filter Design
[sosbp,gbp] = zp2sos(z,p,k); % Convert To Second-Order-Section For Stability
figure(1)
freqz(sosbp, 2^16, Fs) % Filter Bode Plot
filtered_signal = filtfilt(sosbp, gbp, signal); % Filter Signal
I always use the freqz function to analyze the filter passband to be certain I am getting the result I want, and that the filter is stable.
##### 2 CommentsShow 1 older commentHide 1 older comment
Star Strider on 8 Aug 2017
As always, my pleasure!
I always design my own filters because I can control the filter characteristics with greater precision. I sometimes refer people to designfilt for simple filters (such as the 60 Hz notch filter in the documentation on Remove the 60 Hz Hum from a Signal (link)), since it is easier for those who need a simple filter design to use. I post code similar to what I posted here in my Answers, since it is preferable to going into a long discussion on using designfilt.
I haven’t delved into designfilt, so I don’t know what design decisions it makes, although I suspect it uses transfer function representation. I am not certain if it outputs transfer function or second-order-section representation.
Your cascaded approach worked because the individual filters were apparently stable (easier with simpler filters), so cascading them produced a stable result. Bandpass designs, especially with sharp cutoffs, tend to be more difficult to design using transfer function representation rather than zero-pole-gain representation. I suspect the reason is that zero-pole-gain ‘translate’ more easily to second-order-section representation than do transfer function representations.