> restart; with(inttrans);

DVV2 Opgaven

>

Serie 1 - Laplace

Opgave 1 - Convolutie

Bereken (x*v)(t) voor -infinity < t < infinity voor het geval dat

> x:= t-> Heaviside(t-1) + Heaviside(t-2) - 2*Heaviside(t-3):
v:= t-> 2*Heaviside(t) - Heaviside(t-1) - Heaviside(t-2):

Bedenk dat een verschuiving naar negatieve tijden bij de Laplace-transformatie niet geoorloofd is.

=> Ik heb alles dus een naar rechts geschoven.

Met laplace

> Y:= s->laplace(x(t),t,s) * laplace(v(t),t,s): 'Y(s)'=Y(s);

> y:= t->invlaplace(Y(s),s,t): 'y(t)'=y(t);

> plot([subs(t=a+1,y(t)),x(a+1),v(a+1)],a=-1..5,legend=["y","x","v"]);

Opgave 2 - Stapresponsie en uitgangssignaal uit impulsresponsie

Een lineair, tijdinvariant systeem heeft de impulsresponsie

> assume(t>=0): h:= t->exp(-t)+sin (t):

a) Bereken op twee verschillende manieren de stapresponsie g(t) voor alle t>=0

b) Bereken op twee verschillende manieren het uitgangssignaal y(t) voor alle t>=0 wanneer het ingangssignaal

> input:= t->Heaviside(t)-Heaviside(t-2):

en het systeem geen initiele energie bevat op t=0.

>

a) Bereken de stapresponsie.

Integreren

De stapresponsie kan gevonden worn door de convolutie met de eenheidsstapfunctie Heaviside te nemen. Echter,

de functie is alleen gedefinieerd voor t>=0, dus integreren vanaf t=0! (hier gaan de meesten de fout in)

> sr1(a):= int(h(t),t=0..a);

Anderzijds is de stapresponsie de integraal van de impulsresponsie en die hebben we.

> sr2(a):= int(h(t),t=0..a);

Laplace

Of met de laplace-transformatie:

laplacetransformeren, delen door s (=integreren in tijddomein) en dan terugtransformeren.

> H(s):= laplace(h(t),t,s);

> SR3(s):=H(s)/s;

> sr3(t):=invlaplace(SR3(s),s,t);

En dit geeft (gelukkig) hetzelfde als de andere berekeningen.

b) Bereken het uitgangssignaal

Het uitgangssignaal is de convolutie van de impulsresponsie met het ingangssignaal.

Convolutieproduct integreren

> output(t):=simplify(int(h(t-lambda)*input(lambda),lambda=-infinity..infinity));

Eens kijken hoe het eruit ziet.

> plot([output(t),input(t),h(t)],t=-1..30,legend=["output","input","h"]);

Convolutie met laplace

Convolutie is het product van de laplace-getransformeerden

> OUTPUT(s):=laplace(h(t),t,s) * laplace(input(t),t,s);

En nu weer naar t-domein

> output(t):=invlaplace(OUTPUT(s),s,t);

> plot([output(t),input(t),h(t)],t=-1..30,legend=["output","input","h"]);

Opgave 3 - Uitgangssignaal uit impulsresponsie

Een lineair, tijdsinvariant systeem heeft de impulsresponsie

> h:= t->sin(t)*Heaviside(t-2):

Bereken het uitgangssignaal y(t) voor alle t>=0 wanneer het ingangssignaal gelijk is aan

> x:= t->Heaviside(t)-Heaviside(t-1): assume(t>=0):

en het systeem geen initiele energie bevat.

Convolutie door integratie

> y(t):= simplify(int( h(t-lambda)*x(lambda),lambda=-infinity..infinity));

Eventjes een plaatje maken

> plot([y(t),h(t),x(t)],t=-0..15,legend=["y","h","x"]);

Convolutie met Laplace

Convolutie in het t-domein is vermenigvuldiging in het s-domein. Dus transformeren naar s-domein en vermenigvuldigen:

> Y(s):= simplify(laplace(x(t),t,s) * laplace(h(t),t,s));

En dan weer terug naar het t-domein:

> y2(t):= invlaplace(Y(s),s,t);

Nog even grafisch checken

> plot([y2(t),h(t),x(t)],t=-0..15,legend=["y","h","x"]);

Opgave 4 - Laplace-transformaties

