# Why this bandpass butterworth is unstable (while the corresponding low and high pass are stable)?

63 views (last 30 days)
Luca Amerio on 15 May 2019
Commented: Greg Dionne on 17 May 2019
I have a signal sampled at 750Hz
T = 60;
fs = 750;
t = (0:T*fs-1)/fs;
x = sin(6*t) + 0.1*rand(size(t));
I want to filter it between 50 and 80 Hz. If I use an high pass filter followed by a low pass everything is fine.
fHP = 50;
fLP = 80;
order = 3;
[b,a] = butter(order,fHP/(fs/2),'high');
x_hp = filter(b,a,x);
[b,a] = butter(order,fLP/(fs/2),'low');
x_lp = filter(b,a,x);
x_hp_lp = filter(b,a,x_hp);
figure
subplot(2,1,1)
hold on
plot(x ,'DisplayName','original data')
plot(x_hp ,'DisplayName','High pass')
plot(x_lp ,'DisplayName','Low pass')
plot(x_hp_lp,'DisplayName','High pass + Low pass')
However if I use the corresponding bandpass filter, the filter is unstable and the filtered time-history diverges.
[b,a] = butter(order,[fLP, fHP]/(fs/2),'bandpass');
x_bp = filter(b,a,x);
subplot(2,1,2)
plot(x_bp ,'DisplayName','Band pass')
Why does this happens and how can I correct this behaviour?
EDIT: Following Greg Dionne suggestion, I tried to switch to sos coefficients. This however didn't solve the problem.
Below is a code to check this:
[z,p,k] = butter(order,fHP/(fs/2),'high');
[sos,g] = zp2sos(z,p,k);
x_hp = filtfilt(sos,g,x);
[z,p,k] = butter(order,fLP/(fs/2),'low');
[sos,g] = zp2sos(z,p,k);
x_lp = filtfilt(sos,g,x);
x_hp_lp = filtfilt(sos,g,x_hp);
figure
subplot(2,1,1)
hold on
plot(x ,'DisplayName','original data')
plot(x_hp ,'DisplayName','High pass')
plot(x_lp ,'DisplayName','Low pass')
plot(x_hp_lp,'DisplayName','High pass + Low pass')
[z,p,k] = butter(order,[fLP, fHP]/(fs/2),'bandpass');
[sos,g] = zp2sos(z,p,k);
x_bp = filtfilt(sos,g,x);
subplot(2,1,2)
plot(x_bp ,'DisplayName','Band pass')
In this case, x_bp is all nan.

Greg Dionne on 15 May 2019
You will have better results if you use second-order sections.
See the "Limitations" section of https://www.mathworks.com/help/signal/ref/butter.html for an example of how to use them.
Greg Dionne on 17 May 2019
I see, your filter is unstable.
A recursive filter is stable when all the poles lie within the unit circle.
If you look at the poles returned, you'll see they have greater than unit magnitude.
>> [z,p,k] = butter(order,[fLP, fHP]/(fs/2),'bandpass');
>> abs(p)
ans =
1.0768
1.0768
1.1354
1.1354
1.0523
1.0523
You can also see which frequencies are going to "run away"
>> angle(p)*(fs/(2*pi))
ans =
77.3030
-77.3030
61.7558
-61.7558
51.1859
-51.1859
In your case it looks you'll get a band of frequencies that will grow with respect to time (centered at around 60 Hz).
You can also look at the pole-zero plot in fvtool if you find it more helpful to visualize them in a plot:
fvtool(b,a) %( or fvtool(sos,g))
If you click on the pole-zero plot, you'll see the following plot:
As for a simple way to proceed? I probably would use the filter designer which does all the checking for you and lets you make tradeoffs on the pass/stop bands.
filterDesigner
To see how to do this in code you can click "Generate Code" from the file button.
That gives you something like this:
Fs = 750; % Sampling Frequency
Fstop1 = 10; % First Stopband Frequency
Fpass1 = 50; % First Passband Frequency
Fpass2 = 80; % Second Passband Frequency
Fstop2 = 300; % Second Stopband Frequency
Astop1 = 60; % First Stopband Attenuation (dB)
Apass = 1; % Passband Ripple (dB)
Astop2 = 80; % Second Stopband Attenuation (dB)
match = 'stopband'; % Band to match exactly
% Construct an FDESIGN object and call its BUTTER method.
h = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, ...
Astop2, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);
From there you can get the second order sections from the "sosMatrix" and the gain from the "ScaleValues" properties.
Hope this helps!
-Greg

R2018b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!