 Solution to discrete time system

A Discrete Time System represented by a linear difference

equation as

y(n+N) + a1y(n+N-1) + a2y(n+N-2) +….+ aN-1y(n+1) + aNy(n)

= bo x(n+M) + b1x(n+M-1) + b2x(n+M-2) +….+ bM-1 x(n+1) + bMx(n)

Where,
n ---------> running variable (integer values)
N ----------> order of the system -------> number count for
advance operators
OR replacing n with n-N

y(n) + a1y(n-1) + a2y(n-2) +….+ aN-1y(n-N+1) + any(n-N)

= bo x(n-N+M) + b1x(n-N+M-1) + b2x(n-N+M-2) +….+ bM-1 x(n-N+1) +

bMx(n-N)

OR

(EN + a1EN-1 + a2EN-2 +….+ aN-1E + aN) y(n) = (boEm+ b1Em-1 +

b2Em-2 +….+ bM-1E + bM) x(n) ------ (1)
1
Where, E -----> is the advance operator -------- > n+1
ai’s, are the feedback coefficients, bi’s are the feed forward
coefficients ----- the values of these depends on the system
parameter values, in various combinations.
y(n) is the output sequence; x(n) is the input sequence
Generally (practically) M≤N
Equation (1) can be written as , Q(E) y(n)=P(E) x(n)

Eg. 1 : y[n+2] – y[n+1] + 0.24y[n] = x[n+2] – 2x[n+1] OR

y[n] – y[n-1] + 0.24y[n-2] = x[n] – 2x[n-1] OR

(E2 – E + 0.24) y[n] = (E2 – 2E)x[n] ; M=N

Right hand side does not contain instantaneous variable x[n].

Eg. 2 : y[n+2] – 0.16 y[n+1] + 0.5y[n] = x[n+1] + 0.5x[n] OR

y[n] – 0.16 y[n-1] + 0.5y[n-2] = x[n-1] +0.5 x[n-2] OR

(E2 – 0.16E + 0.5) y[n] = (E +0.5)x[n] ; M = N -1

Right hand side contain instantaneous variable x[n].
 Procedural aspects of the solution is same as discussed in the
continuous time case.

 Solving Difference equation using Iterative procedure :

Ex: Given the difference equation,
y[n+2] – y[n+1] + 0.24y[n] = x[n+2] – 2x[n+1] with initial
conditions : y[-1] = 2 ; y[-2] = 1, and with input x[n] = n (for n ≥ 0)
Difference equation can be written as,
y[n+2] = y[n+1] – 0.24y[n] + x[n+2] – 2x[n+1]
Put n =-2, and solve for y(0)the above difference eqn.
Put n =-1, and solve for y(1)the above difference eqn.
e.t.c.
The resulting solution is not in the closed form,
instead in the sequence form.

ex02iterativediscrete1.m and

y= y= y= y= y= y= y=
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
1.7600 1.7600 1.7600 1.7600 1.7600 1.7600 1.7600
0 2.2800 2.2800 2.2800 2.2800 2.2800 2.2800
0 0 1.8576 1.8576 1.8576 1.8576 1.8576
0 0 0 0.3104 0.3104 0.3104 0.3104
0 0 0 0 -2.1354 -2.1354 -2.1354
0 0 0 0 0 -5.2099 -5.2099
0 0 0 0 0 0 -8.6974
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
To find the Zero Input Response of the DT system:
With x[n] = 0 , equation (1) becomes, Q(E) yo[n] = 0 -----> (2)
i.e. yo[n+N] + a1yo[n+N-1] +a2 yo[n+N-2] +…+ aN-1 yo[n+1]+aNy[n] = 0
OR (EN + a1EN-1 + a2EN-2 +….+ aN-1E + aN) yo[n] = 0
Eqn. (2) shows that a linear combination of yo[n] and its N successive
advance operators is zero, not at some value of n, but for all values
of time n, of the response.
Such result is possible if and only if yo[n] and all N successive
advance operations are of the same form. Otherwise their sum can
never be zero, for all values of n.
Exponential function γn , (where γ is constant) has this property.
i.e. Ek{ γn } = γn+k = γn γk
Hence the non trivial solution is of the form
yo[n]=C γn
Now, E yo[n] = C γ1γn;
E2 yo[n] = C γ2 γn; ….
EN yo[n] = C γN γn; all will have the same form.
Substituting these in equation (2),
C(γN +a1 γN-1+ a2 γN -2+….+aN-1 γ +aN) γn = 0 ------------------> (3)
A non trivial solution valid for all value of n is
γN +a1 γN-1+ a2 γN -2+….+aN-1 γ +aN = 0 ------> (4)
That is Q(γ) = 0; or (γ- γ1) (γ- γ2)…. (γ- γN) = 0,
Where γi i = 1,2,3…N are the roots of the polynomial Q(γ)
For each values of γ, equation (3) will be satisfied, so that, total
solution to equation (2) is
yo[n] = C1γ1n + C2 γ2n +…..CN γNn ------------> (5)
Here,
 Q(γ) = γN +a1 γN-1+ a2 γN -2+….+aN-1 γ +aN = 0 is the characteristic