Een tijdcontinue signaal x(t) heeft de laplace-getransformeerde:

> X:= s->(s+1)/(s^2+5*s+7);

Bepaal de laplace-getransformeerde V(s) van de volgende signalen.

Goed, lekker simpel: eerst de functie x bepalen en die verder gebruiken.

> x:=t-> invlaplace(X(s),s,t): 'x(t)'=x(t);

a)

> unassign('v(t)','V(t)');

> v:=t->x(3*t-4)*Heaviside(3*t-4):

> V(s):=simplify(laplace(v(t),t,s));

b)

> unassign('v(t)','V(t)'):

> v(t):=t*x(t):

> V(s):=simplify(laplace(v(t),t,s));

c)

> unassign('v(t)','V(t)'):

> v(t):=diff(x(t),t$2):

> V(s):=simplify(laplace(v(t),t,s));

d)

> unassign('v(t)','V(t)'):
v(t):=int(x(t),t):

> V(s):=simplify(laplace(v(t),t,s));

e)

> unassign('v(t)','V(t)'):
v(t):=x(t)*sin(2*t);

> V(s):=laplace(v(t),t,s);

>

f)

> unassign('v(t)','V(t)'):
v(t):=e^(-3*t)*x(t);

> V(s):=simplify(laplace(v(t),t,s));

g)

Dit weten we uit ons hoofd, de convolutie van x en y is het product van de laplace-getransformeerden.

> unassign('v(t)','V(t)'):

> V(s):=X(s)*X(s);

Opgave 5 - Invlaplace-transformaties

Bepaal de laplace-getransformeerde van de volgende functies

> unassign('x(t)','X(s)'):

a)

> X(s):=(s+2)/(s^2+7*s+12):
x(t)=invlaplace(X(s),s,t);

b)

> X(s):=(s+1)/(s^3+5*s^3+7*s):
x(t)=invlaplace(X(s),s,t);

c)

> X(s):=(2*s^2-9*s-35)/(s^2+4*s+2):
x(t)=invlaplace(X(s),s,t);

d)

> X(s):=(3*s^2+2*s+1)/(s^3+5*s^2+8*s+4):
x(t)=invlaplace(X(s),s,t);

e)

> X(s):=(s^2+1)/(s^5+18*s^3+81*s):
x(t)=invlaplace(X(s),s,t);

f)

> X(s):=(s+exp(-s))/(s^2+s+1):
x(t)=invlaplace(X(s),s,t);

g)

> X(s):=s/(s+1)+(s*exp(-s)+exp(-2*s))/(s^2+2*s+1):
x(t)=invlaplace(X(s),s,t);

Opgave 6 - Invlaplace

a)

> X(s):=(s^2-2*s+1)/(s*(s^2+4)):
x(t)=invlaplace(X(s),s,t);

b)

> X(s):=(s^2-2*s+1)/(s*(s^2+4)^2):
x(t)=invlaplace(X(s),s,t);

c)

> X(s):=(s^2-2*s+1)/(s^2*(s^2+4)):
x(t)=invlaplace(X(s),s,t);

d)

> X(s):=(s^2-2*s+1)/(s^2*(s^2+4)^2):
x(t)=invlaplace(X(s),s,t);

e)

> X(s):=(s^2-2*s+1)/((s+2)^2+4):
x(t)=invlaplace(X(s),s,t);

f)

> X(s):=(s^2-2*s+1)/(s*((s+2)^2+4)):
x(t)=invlaplace(X(s),s,t);

g)

> X(s):=(s^2-2*s+1)/((s+2)^2+4)^2:
x(t)=invlaplace(X(s),s,t);

h)

> X(s):=(s^2-2*s+1)/(s*((s+2)^2+4)^2):
x(t)=invlaplace(X(s),s,t);

Opgave 7 - Convolutie

Bereken, gebruik makend van de Laplace-transformatie, de convolutie (x*v)(t) voor het geval dat

> conv := proc(x::uneval,y::uneval)

> invlaplace( laplace(x(t),t,s)*laplace(y(t),t,s),s,t);

> end proc:

a)

> x:=t-> exp(-t)*Heaviside(t):
v:=t-> sin(t)*Heaviside(t):

> conv(x,v);

b)

> x:=t-> cos(t)*Heaviside(t):
v:=t-> sin(t)*Heaviside(t):

> conv(x,v);

Opgave 8 - ODE met randvoorwaarden

