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

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.

Answer 1

[Greg Dionne]

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.

[Luca Amerio]

Also using second-order sections, I still cannot use the bandpass filter. See my edid

[Greg Dionne]

Out of curiosity, is your data single-precision? Even second-order sections can suffer from numerical concerns, and single-precision often exacerbates this.

If you post your data, I can give it a try.

[Luca Amerio]

My data are double precision.

The issue can be reproduced also using randomly generated data.

I will attach my data as soon as possible

[Greg Dionne]

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