MATLAB Answers

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

17 views (last 30 days)
Luca Amerio
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.

  0 Comments

Sign in to comment.

Accepted Answer

Greg Dionne
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.

  4 Comments

Show 1 older comment
Greg Dionne
Greg Dionne on 16 May 2019
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
Luca Amerio on 17 May 2019
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
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:
PoleZero.png
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
PoleZero.png
PoleZero.png
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

Sign in to comment.

More Answers (0)

Sign in to answer this question.


Translated by