how to implement a non-causal system?

How do you implement a system like:

y(k+2)+2y(k+1)+1/2*y(k)=x(k+1)+x(k)

where x: input to the system and y: output of the system

What if you want to draw a block diagram for this system on paper or in simulink, how will you represent the time advancement?

Can anybody tell me what is the significance of such type of system?

[Dr.Athar Ravish Khan]

You do refer the code to test the causality of system

clc;
close all;
k=2;%delay
n=0:2+k;
x=[10 2 5 zeros(1,k)]; %x(n)
subplot(411)
stem(n,x)
xdelay=[zeros(1,k) x(1:3)]; %x(n-2)
subplot(412)
stem(n,xdelay)
y=x+n.*xdelay;  %y(n)=x(n)+n*x(n-2)
% delayed output y'(n)=x(n-k)+(n-k)*x(n-k-2)
nk=(0:length(n)-1+k)-k;
ydelayed=[xdelay zeros(1,k)]+nk.*[zeros(1,k) xdelay]  
subplot(413)
stem(0:length(ydelayed)-1,ydelayed)
n1=(0:length(n)-1+k);
ydin=[xdelay zeros(1,k)]+n1.*[zeros(1,k) xdelay]  % output due to delayed input 
subplot(414)
stem(0:length(ydin)-1,ydin)
%%ydelayed is not equal system is non causal

Answer 1

[Walter Roberson]

According to Maple's recurrence solver, if you have initial values for y(1) and y(2), then

y(k) = (1/4)*sqrt(2)*(Sum(((1/(-2+sqrt(2)))^(k-k0) * (2+sqrt(2)) + (2+sqrt(2))^(-k+k0) * (2*(-1)^(1+k-k0) + sqrt(2)*(-1)^(k-k0))) * (x(k0-2) + x(k0-1)), k0 = 3 .. k)) + (1/4*((6*y(1) + 4*y(2)) * sqrt(2) - 4*y(2) - 8*y(1))) * (-1-(1/2)*sqrt(2))^k - (3/2)*(-1+(1/2)*sqrt(2))^k * ((y(1) + (2/3)*y(2)) * sqrt(2) + (4/3)*y(1) + (2/3)*y(2))

What it means, I don't know.

[EDIT]

