Should ztrans Work Better when Using heaviside() ?

Recently ran across peculiar behavior with ztrans.

Define some variables

syms n integer
syms T positive

Set sympref so that heaviside[n] is the discrete-time unit step

sympref('HeavisideAtOrigin',1);

Simple function to take the z-tranform of a signal and then the inverse z-transform, which should return the input signal

g = @(f) [ztrans(f) ; simplify(iztrans(ztrans(f)))];

1. Simple case, works as expected

g([cos(n) sin(n)])

ans = 

2. Introduce symbolic sampling period, still works fine.

g([cos(T*n) sin(T*n)])

ans = 

3. Multiply the simple case by heaviside[n], which shouldn't have any material effect because heaviside[0] = 1

g([cos(n) sin(n)]*heaviside(n))

ans = 

It has no material effect, but we can see that ztrans is accounting for the possibility that heaviside[0] might not be 1, so it appears to be using the general rule

y[n] = f[n]*heaviside[n] = f[0]*heaviside[0]*kroneckerDelta[n] + f[n]*u[n-1]

where u[n] is the discrete-time unit step. The z-transform is then

Y(z) = f[0]*heaviside[0] + Z(f[n+1])/z

As an aside, it seems like it would be easier to solve by

y[n] = f[n]*heaviside[n] = f[n] + (f[0]*heaviside[0] - f[0])*kroneckerDelta(n)

Y(z) = Z(f[n]) - f[0]*h[0]

4. Now introduce the sampling period and heaviside. Based on 2 and 3 above, it would appear that the ztrans should be easily found as would the iztrans.

g([cos(T*n) sin(T*n)]*heaviside(n))

ans = 

But, in neither case is the closed form expression returned for the z-transform. Oddly, iztrans returns the closed form expression for cos, but not for sin.

[David Goodmanson]

Hi Paul,

can you comment on

heaviside(0)
  ans =  0.5000
heaviside(sym(0))
  ans = 1/2

[Paul]

Greetings David,

Not sure what you're asking for. After this command

sympref('HeavisideAtOrigin',1);

we have

heaviside(0)

ans = 1

heaviside(sym(0))

ans = 

1

[David Goodmanson]

Yes, my mistake, I spaced out on that.

There are some odd things going on. I tried out complex exponentials and for case 4, the T*n heaviside case, they work. So you could get the correct result with linear combinations of those, but that is not really a satisfactory situation.

I used an arbitrary value p for heaviside(0) in order to show the effect of that. With the output transposed from yours,

syms n integer
syms T positive
syms p
sympref('HeavisideAtOrigin',p);
f = @(q) [cos(q); sin(q); exp(i*q); exp(-i*q)];
fh = @(q) [cos(q); sin(q); exp(i*q); exp(-i*q)]*heaviside(n);
g = @(f) [ztrans(f) simplify(iztrans(ztrans(f)))];
g(f(n))     % case 1
g(f(T*n))   % case 2
g(fh(n))    % case 3
g(fh(T*n))  % case 4
sympref('HeavisideAtOrigin',1/2)
% since the sympref setting is remembered after a session
% (I found that out the hard way).

Case 3 shows the effect of heaviside at the orgin.

 
% case 1
[(z*(z - cos(1)))/(z^2 - 2*cos(1)*z + 1),     cos(n)]
[      (z*sin(1))/(z^2 - 2*cos(1)*z + 1),     sin(n)]
[                        z/(z - exp(1i)),  exp(n*1i)]
[                       z/(z - exp(-1i)), exp(-n*1i)]

% case 2
[(z*(z - cos(T)))/(z^2 - 2*cos(T)*z + 1),     cos(T*n)]
[      (z*sin(T))/(z^2 - 2*cos(T)*z + 1),     sin(T*n)]
[                      z/(z - exp(T*1i)),  exp(T*n*1i)]
[                     z/(z - exp(-T*1i)), exp(-T*n*1i)]

% case 3
[p - (z - (z^2*(z - cos(1)))/(z^2 - 2*cos(1)*z + 1))/z,     cos(n) + p*kroneckerDelta(n, 0) - kroneckerDelta(n, 0)]
[                    (z*sin(1))/(z^2 - 2*cos(1)*z + 1),                                                     sin(n)]
[                               p + 1/(z*exp(-1i) - 1),  exp(n*1i) + p*kroneckerDelta(n, 0) - kroneckerDelta(n, 0)]
[                                p + 1/(z*exp(1i) - 1), exp(-n*1i) + p*kroneckerDelta(n, 0) - kroneckerDelta(n, 0)]
 

Case 4, who knows why sin and cos are doing what they do.

% case 4
[p + ztrans(cos(T*(n + 1)), n, z)/z,     cos(T*n) + p*kroneckerDelta(n, 0) - kroneckerDelta(n, 0)]
[    ztrans(sin(T*(n + 1)), n, z)/z,                iztrans(ztrans(sin(T*(n + 1)), n, z)/z, z, n)]
[          p + 1/(z*exp(-T*1i) - 1),  exp(T*n*1i) + p*kroneckerDelta(n, 0) - kroneckerDelta(n, 0)]
[           p + 1/(z*exp(T*1i) - 1), exp(-T*n*1i) + p*kroneckerDelta(n, 0) - kroneckerDelta(n, 0)]

[Paul]

Hi David,

% since the sympref setting is remembered after a session
% (I found that out the hard way).

I too am a member of that club. I put

sympref('default')

into my startup.m script and then modify as needed.

Your results with the complex exponentials make the (lack of) results with sin and cos even more puzzling.