polynomial of the system
 γ i i = 1,2,3…N are the characteristic roots of the system.
Also called Eigen values, characteristic values
 γ1n , γ2n , … γNn are the characteristic modes of the system
 solution to the equation is nothing but the linear combination of
these characteristic modes
 Constants C1,C2,….CN can be evaluated by the N initial conditions
of the system.

Note : solution eqn.4) assumes N characteristic roots

γ i ; i = 1,2,3…N are distinct

“ If there are repeated roots [some roots occurring more

than once] the form of solution is modified”
Eg: (E – a)2 y[n] = 0 :

the solution is yo[n] = C1γ1n + C2 n γ1n

If (E – a)r y[n] = 0 : the solution is

yo[n] = C1 γ1n + C2 n γ1n + C3 n2 γ1n +….+ Cr nr-1 γ1n
The characteristic modes are, γ1n , n γ1n , n2 γ1n , … , nr-1 γ1n

Note : complex roots :

γ 1, γ 2 will be combined to take complex characteristic modes .
They occur in complex conjugate pairs
If , γ1 =|γ|ejb γ2 = γ1* =|γ|e-jb
yo(t) = C1 |γ|n ejbn + C2 |γ|n e-jbn

where, C1, C2 must be conjugates, say C1 = (C/2) ejθ, C2 = (C/2) e-jθ

then , yo[n] = C |γ|n cos(bn+θ)

Example
 To find the Impulse Response of the DT system:
Consider the difference equation,

Q(E)y[n]=P(E)x[n]

 Impulse response h[n] is the system response to an impulse input

δ[n] applied at n=0 with all initial conditions zero i.e. at n < 0.

Q(E)h[n]=P(E)δ[n] ----------> 1) and

h[-1] = h[-2] = …= h[-N] = 0

The above equation can be solved iteratively also, with x[n] = δ[n]

Eg: y[n] – 0.6y[n-1] – 0.16y[n-2] = 5x[n]+0x[n-1]+0x[n-2] ;

Replacing, y[n] = h[n] ; and x[n] = δ[n] , and h[-1] = h[-2] = …= h[-N] = 0

h[n] = 0.6h[n-1] + 0.16h[n-2] + 5 δ[n]+0x[n-1]+0x[n-2] ;

For n = 0 ; h[0] – 0.6(0) – 0.16(0) = 5 Χ 1

i.e. h[0] = 5

For n = 1 ; h[1] – 0.6(5) – 0.16(0) = 5 Χ 0

i. e h[1] = 3 e.t.c.

To find the Impulse Response – General solution

 When the impulse is struck on system, “it generates energy storage”

 i.e. it creates non-zero initial conditions instantaneously within the

system at n ≥ 0.

 Although, impulse input δ[n] vanishes as n > 0, so that there is no

input for n>0, there will be response generated by these newly
created initial conditions.

 That means the impulse response has characteristic modes valid for
n>0

 At a single moment n = 0, there can be at most an impulse of some

magnitude, Hence,

h [n] = Ao δ[n] + Linear combination of characteristic modes n ≥ 0

= Ao δ[n] + { [C1γ1n + C2 γ2n +…..CN γNn]u[n] }

Let h [n] = Ao δ[n] + {yo[n]u[n]} -------------> 2)

Putting equation 2) in equation 1) i.e Q(E)h[n]=P(E)δ[n]

Q(E){ Ao δ[n] + yo[n]u[n] } =P(E)δ[n]

Since, yo[n] is made up of characteristic modes,

Q(E){ yo[n]u[n] } = 0 (satisfies ZIR) therefore,

Q(E){ Ao δ[n] } =P(E)δ[n] i.e. assuming M = N

Ao { δ(n+N) + a1δ(n+N-1) + a2δ(n+N-2) +….+ aN-1δ(n+1) + aNδ(n) }

= bo δ(n+N) + b1δ(n+N-1) + b2δ(n+N-2) +….+ bN-1 δ(n+1) + bNδ(n)

