Borrar filtros
Borrar filtros

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

41 visualizaciones (últimos 30 días)
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.

Respuesta aceptada

Greg Dionne
Greg Dionne el 15 de Mayo de 2019
You will have better results if you use second-order sections.
See the "Limitations" section of for an example of how to use them.
  4 comentarios
Luca Amerio
Luca Amerio el 17 de Mayo de 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 el 17 de Mayo de 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 =
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!

Iniciar sesión para comentar.

Más respuestas (0)




Community Treasure Hunt

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

Start Hunting!

Translated by