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

69 ビュー (過去 30 日間)
Luca Amerio
Luca Amerio 2019 年 5 月 15 日
コメント済み: Greg Dionne 2019 年 5 月 17 日
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);
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);
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);
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);
plot(x_bp ,'DisplayName','Band pass')
In this case, x_bp is all nan.


Greg Dionne
Greg Dionne 2019 年 5 月 15 日
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.
  4 件のコメント
Greg Dionne
Greg Dionne 2019 年 5 月 17 日
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 =
You can also see which frequencies are going to "run away"
>> angle(p)*(fs/(2*pi))
ans =
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.
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!


その他の回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by