MATLAB: IIR Bandpass Filter – Numerical Precision Causing Instability

bandpassdigital signal processingfilterfilteringfiltfilthighpassiiriir filteringlowpassnumerical precisionpoles

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

Best Answer

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.
Related Question