%
% Routine for calculating a finite length Z-transform in a straight forward way
%--------------------------------------------------------------------------
% Parameters: start
%--------------------------------------------------------------------------
%
% - Choose a value for z (any complex number):
%
z = -0.6 - 1i*0.7;
% -- Some other suggested options for z: uncomment to try:
%z = 0.99*1i;
%z = 1.01*exp(-1i*pi/10);
%z = 0.5;
%z = -1;
%
% - Define a sequence and its corresponding indices:
%
% -- First the indices (any set of unique integers will do):
n = [0:100];
% --- Some other choices for n: uncomment to try:
%n = [0:100] - 50;
%n = [300 -7 12 0 5]; % make sure to change x to only be length 5!
%
% -- Now anything at all for x, as long at it's the same length as n
x = (0.5 + 1i*0.6).^n;
% --- Some other choices for x: uncomment to try:
%x = (-0.95).^n; % uncomment to try a different choice for x
%x = 1./(abs(n)+1);
%x = atan(n);
%x = [1 2 3 4 5]; % note there are only 5 entries: n must have same size.
%
%--------------------------------------------------------------------------
% Parameters: end
%--------------------------------------------------------------------------
% First check to be sure that the x sequence is as long as it's index
% sequence n:
if length(x) ~= length(n)
error('Length of x does not match the length of n');
end
% Let's take a look at what the x sequence looks like
% - find how to scale the y-axis:
xmax = ceil(10*max(abs([real(x) imag(x)])))/10;
% - now for the plot:
figure(1),
subplot(4,1,1)
stem(n,real(x),'r.');
title('real(x[n])');
ylim([-1 1]*xmax);
subplot(4,1,2)
stem(n,imag(x),'b.');
title('imag(x[n])');
ylim([-1 1]*xmax);
% Compute the sequence of z terms and put them in a vector called zvec:
zvec = z.^(-n);
% Plot zvec too
% - find how to scale the y-axis:
zmax = ceil(10*max(abs([real(zvec) imag(zvec)])))/10;
% - now for the plot:
subplot(4,1,3)
stem(n,real(zvec),'r.');
title(['real(z^{-n}) when z = ' num2str(z)]);
ylim([-1 1]*zmax);
subplot(4,1,4)
stem(n,imag(zvec),'b.');
title(['imag(z^{-n}) when z = ' num2str(z)]);
ylim([-1 1]*zmax);
xlabel('n');
% Calculate all the terms of the sum:
terms = x.*zvec;
% Plot the terms as well:
% Plot zvec too
% - find how to scale the y-axis:
tmax = ceil(10*max(abs([real(terms) imag(terms)])))/10;
% - now for the plot:
figure(2),
subplot(2,1,1)
stem(n,real(terms),'r.');
title(['Real part of terms of X(z): real(x[n] z^{-n}) when z = ' num2str(z)]);
ylim([-1 1]*tmax);
subplot(2,1,2)
stem(n,imag(terms),'b.');
title(['Imaginary part of terms of X(z): imag(x[n] z^{-n}) when z = ' num2str(z)]);
ylim([-1 1]*tmax);
xlabel('n');
% Sum the terms to calculate the Z-transform:
Xz = sum(terms);
display(['Z-transform of x = Z{x[n]} = X(z) and X(z=' num2str(z) ') = ' num2str(Xz)]);
% If by chance the user selected a geometric sequence for x[n] we can
% calculate X(z) much easier like this:
if length(x) > 1
a = x(2)/x(1);
Xz_for_geo_x = (terms(1) - terms(end)*a*(z^-1))/(1 - a*(z^-1));
display(['If x[n] was a geometric sequence then we can guess that X(' num2str(z) ') = ' num2str(Xz)]);
d = abs(Xz - Xz_for_geo_x);
if d <= abs(Xz)*1e-10
display(['A relative error of ' num2str(d) ' indicates x[n] was probably a geometric sequence!']);
else
display(['A relative error of ' num2str(d) ' indicates x[n] was NOT a geometric sequence!']);
end
end
% END script
Here's some example output:
Z-transform of x = Z{x[n]} = X(z) and X(z=-0.6-0.7i) = 0.54138-0.0034482i
If x[n] was a geometric sequence then we can guess that X(-0.6-0.7i) = 0.54138-0.0034482i
A relative error of 1.2143e-17 indicates x[n] was probably a geometric sequence!
 |
| Figure 1 |
 |
