recreating in matlab Butterworth Filter filter response

1 view (last 30 days)
There is a manual which presents a filter response. In the video they present a formula and a plot of the response. However when I tried to implement it in MATLAB I get a totally different plot. Where did I go wrong implementing this formula?
plots and code of the blog and my impelentation are attached.
Thanks.
clc
clear all
s=0:0.01:50;
H=1./((s/20.01).^4+2.6131*(s/20.01).^3+3.4142*(s/20.01).^2+2.6131*(s/20.01)+1);
plot(s,20*log10(abs(H)))

Accepted Answer

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 30 Dec 2023
Here is the complete corrected code (figure 2 is from your code which is correct):
w1 = 20.01;
w2 = 24.36;
H1=tf(1,[1/(w1^4), 2.6131/(w1^3), 3.4142/(w1^2), 2.6131/w1, 1]);
H2=tf(1,[1/(w2^4), 2.6131/(w2^3), 3.4142/(w2^2), 2.6131/w2, 1]);
figure
w = linspace(0,50,200);
[MAG1,~] = bode(H1,w);
[MAG2,~] = bode(H2,w);
MAG1 = squeeze(MAG1);
MAG2 = squeeze(MAG2);
plot(w, 20*log10(MAG1), 'b-', w, 20*log10(MAG2), 'r-','LineWidth', 2)
xlabel('\omega, [rad/s]')
ylabel('Freq. Response, [dB]')
legend('@ \omega_1 = 20.01 [rad/s]', '@ \omega_2 = 24.36 [rad/s]', 'Location', 'Best')
grid on
figure
s=0:0.01:50;
H1=1./((s/w1).^4+2.6131*(s/w1).^3+3.4142*(s/w1).^2+2.6131*(s/w1)+1);
H2=1./((s/w2).^4+2.6131*(s/w2).^3+3.4142*(s/w2).^2+2.6131*(s/w2)+1);
plot(s,20*log10(abs(H1)), 'r', s, 20*log10(abs(H2)), 'b', 'LineWidth', 2)
xlabel('\omega, [rad/s]')
ylabel('20*log10|H(\omega)|')
legend('@ \omega_1 = 20.01 [rad/s]', '@ \omega_2 = 24.36 [rad/s]', 'Location', 'Best')
grid on
  3 Comments
Sulaymon Eshkabilov
Sulaymon Eshkabilov on 30 Dec 2023
(1) "~" in [MAG1,~] = bode(H1,w) means skip phase values
(2) Both plots will be the same if s = 1i*w is used:
w1 = 20.01;
w2 = 24.36;
figure
s=0:0.01:50;
H1=1./((1i*s/w1).^4+2.6131*(1i*s/w1).^3+3.4142*(1i*s/w1).^2+2.6131*(1i*s/w1)+1);
H2=1./((1i*s/w2).^4+2.6131*(1i*s/w2).^3+3.4142*(1i*s/w2).^2+2.6131*(1i*s/w2)+1);
plot(s,20*log10(abs(H1)), 'r', s, 20*log10(abs(H2)), 'b', 'LineWidth', 2)
xlabel('\omega, [rad/s]')
ylabel('Freq Response, [dB]')
legend('@ \omega_1 = 20.01 [rad/s]', '@ \omega_2 = 24.36 [rad/s]', 'Location', 'Best')
grid on

Sign in to comment.

More Answers (1)

Chunru
Chunru on 30 Dec 2023
Edited: Chunru on 30 Dec 2023
It seems that there is a confusion in s-domain and omega domain.
The following is the Laplace Transform in s-domain. The plotting is for real value of s.
s=(0:0.01:2*pi)*20.01;
H=1./((s/20.01).^4+2.6131*(s/20.01).^3+3.4142*(s/20.01).^2+2.6131*(s/20.01)+1);
plot(s,20*log10(abs(H)));
xlabel("s")
ylabel("20*log10(abs(H(s))");
The following is in omega domain () which is related to the frequency response and it is what one would expect.
omega = (0:0.01:2*pi)*20.01;
s = 1i*omega;
H=1./((s/20.01).^4+2.6131*(s/20.01).^3+3.4142*(s/20.01).^2+2.6131*(s/20.01)+1);
plot(omega, 20*log10(abs(H)));
xlabel("\omega")
ylabel("20*log10(abs(H(j\omega))")

Community Treasure Hunt

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

Start Hunting!