Discrete filter bode shows wrong results

3 views (last 30 days)
Vorname
Vorname on 4 Dec 2023
Edited: Paul on 4 Dec 2023
I have a discrete Filter consisting of the difference of two FIR filters. One filter is multiplied with 0.5, the other is delayed by z^-1. But at the end the bode diagram seems not correct and also when I let myself display the resulting function, it is not what I have calculated. Where is my fault?
Ts = 1;
FIR1 = tf([-0.1675, 0.1750, 0.3500, -0.2250, -0.1675, -0.1250], 1, Ts);
FIR2 = tf([0.1250, -0.5850, 0.8500, -0.3000, 0.0500, -0.5850, -0.1250], 1, Ts);
H_1 = FIR1 * tf(1, [1 0], Ts);
H_2 = 0.5 * FIR2;
H_diff = H_1 - H_2;
figure;
bode(H_diff)

Answers (2)

Star Strider
Star Strider on 4 Dec 2023
The fault, if there is one, is in not understanding the Control System Toolbox conventions. Addition creates a system with the two systems in parallel (multiplication creates them in series). If you want to subtract the bode results from each other, it is necessary to get the actual results first.
I believe this is what you want.
If it is not, please describe it in detail —
Ts = 1;
FIR1 = tf([-0.1675, 0.1750, 0.3500, -0.2250, -0.1675, -0.1250], 1, Ts);
FIR2 = tf([0.1250, -0.5850, 0.8500, -0.3000, 0.0500, -0.5850, -0.1250], 1, Ts);
H_1 = FIR1 * tf(1, [1 0], Ts);
H_2 = 0.5 * FIR2;
% H_diff = H_1 - H_2;
[mag2,phase2,wout2] = bode(H_2);
[mag1,phase1,wout1] = bode(H_1, wout2);
H_1c = squeeze(mag1) .* exp(1j*deg2rad(squeeze(phase1)));
H_2c = squeeze(mag2) .* exp(1j*deg2rad(squeeze(phase2)));
H_diff = H_1c - H_2c;
figure
tiledlayout(2,1)
nexttile
semilogx(wout1, mag2db(abs(H_diff)))
grid
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')
nexttile
semilogx(wout1, rad2deg(unwrap(angle(H_diff))))
grid
xlabel('Frequency (rad/s)')
ylabel('Phasee (deg)')
sgtitle('Bode Plot of ‘H\_diff’')
.
  1 Comment
Star Strider
Star Strider on 4 Dec 2023
MATLAB is consistent.
Checking with the Symbolic Math Toolbox —
Ts = 1;
FIR1 = tf([-0.1675, 0.1750, 0.3500, -0.2250, -0.1675, -0.1250], 1, Ts);
FIR2 = tf([0.1250, -0.5850, 0.8500, -0.3000, 0.0500, -0.5850, -0.1250], 1, Ts);
H_1 = FIR1 * tf(1, [1 0], Ts)
H_1 = -0.1675 z^5 + 0.175 z^4 + 0.35 z^3 - 0.225 z^2 - 0.1675 z - 0.125 ----------------------------------------------------------------- z Sample time: 1 seconds Discrete-time transfer function.
H_2 = 0.5 * FIR2
H_2 = 0.0625 z^6 - 0.2925 z^5 + 0.425 z^4 - 0.15 z^3 + 0.025 z^2 - 0.2925 z - 0.0625 Sample time: 1 seconds Discrete-time transfer function.
% H_diff = H_1 - H_2
H_1.Numerator
ans = 1×1 cell array
{[-0.1675 0.1750 0.3500 -0.2250 -0.1675 -0.1250]}
H_2.Numerator
ans = 1×1 cell array
{[0.0625 -0.2925 0.4250 -0.1500 0.0250 -0.2925 -0.0625]}
syms z
H_1s = vpa(poly2sym(H_1.Numerator{:},z), 5)
H_1s = 
H_2s = vpa(poly2sym(H_2.Numerator{:},z), 5)
H_2s = 
H_diffs = vpa(H_1s - H_2s, 5)
H_diffs = 
I trust these results.
.

Sign in to comment.


Paul
Paul on 4 Dec 2023
FIR filters are typically specified with coefficients in ascending powers of z^-1, but the the default for tf in the Control System Toolbox specifies the transfer function numerator and denominator in terms of descending powers of z.
Ts = 1;
FIR1 = tf([-0.1675, 0.1750, 0.3500, -0.2250, -0.1675, -0.1250], 1, Ts)
FIR1 = -0.1675 z^5 + 0.175 z^4 + 0.35 z^3 - 0.225 z^2 - 0.1675 z - 0.125 Sample time: 1 seconds Discrete-time transfer function.
FIR2 = tf([0.1250, -0.5850, 0.8500, -0.3000, 0.0500, -0.5850, -0.1250], 1, Ts)
FIR2 = 0.125 z^6 - 0.585 z^5 + 0.85 z^4 - 0.3 z^3 + 0.05 z^2 - 0.585 z - 0.125 Sample time: 1 seconds Discrete-time transfer function.
Perhaps you really meant to specify FIR1 and FIR2 in terms of z^-1, like so
FIR1 = tf([-0.1675, 0.1750, 0.3500, -0.2250, -0.1675, -0.1250], 1, Ts,'Variable','z^-1')
FIR1 = -0.1675 + 0.175 z^-1 + 0.35 z^-2 - 0.225 z^-3 - 0.1675 z^-4 - 0.125 z^-5 Sample time: 1 seconds Discrete-time transfer function.
FIR2 = tf([0.1250, -0.5850, 0.8500, -0.3000, 0.0500, -0.5850, -0.1250], 1, Ts,'Variable','z^-1')
FIR2 = 0.125 - 0.585 z^-1 + 0.85 z^-2 - 0.3 z^-3 + 0.05 z^-4 - 0.585 z^-5 - 0.125 z^-6 Sample time: 1 seconds Discrete-time transfer function.
And then
H_1 = FIR1 * tf(1, [1 0], Ts);
H_2 = 0.5 * FIR2;
H_diff = H_1 - H_2;
figure;
bode(H_diff)
If that's not the issue, it would be helpful to know why you think your Bode diagram is not correct and/or what the difference is between the resulting function (is that H_diff?) and what you calculated.
  2 Comments
Vorname
Vorname on 4 Dec 2023
Edited: Vorname on 4 Dec 2023
Thank you a lot for this answer! Yes, I really mean a form like this: FIR1 = 0.125 - 0.585 z^-1 + 0.85 z^-2 - 0.3 z^-3 + 0.05 z^-4 - 0.585 z^-5 - 0.125 z^-6.
But I am still struggeling. as I calculated by hand, the resulting transfer function should be
(Fir1 delated by z^-1, FIR2 multiplied with 0.5, both get subtracted)
Paul
Paul on 4 Dec 2023
Edited: Paul on 4 Dec 2023
That actually looks like the definition of FIR2, not FIR1. I showed you how to enter that correctly. Presumably you can do the same for FIR1, which i also showed how to enter correctly.
Whatever you wanted to show for the "resulting transfer function" did not show up in your comment, so I don't know what you think the resulting transfer function should be. If it would help, you can insert an image file (.png, .jpeg, etc.) into a comment using the Image (leftmost) icon on the Insert menu.

Sign in to comment.

Categories

Find more on Get Started with Control System Toolbox in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

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

Start Hunting!