Find the impulse response
33 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Marina Petani
am 21 Mai 2022
Kommentiert: Marina Petani
am 23 Mai 2022
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)
2 Kommentare
Paul
am 22 Mai 2022
What exactly doesn't work?
Why is there a coefficient of 0.1429 in the denominator of H? Shouldn't that be 17?
Akzeptierte Antwort
Paul
am 22 Mai 2022
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)
Let's check the roots of the denominator
[num,den] = numden(H(z));
solve(den)
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))
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)
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')
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 = ztrans(eq)
syms Y
eq = subs(eq,[ztrans(y) y(-2:-1)],[Y zeros(1,2)]);
H1(z) = rhs(isolate(eq,Y))
Verify that H1(z) is the same as H(z)
simplify(H1 - H)
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!