IIR Bandpass Filter - Numerical Precision Causing Instability?

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

Answer 1

[Star Strider]

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.

[Tom Andersson]

I am using the designfilt() command with user-chosen frequency and attenuation specifications to design the digital filter. Then I am filtering by passing the digital filter object to filtfilt(). I'm a bit of a MATLAB & DSP newbie so call me out if this is a filter design sin. I've replaced my designfilt() approach with yours and it seems to work perfectly, so thank you very much!

I'm still really curious as to why my 'stable' iirbandpass from designfilt() turned out to be unstable, and why the 'identical' process of cascading an iirlowpass and iirhighpass did not.

[Star Strider]

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.