Sum in two variables for response of NON-LTI discrete time system

Hi everyone.

I'm quite new to MATLAB and I'm trying to write a code that's able to perform the sum written below, to compute the response of a NON-LTI discrete-time system (so it can't be used the function conv()).

I cannot figure how to define the two functions and how to iterate on them to get the correct result.

This is one of the many tries I made, defining two anonymus functions: the result is always zero, which is of course wrong!

delta = 0.01;
nmax= 30;
n = -nmax:1:nmax;
kmax=30;
k=-kmax:1:kmax;
%====================================================================
x = @(n) (1/3).^(n).*heaviside(n+delta);        %Input signal
h = @(n, k) (1/2).^(n-k).*heaviside(n+delta);   %Impulse response
%====================================================================
yn = 3*(0.5).^n.*heaviside(n+delta);            %Expected response

ni=1:length(n);
for ki=1:length(k)
    xa = x(ki);
    ha = h(ni, ki);
    y(ni)=sum(xa.*ha);                          %Response (WRONG!)
end
%====================================================================
stem(n,yn,'k*');
hold on
stem(n, y, 'r:');
hold off

If anyone has an idea to make this code work, I'd be very grateful!

Thanks in advandce.

Answer 1

[Paul]

Hi @ALESSIO MANIÀ

1. The loop should be over n, beause each time through the loop it should be computing y(ni), i.e., a single element of y.

for ni = 1:length(n)

2. Internal to the loop, it should compute xa and ha for all of the values of k.

3. But ha should only be computed for the loop value of n

xa = x(k);
ha = h(n(ni),k)

4. I assume the code uses delta to deal with heaviside(0) = 1/2. FWIW, the value of heaviside(0) can be controlled with sympref, though I think the code is fine as is.

[ALESSIO MANIÀ]

Thanks a lot @Paul, your answer solves all the issues I was trying to manage perfectly!

And I really appreciated also the tip about the heaviside function, which as you've seen I tried to manage in a not very elegant way!

[Paul]

Glad to help. Actually, I thought that was pretty clever how you managed heaviside(). And that's a subtle problem with heaviside() that's easy to miss. For that reason, I avoid using heaviside() in discrete-time problems and instead just define my own function for the discrete-time unit step.

The loop approach very clearly shows the algorithm, but if you're interested it can be replaced with a single line of code using "implicit expansion," which is discussed in the documentation

nmax= 30;
n = -nmax:1:nmax;
kmax=30;
k=-kmax:1:kmax;
%====================================================================
dunitstep = @(n) double(n >= 0);          %Discrete unit step
x = @(n) (1/3).^(n).*dunitstep(n);        %Input signal
h = @(n, k) (1/2).^(n-k).*dunitstep(n);   %Impulse response
%====================================================================
yn = 3*(0.5).^n.*dunitstep(n);            %Expected response

y = sum(x(k).*h(n(:),k),2);               %Computed response
%====================================================================
stem(n,yn,'k*');
hold on
stem(n, y, 'r:');
hold off