Een tijdcontinu systeem wordt beschreven door de ingangs/uitgangs-differentiaalvergelijking:

> unassign('x(t)','y(t)'):

> de_io:= diff(y(t),t$2)+4*diff(y(t),t)+3*y(t) =
2*diff(x(t),t$2)-4*diff(x(t),t)-x(t);
assume(t>=0);

Bereken het uitgangssignaal y (t) van dit systeem voor t>=0 voor de volgende gevallen

a)

De voorwaarden

> ic1:= y(0)=-2: # y(0-)=-2

> ic2:= D(y)(0)=1: # y'(0-)=1

> x:= t-> 0: # x(t)=0

Eerst maar eens handmatig met behulp van de laplace-transformatie.

Naar het laplace-domein

> algb_eq:=laplace(de_io,t,s);

Initial Conditions invullen

> algb_ic:=eval(algb_eq,{ic1,ic2});

Oplossen in laplace domein

> algb_sol:=solve(algb_ic,laplace(y(t),t,s));

En weer terug naar het t-domein.

> y_out:=invlaplace(algb_sol,s,t);

En dan nu in een keer door maple

dsolve({de_io,ic1,ic2},y(t), method=laplace);

b)

De beginvoorwaarden (initial conditions)

> ic1:= y(0)=0:

> ic2:= D(y)(0)=0:

> x:= t->Dirac(t):

En oplossen die hap

> dsolve({de_io,ic1,ic2},y(t),method=laplace);

c)

> ic1:= y(0)=0:

> ic2:= D(y)(0) =0:

> x:= t->Heaviside(t):

> dsolve({de_io,ic1,ic2},y(t),method=laplace);

d)

> ic1:= y(0)=-2:

> ic2:= D(y)(0)=1:

> x:= t->Heaviside(t):

> dsolve({de_io,ic1,ic2},y(t),method=laplace);

>

Serie 2 - Fourier

> restart: with(inttrans):

Opgave 9

Bereken de convolutie van x met zichzelf voor -infinity < t <infinity voor het geval dat

> assume(t,real):x(t):=sin(t)/t;

Convolutie in het t-domein is gelijk aan de vermenigvuldiging in het omega-domein.

> Y(omega):= simplify(fourier(x(t),t,omega)^2);

> y(t):= simplify(invfourier(Y(omega),omega,t));

Opgave 10 en 11

Schrijf elk van de volgende signalen als een reeks van complexe exponenten (Maak hierbij gebruik van trigonometrische identiteiten).

Bereken de (complexe) Fourier-getransformeerde van elk van de signalen. Vergelijk de uitkomsten met de reeksen.

a)

> x:=t-> cos(5*t+Theta);

Complexe reeks; lang leve maple :)

> 'x(t)'=convert(convert(convert(x(t),exp),`+`),exp);

En de Fourier-getransformeerde

> 'X(omega)'=fourier(x(t),t,omega);

b)

> x:=t-> sin(t)+cos(t);
'x(t)'=convert(convert(x(t),exp),`+`);
'X(omega)'=fourier(x(t),t,omega);

c)

> x:=t-> (sin(4*t))^2;
'x(t)'=convert(expand(convert(x(t),exp),exp),`+`);
'X(omega)'=fourier(x(t),t,omega);

d)

> x:=t-> cos(2*t)*sin(3*t);
'x(t)'=convert(combine(convert(x(t),exp),exp),`+`);
'X(omega)'=fourier(x(t),t,omega);

e)

> x:=t-> (cos(5*t))^2;
'x(t)'=convert(combine(convert(x(t),exp),exp),`+`);
'X(omega)'=fourier(x(t),t,omega);

f)

> x:=t-> cos(3*t)+cos(5*t);
'x(t)'=convert(combine(convert(x(t),exp),exp),`+`);
'X(omega)'=fourier(x(t),t,omega);

Opgave 12

Een tijdcontinu signaal x(t) heeft de Fourier-getransformeerde

> x(t):='x(t)':t:='t':X(omega):= 1/(I*omega+3);

Bereken de Fourier-getransformeerde van de volgende signalen:

Eerst maar eens x(t) zelf berekenen.

> x:=t->invfourier(X(omega),omega,t): 'x(t)'=x(t);

a)

De functie die we moeten Fourier-transformeren

> v(t):=x(5*t-4);

En de Fourier-getransformeerde