Setting n = 0 with δ[k] = 0 for k ≠ 0 and δ[k] = 1 for k = 0,

We get,

Ao { 0 + a10+ ….+ aN-10 + aN δ[0] } = bo 0 + b10+ ….+ bN-1 0 + bN δ[0])

i.e. Ao aN δ[0] = bN δ[0]

Ref : y(n+N) + a1y(n+N-1) + a2y(n+N-2) +….+ aN-1y(n+1) + aNy(n)

= bo x(n+M) + b1x(n+M-1) + b2x(n+M-2) +….+ bM-1 x(n+1) + bMx(n)

𝐛
i.e. Ao = 𝐚𝐍
𝐍

aN , bN are the coefficient of following system differential equation

(EN + a1EN-1 + a2EN-2 +….+ aN-1E + aN) y(n) =

m m-1 m-2
(boE + b1E + b2E +….+ bN-1E + bN) x(n)
𝐛
Hence, h[n] = 𝐚𝐍 δ[n] + { [C1γ1n + C2 γ2n +…..CN γNn]u[n] }
𝐍

Constants C1,C2,….Cn of yo[n] can be evaluated from the knowledge

of N values of h[n], after n ≥ 0, which again can be evaluated
iteratively knowing the impulse input at n=0 and difference equation.

Eg. 1 : y[n+2] – y[n+1] + 0.24y[n] = x[n+2] – 2x[n+1] OR

(E2 – E + 0.24E) y[n] = (E2 – 2E)x[n] ; M=N

Term corresponding to instantaneous value of input is absent

Hence, Ao = 0
 Eg. 2 : y[n+2] – 0.16 y[n+1] + 0.5y[n] = x[n+1] + 0.5x[n] OR

(E2 – 0.16E + 0.5) y[n] = (E +0.5)x[n] ; M = N-1

Term corresponding to instantaneous value of input is present,

Hence Ao = 0.5

Example

To find the Zero State Response (ZSR) to the system:

That is the response of the system to an external arbitrary
input, x[n], making the initial conditions zero.

 Here we use superposition principle.

 Expressing an arbitrary input sequence x[n], as a weighted

(strength = x[nT] ) sum of shifted several discrete impulse signals.

 i.e. x[n] = x[0]δ[n] + x[1]δ[n-1] + x[2]δ[n-2] + …

+ x[-1]δ[n+1] + x[-2]δ[n+2] + …

 𝒚[𝒏] = ∑∞
𝒎=−∞ 𝒙[𝒎]𝒉[𝒏 − 𝒎]

Where, m is the summation variable dummy integer variable.

The sequence can be viewed as shown :
System

INPUT -------------> OUTPUT

δ[n] ------------> h[n]

δ[n-m] -------------> h[n- m]

x[m] δ[n - m] --------------> x[m] h[n- m]

∑∞ ∞
𝒎=−∞ 𝒙[𝒎]𝜹[𝒏 − 𝒎] ----------> ∑𝒎=−∞ 𝒙[𝒎]𝒉[𝒏 − 𝒎]

x[n] ------------> y[n] = x[n] *h[n]

“summation of these impulse responses” called

“CONVOLUTION summation” given by

y[n] = x[n] *h[n] = ∑∞

𝒎=−∞ 𝒙[𝒎]𝒉[𝒏 − 𝒎]

= h[n] *x[n] = ∑∞
𝒎=−∞ 𝒉[𝒎]𝒙[𝒏 − 𝒎]

 To get the ZSR, we have to convolve the arbitrary input signal x[n]
with the impulse response of the system h[n]. i.e. ys[n] = x[n] * h[n]

 Applicable to Linear time invariant, non-causal OR causal

signals / systems.

 ZSR and causality

Causality restriction on both signals and systems further simplifies

the limits of summation.

Causal system’s unit impulse response h[n] is a causal signal.

 CONVOLUTION SUMMATION for a causal system and signal is

given by

y[n] = x[n] * h[n] = ∑𝒏𝒎=𝟎 𝒙[𝒎]𝒉[𝒏 − 𝒎]

