Why this bandpass butterworth is unstable (while the corresponding low and high pass are stable)?
19 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Luca Amerio
el 15 de Mayo de 2019
Comentada: Greg Dionne
el 17 de Mayo de 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 comentarios
Respuesta aceptada
Greg Dionne
el 15 de Mayo de 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 comentarios
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 =
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
Más respuestas (0)
Ver también
Categorías
Más información sobre Filter Design en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!