Simulating the impulse response of a higher order filter.

I am trying to simulate an 8th order continuos time filter in MATLAB. The filter transfer function is

I am cascading 8 first order filters byusing the following code for the purpose

clear all;

clc;

wch = 2*pi*1e9;

filter_1 = tf([wch],[1 wch]);

sys = filter_1^8;

bode(sys);

figure

impulse(sys)

Somehow, the bode plot comes out fine, but the impulse response does not show anything.

If I simulate the same code with a second order filter , I can get the impulse response to show up perfectly

sys = filter_1^2

Obviously, thgere is something I am not doing correct. It would great if someone could point out the error. Thanks.

Answer 1

[Naga]

Hello Saqib,

By default, MATLAB might not simulate for a long enough period to capture the full impulse response of a high-order system. Specifying a custom time vector ensures that the simulation runs for an adequate duration. The 'linspace' function is used to create a time vector that defines the start and end times, as well as the number of points in the simulation.

clear all;
clc;

wch = 2*pi*1e9;

filter_1 = tf([wch], [1 wch]);

sys = filter_1^8;

bode(sys);

figure;

% Define a time vector for the impulse response simulation
t = linspace(0, 1e-8, 1000); % Adjust the time vector as needed

% Impulse response using state-space representation
impulse(sys, t);

[saqib]

@Paul @Naga

Thank you both for your inputs, It seems with a properly specified time vector, the impulse response is simulated properly.

I had one more question resgarding this. It seems that for both the and order filters, the time step is of the order of a 0.1 - 1ns. Why is it that the order filter can simulate properly without specifying a time vector but not the order filter?

[Paul]

Maybe you've uncovered a flaw in impulse. If the t vector is not specified, impulse is suppsoed to smartly define a t vector based on the system dynamics. But it seem like there is a problem as the system order inreases

wch = 2*pi*1e9;
filter_1 = tf([wch], [1 wch]);
for ii = 1:8
    sys = filter_1^ii;
    [y,t] = impulse(sys);
    dt(ii) = t(2);
end
format short e
dt

dt = 1×8

7.3294e-12 7.3294e-12 1.0000e-02 1.0000e-02 1.0000e-02 1.0000e-02 1.0000e-02 1.0000e-02

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The dt that impulse() uses to discretize the model is way too large for a filter order greater than 2. Same issue if using zpk or ss forms too.

Answer 2

[Paul]

Hi @saqib,

You're probably running into issues with numerical accuracy because the coefficients in the system are so large.

Generally speaking, zpk is preferred over tf.

And/or, consider rescaling the time variable, maybe from seconds to nanoseconds.

And/or, consider using the time argument to impulse to make sure it's using a small enough time step.

And/or, consider computing the impulse response of filter_1 and then get the impulse response of the system using convolution.

[Paul]

Upon further inspection, it appears that the logic that impulse uses to select the step size, if not specified by the user, is flawed when the plant has only a repeated, real pole (did not explore if there is a problem if the plant has only a repeated, complex pair of poles).

For example, with a simple, second order system

p = 300;

h = zpk([],[-p, -p],1);

[y,t] = impulse(h);

syms s
figure
fplot(ilaplace(1/(s+p)^2),[0 0.1]); % exact answer
hold on
plot(t,y),grid
xlim([0,0.1])
legend('exact','impulse');

The logic can deal with two repeated poles, as long as there is at least one simple pole

h = zpk([],[-p, -p, -1],1);

[y,t] = impulse(h);

figure
fplot(ilaplace(1/(s+p)^2/(s+1)),[0 0.1])
hold on
plot(t,y),grid
legend('exact','impulse');
xlim([0 1])

In this case, impulse() is still smart enough to compute the time step based on the two fast, repeated poles

t(2)

ans = 3.0701e-04

I'm pretty sure that step will have the same problem in selecting the step size when the plant has only a repeated, real pole.

Edit: The problem can show up when the plant has only repeated poles, whether they be real, complex, or a combination of both. There may be cases, as shown above, where the problem does not arise becasue of a successful test against a numeical tolerance.

Answer 3

[Andrew Ouellette]

Hello,

You are encountering this behavior because your cascaded system is poorly conditioned numerically. Hopefully this is not surprising to you, as your transfer function has coefficients that vary in magnitude from the order of 10^0 to 10^78. In this case, I believe the best solution would be for you to select a different time unit for your system. Let's try with nanoseconds:

wch = 2*pi; %rad/nanosecond
filter_1 = tf([wch],[1 wch],TimeUnit="nanoseconds");
sys = filter_1^8;
f = figure();
t = tiledlayout(f,1,2);
bp = bodeplot(t,sys);
bp.Layout.Tile = 1;
ip = impulseplot(t,sys);
ip.Layout.Tile = 2;

As you can see, both the bode response and impulse response are computed correctly. If you wish to view the charts in seconds, you can then alter the chart units.

f2 = figure();
t2 = tiledlayout(f2,1,2);
bp2 = bodeplot(t2,sys);
bp2.Layout.Tile = 1;
ip2 = impulseplot(t2,sys);
ip2.Layout.Tile = 2;
bp2.FrequencyUnit = "rad/s";
ip2.TimeUnit = "seconds";