Plotting the convolution of u(t) * dirac(t)

Hi, I am trying to plot the convolution of u(t) "Heaviside function" and d(t) "dirac(t), or delta function".

According to the definition of convolution, u(t) * d(t) = integrate { (u(t) x d(t - Tau) ) dTau } from -oo to t.

Based on this information, I came up with below code but it is not working, any ideas how I can fix it?

t = [-1:0.01:1];
u = @(t) t>= 0;      % unit step function
r = @(t) t.*u(t);    % ramp function
d = @(t) dirac(t);


syms tau
y = int(u(t) .* d(t-tau), tau,  -Inf, t);
fplot(y, 'c');

Edit: adding the macro functions' definition for u(t), d(t), r(t)

Answer 1

[Paul]

Hi Jeremy,

The equation for the convolution integral in the Question is incorrect. The argument of u inside the integral should be tau, not t. Also, when dealing with the delta function it's better to evaluate the convolution integral over the entire real line to ensure we integrate through the delta funciton, i.e., if the uppler limit is t, then the integrand "stops" at dirac(0), which is problemlatic.

So

syms tau t real
syms u(t) d(t) y(t)
y(t) = int(u(tau)*d(t-tau),tau,-inf,inf,'Hold',true)

y(t) = 

Nwo define the signals u(t) and d(t) and evaluate y(t)

u(t) = heaviside(t);
d(t) = dirac(t);
y(t) = (release(subs(y(t))))

y(t) = 

That's a funny looking expression, but we can see that y(t<0) = 0, y(0) = 1/2, and y(t>0) = 1, which is heaviside(t) (using the default definiton for the heaviside(0)). Plot it to verify y(t) is the step function.

fplot(y(t))
ylim([-0.5 1.5])

Recall that the general rule for convolution with the dirac function is that the convolution just returns the other signal, as just shown for example, but this can get a little sticky when the other signal has discontinuous jumps, like the step function.

We can try the sympref, as @Walter Roberson did, to assume that heaviside(0) = 1 as typical for the unit step function, but I don't think it will matter. When t = 0, the convolution integral is

y0 = int(u(tau)*d(0-tau),tau,-inf,inf,'Hold',true)

y0 = 

I think the rule, at least in engineering math, for this situation (heaviside discontinuous at t =0) is that the integral evaluates to the average of heaviside(0-) and heaviside(0+), i.e., the values of heaviside(t) just to the left and right of the origin, which will always be 1/2.

[Walter Roberson]

Except that u(t) is not heaviside(t) with default sympref, u(t) is defined as 1 at t=0

[Paul]

Sure, but does it matter?

syms t tau real
sympref

ans = struct with fields:
           FourierParameters: [1    -1]
           HeavisideAtOrigin: 1/2
            AbbreviateOutput: 1
               TypesetOutput: 1
         FloatingPointOutput: 0
      PolynomialDisplayStyle: 'default'
    MatrixWithSquareBrackets: 0

int(heaviside(tau)*dirac(t-tau),tau,-inf,inf)

ans = 

sympref('HeavisideAtOrigin',1);
sympref

ans = struct with fields:
           FourierParameters: [1    -1]
           HeavisideAtOrigin: 1
            AbbreviateOutput: 1
               TypesetOutput: 1
         FloatingPointOutput: 0
      PolynomialDisplayStyle: 'default'
    MatrixWithSquareBrackets: 0

int(heaviside(tau)*dirac(t-tau),tau,-inf,inf)

ans = 

sympref('HeavisideAtOrigin',0.12345);
sympref

ans = struct with fields:
           FourierParameters: [1    -1]
           HeavisideAtOrigin: 2469/20000
            AbbreviateOutput: 1
               TypesetOutput: 1
         FloatingPointOutput: 0
      PolynomialDisplayStyle: 'default'
    MatrixWithSquareBrackets: 0

int(heaviside(tau)*dirac(t-tau),tau,-inf,inf)

ans = 

So we get the same answer in Matlab regardless of the specific value of h(0). Again, it comes down to how to integrate the product of diract(t) with a function that takes a step at t=0, if that integration is even valid at all, as briefly pointed out above.

As I said above, we could just use the general sifting property of the delta function to claim that the convolution of dirac(t) and u(t) (unit step with u(0) = 1) is u(t), but the sifting property only works when the fucntion to be sifted is contninuous at the point where the delta function is applied, so there will be problems at u(0) anyway.

[Paul]

Also, I want to point out that even though the value of heaviside(0) doesn't really matter much, it at all, for continuous time analysis, it can matter greatly if heaviside() is used to geneate a discrete-time unit step function, which must have u[n=0] = 1. However, it's not uncommon to see code like this:

t = 0:.1:1;
u = heaviside(t);
y = lsim(continuoustimetf,u,t);

This code is in error becasue lsim transforms continuoustimetf into a discrete approximation and then computes the solution in discrete time. But the result is incorrect because u(1) = 0.5 (unless the user knew to change the sympref option).

Answer 2

[Star Strider]

I am not certain what you are doing, however it is necessary to call the correct functions in order to evaluate them.

Try this —

tv = [-1:0.1:1];
syms tau t 
u(t) = heaviside(t)

u(t) = 

d(t) = dirac(t)

d(t) = 

y(t,tau) = int(u(t) .* d(t-tau), tau,  -Inf, t)

y(t, tau) = 

figure
hold on
for k = 1:numel(tv)
    t = tv(k);
    hfp(k) = fplot(y(t,tau), [-1 1], '.');
end
hold off
grid
ylim(ylim+[-1 1]*0.1)

% hfp

.

Answer 3

[Walter Roberson]

syms tau t
sympref('heavisideatorigin', 1);

u(t) = heaviside(t);
d(t) = dirac(t);
y = int(u(t) .* d(t-tau), tau,  -Inf, t)
fplot(y, [-1 1], 'c');
xlim([-1.1 1.1])

It is a pulse that is 1 at 1 and 0 elsewhere