> V(omega)=simplify(fourier(v(t),t,omega));

b)

> v(t):=t^2*x(t);
V(omega)=simplify(fourier(v(t),t,omega));

c)

> v(t):=x(t)*exp(2*I*t):
V(omega)=simplify(fourier(v(t),t,omega));

d)

> v(t):=x(t)*cos(4*t):
V(omega)=simplify(fourier(v(t),t,omega));

e)

> v(t):=diff(x(t),t$2):
V(omega)=simplify(fourier(v(t),t,omega));

f)

Hier wordt de convolutie met zichzelf gevraagd en dat is gewoon het product met zichzelf in het frequentie-domein.

> V(omega)=simplify(X(omega)^2);

g)

> v(t):=(x(t))^2:
V(omega)=simplify(fourier(v(t),t,omega));

h)

> v(t):=1/(I*t-3):
V(omega)=simplify(fourier(v(t),t,omega));

Opgave 13

Bereken de Fourier-getransformeerde van de volgende signalen:

a)

> x(t):=exp(-t)*cos(4*t)*Heaviside(t):
'X(omega)'=simplify(fourier(x(t),t,omega));

b)

> x(t):=t*exp(-t)*Heaviside(t):
'X(omega)'=simplify(fourier(x(t),t,omega));

c)

> x(t):=cos(4*t)*Heaviside(t):
'X(omega)'=collect(simplify(fourier(x(t),t,omega)),Dirac);

d)

> x(t):=exp(-abs(t)):
'X(omega)'=fourier(x(t),t,omega);

e)

> x(t):=exp(-t^2):
'X(omega)'=fourier(x(t),t,omega);

Opgave 14

Een ideaal bandfilter heeft de overdrachtsfunctie

> H(omega):= piecewise( abs(omega)>2 and abs(omega)<4, 10*exp(-4*I*omega));

Bereken het uitgangssignaal y(t) van het filter voor het volgende ingangssignaal

> assume(t,real):x(t):=sin(2*t)/(2*t) *cos(t);

Het uitgangssignaal is de convolutie van het ingangssignaal met de overdrachtsfunctie en in het Fourier-domein is dat gewoon vermenigvuldiging.

> X(omega):=simplify(fourier(x(t),t,omega));

We moeten de overdrachtsfunctie naar een Heaviside-vorm converteren, omdat invfourier er anders niet mee om kan gaan.

> Y(omega):= simplify(X(omega) * convert(H(omega),Heaviside));

En nu weer terug naar het t-domein

> y(t):=simplify(invfourier(Y(omega),omega,t));

> plot([y(t),x(t)],t=-10..10,legend=["y(t)","x(t)"]);

Opgave 15

Een ideaal laagdoorlaat-filter heeft de overdrachtsfunctie

> H(omega):= piecewise(-2<omega and omega<2,exp(-I*omega));

Bereken het uitgangssignaal y(t) van het filter voor de volgende ingangssignalen x(t) :

a)

> assume(t,real):x(t):= 5*sin(t/2)/(t/2) *cos(2*t);

Dit is weer convolutie van de twee en dat is vermenigvuldiging in het frequentie-domein.

> Y(omega):= simplify(convert(H(omega),Heaviside) * fourier(x(t),t,omega));

> y(t):= simplify(invfourier(Y(omega),omega,t));

> plot([y(t),x(t)],t=-10..10,legend=["y(t)","x(t)"]);

b)

> assume(t,real): x(t):=Sum(cos(n*Pi*t/2+Pi/6)/n,n=1..infinity);

SORRY, DEZE UITWERKING IS NIET CORRECT NOCH COMPLEET !

> X(omega):=fourier(x(t),t,omega);

Als je hiermee met maple verder gaat, kom je niet op iets bruikbaars. Maar we weten dat het een laagdoorlaatfilter betreft waarbij de frequentie tussen -2 en 2 zit. Aangezien overal die Dirac-functie in zit, blijft na vermenigvuldiging met de overdrachtsfunctie slechts n=1 over. Dit doen we eventjes handmatig:

> Y(omega):= fourier(cos(n*Pi*t/2+Pi/6),t,omega) * exp(-I*omega);

En dan weer terug naar het t-domein:

>

>

>

Serie 3 - Z-transformaties

> restart; with(inttrans): readlib(ztrans):
assume(n,integer); # We hebben hier te maken met discrete signalen

Ter beschikking staan met ztrans:

* ztrans(xn, n, z)

* invztrans(Xz,z,n)

# Hiervoor moet je (complex vaak) soms breuksplitsen met convert(Xz,parfrac,z)=breuksplitsen

* charfcn[i](n) = Delta(n-i)

Simplify ideeen:

- convert(xn,radical)

>

Opgave 16 - Z-transformeren

Een tijddiscreet signaal x(kT) heeft de z -getransformeerde

> X:= z-> z/(8*z^2-2*z-1);

Bereken de z -getransformeerde V(z) van de volgende signalen.

Goed, eerst maar eens x(kT) uitrekenen.

> x:= a-> subs(q=a,invztrans(X(z),z,q)): 'x(a)'=x(a);

a)

in v(kT) worden discrete stapjes kT genomen. kT+2*T = T*(k+2) en T constant

Het is een beetje lastig met die kT als argument, maar we nemen gewoon n in de plaats daarvan. Ik snap het ook niet helemaal precies, maar dat is mij verteld...

> v(n):=x(n+2)*Heaviside(n+2);

En dan de functie Z-transformeren

> V(z)=simplify(ztrans(v(n),n,z));

b)

> v(n):=exp(3*n)*x(n);
V(z)=simplify(ztrans(v(n),n,z));

c)

> v(n):=cos(2*n)*x(n):
V(z)=simplify(ztrans(v(n),n,z));

d)

> v(n):=n^2*x(n);
V(z)=simplify(ztrans(v(n),n,z));

Opgave 17 - Z-terugtransformeren

Bereken de inverse z-getransformeerde x[k] van de volgende signalen voor alle integers .

a)

> X:= z-> z/(z^2+1);

Alle verschillende inverse z-transformatie uitkomsten berekenen

> pre_x:= allvalues(invztrans(X(z),z,n));

Goed, en dan nu de uiteindelijke complex uitgewerkte functie is de som van de uitkomsten:

> x(n):= simplify(evalc(pre_x[1]+pre_x[2]));

b)

> X:=z-> z^2/(z^2+1);

> pre_x:= allvalues(invztrans(X(z),z,n));
x(n):= simplify(evalc(pre_x[1]+pre_x[2]));

c)

> X:=z-> 1/(z^2+1) + 1/(z^2-1);
pre_x:= invztrans(X(z),z,n);

Hier moeten we even opletten. Zoals je ziet is pre_x opgebouwd uit een sommatie-deel met een RootOf-term (die voor twee uitkomsten zorgt), en een ander deel. Als we hierop allvalues gaan toepassen en de termen bij elkaar optellen, wordt het constante deel er twee keer bij opgeteld. Dus we moeten dat constante deel even afsplitsen, allvalues toepassen en het er aan het einde weer bij optellen.

> pre_x_ro:= select(has,pre_x,RootOf); # het deel met RootOf
pre_x_c:= remove(has,pre_x,RootOf); # het constante deel

Goed, nu allvalues toepassen op pre_x_ro

> pre_x_ro_allv:= allvalues(pre_x_ro);

En nu alles bij elkaar optellen

> x(n):=simplify(evalc(pre_x_ro_allv[1]+pre_x_ro_allv[2]) + pre_x_c);

d)

> X:=z-> (z^2-1)/(z^2+1);
pre_x:= invztrans(X(z),z,n):
pre_x_ro_allv:= allvalues(select(has,pre_x,RootOf));
pre_x_c:= remove(has,pre_x,RootOf);
x(n):=simplify(evalc(pre_x_ro_allv[1]+pre_x_ro_allv[2]) + pre_x_c);

Even een opmerking over charfcn . charcfn[k](n) is eigenlijk de discrete equivalent van Dirac(n-k) . Dus deze functie is 1 als n=k en anders 0.

e)

> X:=z-> (z+2)/( (z-1)*(z^2+1) );
pre_x:= invztrans(X(z),z,n):
pre_x_ro_allv:= allvalues( select(has,pre_x,RootOf) );
pre_x_c:= remove(has,pre_x,RootOf);
x(n):=simplify(evalc(pre_x_ro_allv[1]+pre_x_ro_allv[2])+pre_x_c);

f)

> X:=z-> ln( (2*z-1) / (2*z) );
pre_x:= invztrans(X(z),z,n);

