How to find response of a system through convolution?

So, i've been trying to solve a frequency domain problem in the time domain, to better understand how it works and because it offers simulation perfomance for some of my algorithms.

The algorithm i wrote gets results from even complicated systems like three or more orders, but for some reason, it's always slightly different from the correct response. To measure how accurate my algorithm is, i've been using the lsim function.

s = tf('s');
t = (0:0.01:10)';

unit_impulse = zeros(size(t));
unit_impulse(t==0) = 1;
unit_step = ones(size(t));

u = unit_step;
freqFunc = 4 / (s+1);
timeFunc = 4 * exp(-t);

% Calculates reference values
y_ref = lsim(freqFunc, u, t);
h_ref = impulse(freqFunc, t);

h = conv(unit_impulse, timeFunc);
h = h(1:1001);

y = median(diff(t)) * conv(u, h);
y = y(1:1001);

error_impulse = er(h_ref, h, t);
error_response = er(y_ref, y, t);

The error in the impulse response is zero, so i must be doing something wrong when calculating the input response. The mean error is very small, like 0.02, but still enough for it to be noticeable. I divided the y_ref by y elementwise, and noticed an interesting pattern. It goes like this: [Inf, 2, 1.5, 1.333, ... , ..., 1.005]

Now that is the (X+1) / X sequence, and normally used by filters i think.

What am i doing wrong?

Answer 1

[Paul]

I don't think you're doing anything wrong. Rather, what you're seeing is the effect of approximating a convolution integral with a convolution sum scaled by dt, which is basically approximating the convolution integral with a rectangular rule. I don't know what the er() function is, but I suspect that the errors, however they are computed, will get smaller if you make the increment smaller in the time vector. Also, keep in mind that lsim() and impulse() are themselves approximations, though in this case very good approximations. Also, the step response can be computed using step(); no need to use lsim(). lsim(), step(), and impulse() aren't really needed at all. Closed form expressions for the step and impulse responses can be easily obtained and evaluated.

[AlexRD]

Oh wow, you're correct! Can't believe i missed that, thank you very much!

I guess i thought that lsim calculated the response the same way i did, so any dt errors would be the same for both. Do you reckon what lsim does is change the dt to an even smaller value and returns the interpolated results back?

Also, the er function just computes the error and plots both of the outputs.

[Paul]

You're very welcome.

See the Algorithm section on

doc lsim

for a brief discussion on how lsim works.