IIR Bandpass Filter - Numerical Precision Causing Instability?
14 views (last 30 days)
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?
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
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.