Hee, dit werkt niet :( Dan moeten we nog iets anders proberen wat maple niet weet, en dat is het volgende:

<=>

Dus we kunnen X differentieren en met -z vermenigvuldigen en misschien kan maple er dan wel mee overweg. Aan het einde hoeven we dan alleen het antwoord door n te delen.

> X2:= simplify(-z*diff(X(z),z));

Dit ziet er leuk uit :)

> x2:= invztrans(X2,z,n);

Nu nog corrigeren voor wat we in het begin gedaan hebben:

> x:= simplify( x2/n );

Ehm, dit komt niet overeen met het antwoord op het antwoordblad; ik moet iets over het hoofd gezien hebben...

Opgave 18 - Z-terugtransformeren

Bereken de inverse z-getransformeerde x[n] van

de volgende signalen voor alle integers op de reele as.

a)

> X(z):= (z+0.3) / (z^2+0.75*z+0.125);
'x(n)' = invztrans(X(z),z,n);

b)

> X(z):= (5*z+1) / (4*z^2+4*z+1);
'x(n)' = invztrans(X(z),z,n);

c)

Probeer de volgende regel maar eens als je zin hebt; maar het gaat erg lang duren; 'too many levels of recursion'.

> X(z):= (4*z+1) / (z^2 - z + 0.5);
invztrans(X(z),z,n);

Warning, computation interrupted

Dus we gaan het maar eens symbolisch proberen ipv numeriek (wat maple doet als 'ie een float, kommagetal tegenkomt).

> X(z):= (4*z+1) / (z^2-z+1/2);
pre_x:= invztrans(X(z),z,n);
pre_x_ro_allv:= allvalues(select(has,pre_x,RootOf));
pre_x_c:= remove(has,pre_x,RootOf);
x(n):=simplify(evalc(pre_x_ro_allv[1]+pre_x_ro_allv[2]) + pre_x_c);

d)

> X(z):= z/(16*z^2+1);
pre_x:= allvalues(invztrans(X(z),z,n));
x(n):= simplify(evalc(pre_x[1]+pre_x[2]));

e)

> X(z):= (2*z+1) / (z*(10*z^2-z-2));
x(n):= invztrans(X(z),z,n);

f)

Zo te zien geeft deze ook problemen met de numerieke methode.

> X(z):= (z+1) / ( (z-0.5)*(z^2-0.5*z+0.25) );
invztrans(X(z),z,n);

Warning, computation interrupted

Dus naar symbolische methode met breuken.

> X(z):= (z+1) / ( (z-1/2)*(z^2-z/2+1/4) );
pre_x:= invztrans(X(z),z,n);
pre_x_ro_allv:= allvalues(select(has,pre_x,RootOf));
pre_x_c:= remove(has,pre_x,RootOf);
x(n):=simplify(evalc(pre_x_ro_allv[1]+pre_x_ro_allv[2]) + pre_x_c);

g)

Als je het probeert, zie je dat ook deze numeriek niet echt wil.

> X(z):= (z+1)/( z^2*(z^2-z/2+1/4) );
invztrans(X(z),z,n);

En je ziet dat je rekening moet houden met 'constante' factoren dus moeten we even splitsen enzo.

> pre_x_ro_allv:= allvalues(select(has,pre_x,RootOf));
pre_x_c:= remove(has,pre_x,RootOf);
x(n):=simplify(evalc(pre_x_ro_allv[1]+pre_x_ro_allv[2]) + pre_x_c);

Opgave 19 - Bereken impulsresponsie met de overdrachtsfunctie

Bereken de impulsresponsie h(kT) van het lineaire tijdinvariante tijddiscrete systeem met de overdrachtsfunctie

> H(z):= z / ( (z-0.5)^2*(z^2+0.25) );

Eerst gewoon maar proberen

> invztrans(H(z),z,kT);

Warning, computation interrupted

Nou, dit duurt ook weer te lang (of too many levels of recursion). Dan maar mbv breuken.

> H(z):= z/ ( (z-1/2)^2 * (z^2+1/4) );
pre_h:=invztrans(H(z),z,n);

Goed, even rekening houden met 'constante' factoren en het antwoord komt eruit:

> pre_h_ro_allv:= allvalues(select(has,pre_h,RootOf));
pre_h_c:= remove(has,pre_h,RootOf);
h(n):=simplify(evalc(pre_h_ro_allv[1]+pre_h_ro_allv[2]) + pre_h_c);

This file was created by Willem van Engen.