Find the impulse response

its given: 𝑦(𝑛)=𝑥(𝑛−1)−2.2𝑥(𝑛−2)+𝑥(𝑛−3)−17𝑦(𝑛−1)+3.123𝑦(𝑛−2) and 𝑥(𝑛)=(0.9)𝑛[u(n)-u(n-8)]

find h(n) as an expression.

it doesnt work when i try this:

syms z

H = (z^-1 - 2.2 * z^-2 + z^-3)/(1 + 0.1429 * z^-1 - 3.123 * z^-2);

iztrans(H)

[Paul]

What exactly doesn't work?

Why is there a coefficient of 0.1429 in the denominator of H? Shouldn't that be 17?

[Marina Petani]

Thank you. That was a typo. I should have typed 1/7.

The impulsive response turned out to be a very long expression and I am not sure if that is a correct answer. I was wondering if there was another solution for this problem without finding H(z) in the first place?

i dont understand where do i use the expression of x(n).

Answer 1

[Paul]

Hi Marina,

Sometimes the Symbolic Math Toolbox works in mysterious ways and it is difficult to manipulate expressions into simpler forms. I do think that the result of your code is correct. However, a simpler form of the result can be obtained as follows,

Find the z-transform of the unit pulse response by hand

syms z
syms n integer
syms y(n) x(n)
H(z) = (z^-1 - 2.2 * z^-2 + z^-3)/(1 + 0.1429 * z^-1 - 3.123 * z^-2)

H(z) = 

Let's check the roots of the denominator

[num,den] = numden(H(z));
solve(den)

ans = 

The roots are all real, so we should expect the solution have terms in a^n, where a is root of den (the root at 0 becomes a unit pulse). Note that the roots are outside the unit circle so the pulse response is unstable.

Break up H(z) into partial fractions

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

Now take the inverse z-transform

h(n) = iztrans(H(z))

h(n) = 

The result is hard to read, but we can use vpa to get an approximation of the exact solution to see its form

vpa(h(n),5)

ans = 

The solution has the expected form.

We can verify the exact solution by plotting the first 11 outputs against the solution from the Control Systems Toolbox

htf = tf([0 1 -2.2 1],[1 0.1429 -3.123],-1,'Variable','z^-1')

htf =
 
     z^-1 - 2.2 z^-2 + z^-3
  ----------------------------
  1 + 0.1429 z^-1 - 3.123 z^-2
 
Sample time: unspecified
Discrete-time transfer function.

figure
impulse(htf,10);
hold on
stem(0:10,double(h(0:10)),'r')

As an aside, instead of computing the transfer function by hand, we can derive it using the tools at hand. It takes a few more steps, but mitigates the possibility of errors

x(n) = kroneckerDelta(n,0);
eq = y(n) == x(n-1) - 2.2*x(n-2) + x(n-3) - 0.1429*y(n-1) + 3.123*y(n-2)

eq = 

eq = ztrans(eq)

eq = 

syms Y
eq = subs(eq,[ztrans(y) y(-2:-1)],[Y zeros(1,2)]);
H1(z) = rhs(isolate(eq,Y))

H1(z) = 

Verify that H1(z) is the same as H(z)

simplify(H1 - H)

ans(z) = 

0

[Marina Petani]

Thank you a lot. That was really helpfull.