Use of syms, piecewise & conv commands

I'm trying to define the values of x(t) at various values of t to plot the input h(t) and x(t) vs the output y(t) after performing convolution. I am getting the a error when trying to run my code and am hoping one of the experts on here may be able to assist. I'm still quite new to MATLAB, I've looked at the associated help documetation, however the error code points me to the inbuilt conv() function which is way over my head. My code is below, along with the error message I am receiving.

Thanks in advance for any help.

Code

t = [0:1:10] % define range of t values

t = 1×11

0 1 2 3 4 5 6 7 8 9 10

h = (exp(-(t-2)).*heaviside(t-2)); % define function h(t)
%x = t .* (t >= 0) .* (t <= 1); % define function x(t) (commented out to attempt piecewise command)
syms x(t)
x(t) = piecewise((0 <= t) & (t <= 1), t, 0)

x(t) = 

y = conv(h, x); % carry out convolution

Error using indexing
Invalid argument at position 2. Symbolic function indexing evaluates the function at the input arguments. To perform colon indexing call FORMULA before indexing.

Error in conv (line 32)
c = conv2(a(:),b(:));

figure;
plot(t, h, 'r', 'LineWidth', 1, 'DisplayName', 'h(t) = e^{-(t-2)} * u(t-2)'); % plotting h(t)
hold on; % holf to allow multiple plots on one axes
plot(t, x, 'b', 'LineWidth', 1, 'DisplayName', 'x(t) = t (0 <= t <= 1)'); % plotting x(t)
plot(t, y, 'g', 'LineWidth', 1, 'DisplayName', 'y(t) = x(t) * h(t)'); % plotting y(t)
xlabel('t'); % add x-label
ylabel('Amplitude'); % add y-label
title('Plot of h(t), x(t), and y(t)'); % add title
legend('Location', 'northeast'); % define legend location
grid on; % add grid
hold off; % hold off to allow plotting of figure

Answer 1

[Paul]

Hi Ben,

First you have to decide if you want to solve the problem for a closed form expression using symbolic tools, like syms and piecewise and int, or if you want to approximate the solution numerically using ordinary numeric arrays and conv (you could also go for a numeric solution using integral). The code as posted is mixing numeric and symbolic in ways that can't be done.