If you want to start with y(0) (you didn't say so I had to guess) then the formula becomes

y(k) = (1/4)*sqrt(2)*(Sum(((1/(-2+sqrt(2)))^(k-k0)*(2+sqrt(2))+(2+sqrt(2))^(-k+k0)*(2*(-1)^(1+k-k0)+sqrt(2)*(-1)^(k-k0)))*(x(k0-2)+x(k0-1)), k0 = 2 .. k))+(1/4*((-2*y(0)-2*y(1))*sqrt(2)+2*y(0)))*(-1-(1/2)*sqrt(2))^k+(1/2*((y(0)+y(1))*sqrt(2)+y(0)))*(-1+(1/2)*sqrt(2))^k

which is essentially identical in the summation but differs in the initialization.

According to this formula, the first several y(k) values starting at y(0) are:

y(0), 

y(1), x[0]+x[1]-(1/2)*y(0)-2*y(1), 

-2*x[0]-x[1]+x[2]+(7/2)*y(1)+y(0), 

(7/2)*x[0]+(3/2)*x[1]-x[2]+x[3]-(7/4)*y(0)-6*y(1), 

-6*x[0]-(5/2)*x[1]+(3/2)*x[2]-x[3]+x[4]+3*y(0)+(41/4)*y(1), 

(41/4)*x[0]+(17/4)*x[1]-(5/2)*x[2]+(3/2)*x[3]-x[4]+x[5]-(41/8)*y(0)-(35/2)*y(1),

-(35/2)*x[0]-(29/4)*x[1]+(17/4)*x[2]-(5/2)*x[3]+(3/2)*x[4]-x[5]+x[6]+(35/4)*y(0)+(239/8)*y(1), 

(239/8)*x[0]+(99/8)*x[1]-(29/4)*x[2]+(17/4)*x[3]-(5/2)*x[4]+(3/2)*x[5]-x[6]+x[7]-51*y(1)-(239/16)*y(0), 

-51*x[0]-(169/8)*x[1]+(99/8)*x[2]-(29/4)*x[3]+(17/4)*x[4]-(5/2)*x[5]+(3/2)*x[6]-x[7]+x[8]+(51/2)*y(0)+(1393/16)*y(1), 

(1393/16)*x[0]+(577/16)*x[1]-(169/8)*x[2]+(99/8)*x[3]-(29/4)*x[4]+(17/4)*x[5]-(5/2)*x[6]+(3/2)*x[7]-x[8]+x[9]-(1393/32)*y(0)-(1189/8)*y(1), 

-(1189/8)*x[0]-(985/16)*x[1]+(577/16)*x[2]-(169/8)*x[3]+(99/8)*x[4]-(29/4)*x[5]+(17/4)*x[6]-(5/2)*x[7]+(3/2)*x[8]-x[9]+x[10]+(1189/16)*y(0)+(8119/32)*y(1), 

(8119/32)*x[0]+(3363/32)*x[1]-(985/16)*x[2]+(577/16)*x[3]-(169/8)*x[4]+(99/8)*x[5]-(29/4)*x[6]+(17/4)*x[7]-(5/2)*x[8]+(3/2)*x[9]-x[10]+x[11]-(8119/64)*y(0)-(3465/8)*y(1), 

-(3465/8)*x[0]-(5741/32)*x[1]+(3363/32)*x[2]-(985/16)*x[3]+(577/16)*x[4]-(169/8)*x[5]+(99/8)*x[6]-(29/4)*x[7]+(17/4)*x[8]-(5/2)*x[9]+(3/2)*x[10]-x[11]+x[12]+(3465/16)*y(0)+(47321/64)*y(1), 

(47321/64)*x[0]+(19601/64)*x[1]-(5741/32)*x[2]+(3363/32)*x[3]-(985/16)*x[4]+(577/16)*x[5]-(169/8)*x[6]+(99/8)*x[7]-(29/4)*x[8]+(17/4)*x[9]-(5/2)*x[10]+(3/2)*x[11]-x[12]+x[13]-(47321/128)*y(0)-(40391/32)*y(1), 

-(40391/32)*x[0]-(33461/64)*x[1]+(19601/64)*x[2]-(5741/32)*x[3]+(3363/32)*x[4]-(985/16)*x[5]+(577/16)*x[6]-(169/8)*x[7]+(99/8)*x[8]-(29/4)*x[9]+(17/4)*x[10]-(5/2)*x[11]+(3/2)*x[12]-x[13]+x[14]+(275807/128)*y(1)+(40391/64)*y(0)

You can see that each of these is essentially a linear combination of y(0) with y(1).

Answer 2

[Sample Hold]

horrible what is that equation......

consider this one:

If you have y(1) and y(2) then you can get a current output for y(0) from x(0) and x(1). After that to calculate y(1), you need y(2) and y(3). What is the guarantee that the y(1) this time will not be a different from the one used before to calculate y(0)?

And, how will you draw a block diagram for this type of system?

Answer 3

[Daniel Shub]

Add a two smaple delay and your system becomes causal.

Answer 4

[Bilawal]

y(k) = (1/4)*sqrt(2)*(Sum(((1/(-2+sqrt(2)))^(k-k0) * (2+sqrt(2)) + (2+sqrt(2))^(-k+k0) * (2*(-1)^(1+k-k0) + sqrt(2)*(-1)^(k-k0))) * (x(k0-2) + x(k0-1)), k0 = 3 .. k)) + (1/4*((6*y(1) + 4*y(2)) * sqrt(2) - 4*y(2) - 8*y(1))) * (-1-(1/2)*sqrt(2))^k - (3/2)*(-1+(1/2)*sqrt(2))^k * ((y(1) + (2/3)*y(2)) * sqrt(2) + (4/3)*y(1) + (2/3)*y(2))

y(k) = (1/4)*sqrt(2)*(Sum(((1/(-2+sqrt(2)))^(k-k0)*(2+sqrt(2))+(2+sqrt(2))^(-k+k0)*(2*(-1)^(1+k-k0)+sqrt(2)*(-1)^(k-k0)))*(x(k0-2)+x(k0-1)), k0 = 2 .. k))+(1/4*((-2*y(0)-2*y(1))*sqrt(2)+2*y(0)))*(-1-(1/2)*sqrt(2))^k+(1/2*((y(0)+y(1))*sqrt(2)+y(0)))*(-1+(1/2)*sqrt(2))^k

Answer 5

[Paul]

In general (but not always), causality of a linear, time invariant (LTI) system cannot be determined by inspection of the difference equatin that defines the system's input/output relationship.

Causality can always be determined by the system's unit pulse response. The system is causal if the unit pulse response is zero for k < 0 and noncausal otherwise.

Define the variables and signals

syms k integer
syms y(k) x(k) z

Define the discrete-time unit step, assuming the default sympref for heaviside (i.e., heaviside(0) = 1/2

sympref('default');
u(k) = heaviside(k) + kroneckerDelta(k)/2;

Consider the system at hand

deqn = y(k+2) + 2*y(k+1) + 1/2*y(k) == x(k+1) + x(k)

deqn = 

If x[k] is the unit pulse, then we have

x(k) = kroneckerDelta(k);
deqn = subs(deqn)

deqn = 

Taking the bilateral z-transform of both sides, and solving for Y(z) (and renaming it to H(z) as is common of the unit pulse response) we have

H(z) = (z+1)/(z^2 + 2*z + 1/2)

H(z) = 

Compute the poles of H(z) to use later

[~,den] = numden(H(z));
poles = solve(den);
[~,index] = sort(abs(double(poles)));
poles = double(poles(index)).';

The partial fraction expansion of H(z) is

H(z) = partfrac(H(z),'FactorMode','full')

H(z) = 

Plot the poles in the z-plane. The colored regions are three possible regions of convergence (ROCs) for H(z) where the blue region extends to infinity and the boundaries between red and green and between green and blue are not included, i.e., the ROCs are open (the origin is included in the red).

figure
pr = radiusregion(abs(double([0,poles;poles,3])));
pr(1).FaceColor = 'r';pr(2).FaceColor = 'g'; pr(3).FaceColor = 'b';
pr(2).EdgeColor = 'm';pr(2).LineWidth = 1.5;pr(2).LineStyle = '--';
hold on
polarplot(angle(poles),abs(poles),'kx','MarkerSize',10)

The unit pulse response, h[k], can be one of three possbilities depending on what we assume for the ROC.

c = children(H(z));
Ha(z) = c{1}, Hb(z) = c{2}

Ha(z) = 

Hb(z) = 

h1(k) = iztrans(Ha(z),z,k)*u(k)        + iztrans(Hb(z),z,k)*u(k);              % blue
h2(k) = subs(iztrans(Ha(1/z),z,k)*u(k) + iztrans(Hb(1/z),z,k)*u(k),k,-k);      % red
h3(k) = iztrans(Ha(z),z,k)*u(k)        + subs(iztrans(Hb(1/z),z,k)*u(k),k,-k); % green

If we substitute h[k] into the LHS of the difference equation, we should obtain the RHS

deqnlhs = [
    subs(lhs(deqn),[y(k+2),y(k+1),y(k)],[h1(k+2),h1(k+1),h1(k)]);
    subs(lhs(deqn),[y(k+2),y(k+1),y(k)],[h2(k+2),h2(k+1),h2(k)]);
    subs(lhs(deqn),[y(k+2),y(k+1),y(k)],[h3(k+2),h3(k+1),h3(k)])];

deqnlhs doesn't simplify well symbolically (at least for me), but we can show via evaluation that all three solutions satisfy the difference equation

kval = -10:10;
[kval;
 simplify(subs(deqnlhs,k,kval)) - subs(rhs(deqn),k,kval)]

ans = 

Plot all three unit pulse responses

figure
stem(subplot(311),kval,h1(kval));ylabel('h1[k]');axis padded
stem(subplot(312),kval,h2(kval));ylabel('h2[k]');axis padded
stem(subplot(313),kval,h3(kval));ylabel('h3[k]');axis padded

Only h1[k] = 0 for k < 0, so the LTI system described by the original difference equation is causal only if we take h1[k], which corresponds to the blue ROC, as its unit pulse response.

If we assume the system is causal, then we solve the system iteratively forward in time given two auxilliary conditions. For example, we can specify y[0] and y]1], then iterate forward starting from y[2] for some input x[k]