properties ?
width property : If x1[n] and x2[n] have finite widths of W1 and W2,
respectively, then the width of x1[n] * x2[n] is W1 + W2. { The width in
discrete time signal is one less than the number of its elements
(length).
 Alternatively, If x1[n] and x2[n] have finite lengths of L1 and L2,
elements respectively, then the length of x1[n] * x2[n] is L1 + L2-1
elements.

Example

Total solution (response of the system) is the sum of ZIR and ZSR,
computed separately

i.e y[n] = yo[n] + ys[n] = ZIR + ZSR

 Graphical convolution
Explain with MATLAB : ex02graphicalconvolutiondiscrete1.m
t = -4 No overlap between x(τ) and h(t- τ) : hence y(-4) = 0

t = -3 there is one sample overlap between x(τ) and h(t- τ) : hence

y(-3) = 1 X 2 = 2
t = 0 there is 4 samples overlap between x(τ) and h(t- τ) : hence
y(-3) = 1 X 2 +1 X 1.41 + 1X1 + 1X0.7 = 5

t = 5 there is 9 samples overlap between x(τ) and h(t- τ) : hence

y(-3) = 1 X 2 +1 X 1.41 + 1X1 + 1X0.7 + ........ = 6.527

E.T.C.
 ex02graphicalconvolutiondiscrete2.m

 Classical solution to differential equations

Total solution, (response of the system) can be viewed as another set
of two parts

1) Natural response ---> Homogeneous solution (or Complementary

solution), which contains linear combination of Characteristic modes.

2) Forced component of response -----> Particular solution, which

contains Non – Charecteristic modes, generally response is of the
same form as forcing response (input response)

y[n] = yn[n] + yf[n]

Given the difference equation,

Q(E){yn[n] + yf[n]} = P(E)x[n]

i.e. Q(E)yn[n] + Q(E)yf[n] = P(E)x[n]

Since, yn[n] has only Linear combination of characteristic modes,

and has the same form as ZIR.

The first term in L.H.S of diff. eqn. Q(E)yn[n] = 0,

But, the arbitrary constants (coefficients of characteristic modes), are

different, from that of ZIR and is evaluated from the auxiliary
conditions.

Note: Auxiliary conditions are defined /evaluated at n ≥ 0 iteratively

i.e. given as y[0], y[1], …., y[N-1]

Now , The second term in L.H.S of diff. eqn. Q(E)yf[n] = P(E)x[n]

 To find the forced response of the system: yf[n]

The method of undetermined coefficient with following tables listing the
inputs and the corresponding forms of forced function with undetermined
coefficients.

The undetermined coefficients, are evaluated by substituting this

expression for yf[n] in Q(E)yf[n] = P(E)x[n] and equating the
coefficient of similar terms.
 For other class of inputs:

INPUT FORM OUTPUT FORM

1) rn r ≠ γi(i=1,2,3….N)(characteristic roots) βrn

2) rn r = γi(i=1,2,3….N)(characteristic roots) βnrn

3) K (a constant) β (a constant)

4) Cos(Ωn+θ) β Cos(Ωn+φ)

5) (nr+αr-1 nr-1+…..+α1 n+ αo) rn (βr nr+ β r-1 nr-1+…..+ β1n+ βo) rn

NOTE : In the classical method Zero Input and Zero State components
cannot be separated. Consequently, the initial conditions must be
applied to the total response, which begins at n= 0.

Hence, Auxiliary conditions as y[0], y[1], …., y[N-1], are defined

/evaluated at n ≥ 0 iteratively from the given the initial conditions
as y[-1], y[-2], …., y[-N]

Example:

Case 1) Exponential inputs: x[n] = rn

Forced response yf[n], to such input can be expressed in the form β rn

Hence, Q(E)yf[n] = P(E)x[n]

i.e. Q(E) β rn = P(E) rn

E β rn+1 = r rn ;

E2 βrn = r2 rn; …. ;

Er β rn = rr rn; and then,

Q(E) rn = Q(r) rn and P(E) rn = P(r) rn

Substituting these in above equation,

β Q(r) rn = P(r) rn OR
 𝐏(𝐫))
β= = H(r) Hence,
𝐐(𝐫)

Forced response yf[n] = H(r) rn for n ≥ 0 and

Total response y[n] = yn[n] + yf[n]

N
i.e. y[n] = K 
j 1
j
n
j  H (r )r n

Case 2) Constant inputs: x[n] = C

It is a special case of exponential rn with r = 1, hence,

Forced response yf[n] = C H(r) rn with r = 1 for n ≥ 0
i.e. yf[n] = C H(1)

Case 3) Exponential inputs: x[n] = ejΩn

It is a special case of exponential rn with r = ejΩ, hence,

𝐏(ejΩ ))
Forced response yf[n] = H(ejΩ) ejΩn = for n ≥ 0
𝐐(ejΩ )

Case 4) Sinusoidal inputs: x[n] = cosΩn

We know that, cosΩn = ½{ ejΩn + e-jΩn }