Considering that you already have x(t) defined symbolically, and I'm sure you can define h(t) symbolically, I'd encourage you to use int to come up with the expression for the convolution x(t)*h(t) (unless you're required to use conv?).

If you try that and run into problems, feel free to respond with the updated code and explain what problem you're having.

[Walter Roberson]

In particular, conv() is for discrete convolution and is not able to calculate convolution integrals.

[Paul]

However, conv can be used to approximate (as stated above) samples of a convolution integral and is used to do so quite a bit on this forum (much to my disappointment, unless the purpose of a specific problem is to demonstrate that approximation), though often incorrectly.

[Ben]

Hi Paul,

Thanks for your response, the aim of the exercise is to use the conv() command, however I was struggling to work out if my original code had correctly defined x(t) and the MATLAB help documentation seemed to point towards using piecewise as the best way to define x(t) in the correct format that I needed.

I believe we are only attempting to approximate the solution as in further parts of the question we increase the number of data points in the time vector to visualise the effect is has on the output plot.

I think I may have been going down a rabbit hole trying to introuduce symbolic tools this early in our course perhaps.

Appreciate your response.

Cheers

[Paul]

Using conv, your original code was on the right track, except that the time vector and x(t) and h(t) likely needs to be sampled at a smaller spacing (though I gather you're going to play around with the spacing). Also, there's a minor issue with using heaviside for numerical solutions, so easier to just avoid it like you did in the commented version of x.

t = 0:.01:10; % define range of t values, small spacing
%h = (exp(-(t-2)).*heaviside(t-2)); % define function h(t)
h = exp(-(t-2)).*(t >= 2); % define function h(t)
x = t .* (t >= 0) .* (t <= 1); % define function x(t) (commented out to attempt piecewise command)
figure
plot(t,x,t,h),grid

Another issue to consider here is that h(t) only approaches zero as t gets large, but because we don't have infinite memory we'll have to cut it off at some value of t where h(t) is small enough. For our time vector, the last value of h used for computation is

h(end)

ans = 3.3546e-04

That looks pretty small, but if it's not small enough you can extend the t vector to a larger value.

If those plots look correct, then what's the next step?

[Ben]

Thanks Paul,

Seems to be running fine now and your plot above looks exactly like the output I am getting, I am interested to know what the issue with heaviside is though, as I can't seem to see a difference between using it and not, I imagine it will be one of those things that will only show up on occasion?

Cheers

t = [0:0.01:10] % define range of t values

t = 1×1001

0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700 0.1800 0.1900 0.2000 0.2100 0.2200 0.2300 0.2400 0.2500 0.2600 0.2700 0.2800 0.2900

%h = (exp(-(t-2)).*heaviside(t-2)); % define function h(t)
h = exp(-(t-2)).*(t>=2); % define function h(t)

x = t .* (t >= 0) .* (t <= 1); % define function x(t)

y = conv(h,x,'same'); % carry out convolution

figure;
plot(t, h, 'r', 'LineWidth', 1, 'DisplayName', 'h(t) = e^{-(t-2)} * u(t-2)'); % plotting h(t)
hold on; % holf to allow multiple plots on one axes
plot(t, x, 'b', 'LineWidth', 1, 'DisplayName', 'x(t) = t (0 <= t <= 1)'); % plotting x(t)
plot(t, y, 'g', 'LineWidth', 1, 'DisplayName', 'y(t) = x(t) * h(t)'); % plotting y(t)
xlabel('t'); % add x-label
ylabel('Amplitude'); % add y-label
ylim([0,2]); % set y-axis limits
title('Plot of h(t), x(t), and y(t)'); % add title
legend('Location', 'northeast'); % define legend location
grid on; % add grid
hold off; % hold off to allow plotting of figure

[Paul]

Hi Ben,

The solution you've posted should be reconsidered for a couple of reasons. The use of the 'same' option results in a partial solution, at best. Also, if using 'same', the output values of conv do not correspond to the values in t. Finally, assuming that the goal is approximate the convolution integral, the output of conv has to be multiplied by the sampling period.

First, let's get the exact solution.

Define the signals

syms t real
syms h(t) x(t)
h(t) = exp(-(t-2))*heaviside(t-2);
x(t) = t*rectangularPulse(0,1,t);

Plot them

figure
fplot([x(t),h(t)],[0 10]),grid
xlabel('t')

Compute the convolution integral, taking advantage of the fact that h(t) = x(t) = 0 for t < 0

syms tau real
syms y(t)
y(t) = int(h(tau)*x(t-tau),tau,0,t)

y(t) = 

Plot the result

figure
fplot(y(t),[0 20]),grid
xlabel('t')
hold on

Now take samples of x(t) and h(t) and use conv with the default option 'full' and multiply by the sampling period

t = [0:0.01:10]; % define range of t values
T = t(2);  % sampling period
%h = (exp(-(t-2)).*heaviside(t-2)); % define function h(t)
h = exp(-(t-2)).*(t>=2); % define function h(t)
x = t .* (t >= 0) .* (t <= 1); % define function x(t)
y = conv(h,x)*T; % carry out convolution

y has more samples than both x and h, so we have to create an associated time vector

ty = (0:numel(y)-1)*T;

Add to plot (show only a subset of the points in so we can clearly see the overlay).

plot(ty,y,'o','MarkerIndices',1:10:numel(y))

That looks like a pretty good approximation. But, if you zoom in tight you'll see it's not exact, though it could be improved with a smaller sampling period. Also, if you zoom in around t = 10 seconds, you'll see the effect of cutting off the samples of h(t) at t = 10. You could extend the t vector to a larger value, but there will always be a price to pay for using a finite number of samples of h to approximate a signal that lasts forever.

If we use the 'same' option, we get an output that has the same number of elements as the first argument to conv, which in this case is h, which in this problem has the same number of elements as t.

ysame = conv(h,x,'same')*T;
[numel(ysame) numel(h) numel(t)]

ans = 1×3

1001 1001 1001

But the conv output with 'same' doesn't correspond to the values in t. Rather, ysame is the 'central' part of the convolution (unfortunately the doc page conv doesn't define 'central'), which in this case is basically a 10 second interval in the middle of 0 <= t < = 20, which is 5 <= t <= 15

plot(t+5,ysame,'kx','MarkerIndices',1:20:numel(ysame))

legend('h(t)*x(t)','conv(full)','conv(same)')

Finally, the problem with heaviside is that, with the default sympref, it evaluates to 1/2 when the input evaluates to 0

heaviside(0)

ans = 0.5000

Typically, that's not what we want when using discrete-time processing to approximate continuous time results (though sometimes it may be).