For example, suppose x[k] = 0.8^k * u[k]

xval = 0.8.^(0:10);

with auxilliary conditions y[0] = 5 and y[1] = 10

yval(0 + 1) = 5; yval(1 + 1) = 10;  % explicitly showing the compensation for 1-based indexing in Matlab

The iteration loop to solve for y[2] to y[10] is (w/o preallocating y for this simple problem)

for kk = 0:8
    yval(kk+2 + 1) = -2*yval(kk+1 + 1) - 1/2*yval(kk + 1) + xval(kk+1 + 1) + xval(kk + 1);
end

Plot the output and show that it satisfies the auxilliary conditions and the difference equation

kval = 0:10;
figure
stem(kval,yval),xlabel('k');ylabel('y[k]');

yval((0:1) + 1)

ans = 1×2

5 10

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

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

yval((0:8)+2 + 1) + 2*yval((0:8)+1 + 1) + 1/2*yval((0:8) + 1) - xval((0:8)+1 + 1) - xval((0:8) + 1)

ans = 1×9

1.0e-13 * 0.0067 -0.0244 -0.0455 -0.0577 -0.0178 -0.6400 0.9653 0.0333 0.3675

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

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

We can also solve this causal system symbolically via the z-transform

syms x(k)
deqn = lhs(deqn) == x(k+1) + x(k);
x(k) = 0.8^k*u(k);
deqn = subs(deqn);
zeqn = ztrans(deqn);
syms Y(z)
zeqn = subs(zeqn,[y(0),y(1),ztrans(y(k))],[5,10,Y(z)]);
Y(z) = rhs(isolate(zeqn,Y(z)))

Y(z) = 

y(k) = iztrans(Y(z),z,k)

y(k) = 

Verify that the symbolic solution matches the numerical solution found through iteration

double(y(0:10))-yval

ans = 1×11

1.0e-12 * 0 0 0 0.0071 0 0.0142 -0.0284 0.1137 -0.2274 0.4547 -0.9095

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

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