| Figure 2 |
Here's something else you can try: right after the parameters block add a line to the script like this:
x[1] = 10*x[1];
Now unless x[n] = 0 for all n, you no longer have a geometric sequence and the "quick and dirty" method of calculation will fail (which the script should detect). Other things you can do include replacing the quick and dirty method with some of those expressions from Wikipedia's Z-transform table, and then replacing x with the sequence from the table as well (a version of x with a finite number of terms of course), and then see if the quick and dirty method comes up with the same answer as does summing each and every term in the series.
Note that what I call the quick and dirty method actually looks more complicated in Octave or Matlab script, but behind the scenes many many more calculations are being done in the "straightforward method" (i.e. computing all the terms of the sum, and summing them) for any sequence x which is non-zero for any appreciable length.
Also, this line which implements the "quick and dirty" method for geometric sequences x:
Xz_for_geo_x = (terms(1) - terms(end)*a*(z^-1))/(1 - a*(z^-1));
came directly from this comment I left for you, especially the equation at the bottom and the one immediately above it.
Some things you might try:
- Change figure 2 so that you also plot the cumulative sum (both the real and imaginary components) of the terms of X(z) starting from the term at the smallest index.
- Add this asymptote to these cumulative sum curves as a dashed horizontal line (if there is an asymptote).
- For all the curves add an +/- envelope with a dashed line (the envelope is just the magnitude and the negative of the magnitude of the complex numbers)
- For geometric series, add a calculation of the Z-transform if the series were infinite in length rather than finite, and then show the difference between the Z-transform of the finite series and its infinite cousin (which will often be very very small, if z was selected to be inside the ROC.). Hint: the formula will be similar to the formula I link to above, and only a small modification of the line of script I copy right above that link.
- Add plots of the phase and magnitude of the complex numbers (angle() and abs() give those) rather than (or in addition to) the real and imaginary components.
Wow. This is cool. I will come back to it soon.
ReplyDeleteI added some suggestions to try, and I edited the parameters block to add some suggested parameter settings (currently commented out). I also edited one of the display() statements, and the example output to match. I hope I didn't break the script in the process!
DeleteNotice that the actual calculation of the Z-transform is very simple and can be done in one line:
Xz = sum(x.*(z.^(-n)));
... or even more compactly, as a dot product (AKA inner product) of two vectors:
DeleteXz = (z.^(-n))*(x.');
Where x' is the Hermitian (conjugate) transpose of x and x.' means to skip the conjugation and only do the transpose. Thus x' = conj(x.') = conj(x).'
x' and x.' are equivalent for real valued vectors, so for real valued vectors and matrices, most people skip the ".". There's a very good reason why a conjugate is normally done as well, when you transpose complex arrays, but I won't get into that here: suffice it to say, we do not want the conjugation in this case.
Also note that in addition to suppressing the default conjugation operation of the transpose operator ' the period operator "." also generally suppresses the normal rules for multiplying matrices or vectors and matrices or vectors and vectors, and a few other related things. It's meaning is "do the following operation element by element." Thus if A and B are equal dimension matrices or vectors (or arrays of any kind: 3D, 4D, etc), then A.*B does not mean a matrix multiply, it means form a new matrix (call it C and is C = A.*B) where each element of C is just a scalar product of the corresponding elements of A and B. For example if the elements of A are aij and the elements of B are bij, then the elements of C are cij = aij*bij.
The other related operation of the period "." operator is to do something like this: suppose A is a matrix, then 2^A makes no sense, however, B = 2.^A means to form a new matrix B, the same size as A, but with each element bij = 2^aij. It works the other way around as well, so B = A.^2 means bij = aij^2.
And of course I've already covered the special case of the natural log base: exp(A) means element by element, so exp(A) = e.^A. But expm(A) is the matrix exponential, which is a different beast. There's also a corresponding logm() with is the matrix logarithm, and it's the functional inverse of expm(). Likewise, it does not mean element by element logarithm of A (I think log(A) would do that).
It turns out you don't need all those parentheses: this works too:
DeleteXz = z.^-n*x.'
... and remember, for a given complex sequence (or vector) x, the Z-transform is just a mapping from scalars in the complex plane to other scalars in the complex plane, so the above line only produces a single complex number as a result, thus you can safely leave the semicolon off the end, and it will just print like that on your console. Also, I hope you realize that Octave or Matlab scripts are just a way to play back command line instructions, so once you run one all those variables are still available for you to examine with command line instructions, such as:
DeleteXz
which will print the result on the console, even after the script completes.
I was trying to match Wikipedia's notation for the Z-transform, which is something like this:
Z-transform of sequence x = Z{x[n]} = X(z)
However, I don't think that's the best notation, because x[n] implies the just the nth element of the sequence, and it fails to account for the z argument. So I prefer something like this:
Z-transform-of-sequence-x = Z(x,z) = X(z)
That lets you know exactly what the transform is a function of: the whole sequence x, and the complex scalar z.
Once you understand a Z-transform, then understanding a discrete time Fourier transform is a breeze, and from there, understanding a continuous time Fourier transform is pretty easy, and then a Laplace transform follows from that. You can think of all four of them as transformations from the time domain to the frequency domain (where frequencies are expressed as complex numbers), and they all have an inverse as well. And then you'll realize that there's really only two kinds: Z and Laplace (the two Fourier transforms are just special cases of those), ... and of course their inverses.
I guess Z{x,z} or Z{x;z} works too.
DeleteThis is good too, but you need MathJax for it: xZ→X
DeleteThen for Fourier and Laplace you'd have:
DeletexF→X (Fourier)
xL→X (Laplace)
XF−1⟶x (inverse Fourier)
xF⟺X (are a Fourier transform pair)
...
etc.