Hence, Forced response,
yf[n] = ½{ H(ejΩ) ejΩn + H(e-jΩ) e-jΩt } for n ≥ 0
The two terms in the RHS, which are conjugates gets cancel each
other, resulting, yf[n] = Re{ H(ejΩ) ejΩn }
j
But, H(ejΩ) = |H(ejΩ)|e jH ( e )

j
Hence, yf[n] = Re{ |H(ejΩ)| e jH ( e )  jn
}
= |H(ejΩ)| cos [Ωn+ H (e j ) ]
Hence, for generalised input, x[n] = cos(Ωn+θ) ,
yf[n] = |H(ejΩ)| cos [Ωn+ θ + H (e j ) ]

Example:
Following is the MATLAB code for iterative procedure:
for total solution (ZSR +ZIR)

n=[-2:10]';

y= [1;2;zeros(length(n)-2,1)];

%Predefining the o/p array with initial condition and

zeros

x=[0;0;n(3:end)]; % defining input array

for k=1:length(n)-2

y(k+2) = y(k+1) -0.24*y(k) + x(k+2) - 2*x(k+1);

end;

stem(n,y,'k')

O/P -----> y=[1.0000 2.0000 1.7600 2.2800 1.8576

0.3104 -2.1354 -5.2099 -8.6974 -12.4470 -16.3597 -
20.3724 -24.4461]

For finding the zero input response, do not consider the feed forward
coefficients (terms corresponding to the coefficients of x[n] and its advance
operations), in the difference equation.

Use “filter” , “filtic” and “conv” functions:

MATLAB’s filter command provides an efficient way to evaluate the system
response of a constant coefficient linear difference equation represented in
N N
delay form as  a k y[n  k ]   bk x[n  k ]
k 0 k 0

Filter requires three input arguments:

1) A length N+1 vector of feed forward coefficients [b o,b1,b2,…bN]

2) A length N+1 vector of feedback coefficients [ao,a1,a2,…aN] and
3) An input vector.

Optional 4th argument:

If a(1) is not equal to 1, FILTER normalizes the filter coefficients by a(1).

Finding Impulse Response:

For the same above example,

Following is the MATLAB code:

b= [1 -2 0]; a= [1 -1 0.24]; %defining difference

equation

n=[0:30]'; delta = inline('n==0','n'); %defining

impulse function

h= filter(b,a,delta(n)); stem(n,h,'k')

Finding Zero –State Response:

F
or the same above example,

Following is the MATLAB code:

b= [1 -2 0]; a= [1 -1 0.24]; %defining difference

equation

n=[0:30]'; x=inline('n','n') %defining input function

y= filter(b,a,x(n)); stem(n,y,'k')

Finding Zero –Input Response:

Filtic converts the traditional y[-1], y[-2],…y[-N] initial conditions for use
with the filter command

TRANSPOSED DIRECT FORM II filter structure.

An input zero is created with the zeros command

For the same above example, with same initial conditions,

zi = filtic(b,a,[2 1]);

yo= filter(b,a,zeros(size(n)),zi); stem(n,yo,'k')

Finding Total Response:

For the same above example, with same initial conditions, and with same
input,

yt= filter(b,a,x(n),zi); stem(n,yt,'k')

Convolution of two signals : (using conv command)

Convolution of two finite-duration discrete time signals is accomplished by

using the conv command.
 Eg: a) (u[n]-u[n-4)* (u[n]-u[n-4)
>
> conv([1 1 1 1],[1 1 1 1])

O/P ----> 1 2 3 4 3 2 1 Regions of support is (0≤ n≤ 6) (7

elements = 4+4-1)

b) (u[n+4]-u[n)* (u[n]-u[n-4)
>
> conv([1 1 1 1],[1 1 1 1])

O/P ----> 1 2 3 4 3 2 1 Regions of support is (-4≤ n≤ 2)

(graphically supported)

In general, the conv command cannot properly convolve infinite – duration

signals.

Conv can correctly compute a portion of such infinite duration signals.

By passing the first N samples of each, conv retuns a length (2N-1)

sequence. The first N sample are correct; the remaining N-1 samples are
not.

For the same above example, the impulse response (calculated

analytically) of the system is h[n] = {-7 (0.6)n + 8 (0.4)n}u[n], hence
defining the inline function h[n],

>>h1 = inline('(-7*0.6.^n+8*0.4.^n).*(n>=0)','n');

>>ytc= conv(h1(n),x(n));

>>stem([0:60],ytc,'k') Compare the result with

earlier method.

The results are correct over (0≤n≤30), the remaining values are clearly
incorrect; the output envelope should continue to grow, not decay.
Normally incorrect values are not displayed.

>> stem(n,ytc(1:31),'k') will be identical to the

earlier result.