Z-transform of a downsampler

Question

In this paper or multirate filtering, the author establishes the following mathematical relationship. Let $y_D$ be the output of a downsampler such that

$$y_D[n] = x[Mn]$$

where $M$ is the downsampling factor. In other words, we keep every $M$-th sample of the original signal. The author then goes on to state the following:

...the z-transform of $y_D[n]$ is given by

$$Y_D[z]=\frac{1}{M}\sum_{k=0}^{M-1}X[z^{1/M}W^k]$$

where $W^k$ is the $M$-point Discrete Fourier Transform kernel, namely $e^{(-j2\pi k)/M}$.

How can we go from the former expression to the latter? What is the relationship between DFT and the Z-transform that allows for such transition?

Answer 1

This derivation is a tricky one. The approach suggested before has a flaw. Let me demonstrate this first; then I will give the correct solution.

We wish to relate the $\mathcal{Z}$-transform of the downsampled signal, $Y_D(z) = \mathcal{Z}\{x[Mn]\}$, to the $\mathcal{Z}$-transform of the original signal $X(z) = \mathcal{Z}\{x[n]\}$.

The wrong way

One could think of simply plugging the expression for the downsampled signal into the expression of the $\mathcal{Z}$-transform:

$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{x[Mn] z^{-n}}$

A change of variable $n'=Mn$ seems obvious:

$Y_D(z) = \sum\limits_{n' \in M\mathbb{Z}}{x[n'] z^{-n'/M}}$

However, it is important to realize that even though the new summation index $n'$ still runs from $-\infty$ to $\infty$, the sum is now over 1 out of M integer numbers. In other words,

$n' \in M\mathbb{Z}=\{...,-2M,-M,0,M,2M,...\}$,

while the definition of the $\mathcal{Z}$-transform requires

$n \in \{...,-2,-1,0,1,2,...\}$.

Since this is no longer a $\mathcal{Z}$-transform, we cannot write:

$Y_D(z) = X(z^{1/M})$

The right way

Let us first define a 'helper' impulse train signal $t_M[n]$ as:

$\begin{align*} t_M[n] & = \sum\limits_{k=-\infty}^{+\infty}{\delta[n-kM]} \\ & = \left\{ \begin{array}{lr} 1 & : n \in M\mathbb{Z}\\ 0 & : n \notin M\mathbb{Z} \end{array} \right. \end{align*} $

This function is $1$ at one out of every $M$ samples, and zero everywhere else.

Equivalently, the pulse train function can be written as:

$ t_M[n] = \frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}} $

Proof: We need to consider separately the cases $n \in M\mathbb{Z}$ and $n \notin M\mathbb{Z}$:

$$\begin{align*} t_M[n] & = \frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}} \\ & = \left\{ \begin{array}{lr} \frac{1}{M} \sum\limits_{k=0}^{M-1}{1} && : n \in M\mathbb{Z}\\ \frac{1}{M} \frac{1-e^{j2\pi n}}{1-e^{j2\pi n / M}} && : n \notin M\mathbb{Z} \end{array} \right. \\ & = \left\{ \begin{array}{lr} \frac{1}{M} M && : n \in M\mathbb{Z}\\ \frac{1}{M} \frac{1-1}{1-e^{j2\pi n / M}} && : n \notin M\mathbb{Z} \end{array} \right. \\ & = \left\{ \begin{array}{lr} 1 && : n \in M\mathbb{Z}\\ 0 && : n \notin M\mathbb{Z} \end{array} \right. \end{align*}$$ In the case $n \notin M\mathbb{Z}$, we used the expression for the finite sum of a geometric series.

Now let's go back to our original problem of finding the $\mathcal{Z}$-transform of a downsampler:

$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{x[Mn] z^{-n}}$

We apply the substitution $n'=Mn$, keeping in mind that this makes the summation run only over integer multiples of M:

$Y_D(z) = \sum\limits_{n' \in M\mathbb{Z}}{x[n'] z^{-n'/M}}$

We can now use the above impulse train function to safely rewrite this as a summation over all $n \in \mathbb{Z}$:

$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{t_M[n] x[n] z^{-n/M}}$

Using the above formulation for the impulse train function as a finite sum of exponentials, we get:

$\begin{align*} Y_D(z) &= \sum\limits_{n=-\infty}^{+\infty}{(\frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}}) x[n] z^{-n/M}} \\ &= \frac{1}{M} \sum\limits_{k=0}^{M-1}{\sum\limits_{n=-\infty}^{+\infty}{e^{j2\pi k n / M} x[n] z^{-n/M}}} \\ &= \frac{1}{M} \sum\limits_{k=0}^{M-1}{\sum\limits_{n=-\infty}^{+\infty}{x[n] (e^{-j2\pi k / M} z^{1/M})^{-n}}} \end{align*}$

The summation on the right is a summation over all integers, and is therefore a valid $\mathcal{Z}$-transform in terms of $z'=e^{-j2\pi k / M} z^{1/M}$. Therefore, we can write:

$ Y_D(z) = \frac{1}{M} \sum\limits_{k=0}^{M-1}{X(e^{-j2\pi k / M} z^{1/M})} $

This is the formula for the $\mathcal{Z}$-transform of a downsampler.

Answer 2

I've not seen this notation before. However, it does seem to make sense. The $M$-downsampler is defined by the equation:

$$ y_D[n] = x[Mn] $$

Its $z$ transform is defined by the equation:

$$ \begin{align} Y_D(z) &= \sum_{n=-\infty}^\infty y_D[n]z^{-n} \\ &= \sum_{n=-\infty}^\infty x[Mn] z^{-n} \end{align} $$

Apply a change of variable, letting $n' = Mn$. The ranges of the summation are unaffected by the change of variable since they extend to infinity.

$$ Y_D(z) = \sum_{n'=-\infty}^\infty x[n']z^{-n'/M} $$

This looks similar to the $z$ transform of $x[n]$ itself. Recall that it is defined as:

$$ X(z) = \sum_{n=-\infty}^\infty x[n]z^{-n} $$

By inspection, we can therefore conclude the following relationship between the $z$ transforms of $x[n]$ and $y_D[n]$:

$$ Y_D(z) = X\left(z^{1/M}\right) $$

Therefore, the $z$ transform of the downsampler output is closely related to the $z$ transform of the input signal, which is to be expected. In the frequency domain, this results in an $M$-fold stretching of the signal's frequency content.

But how do you go from the above equation to the one you referenced in the paper? It gives a definition of $Y_D(z)$ in terms of $z$ only, while the expression we derived is a function of $z^{1/M}$. So for a particular value of $z$ that you would like to evaluate $Y_D(z)$ at, you would first calculate $z^{1/M}$ (i.e. take the $M$-th root of $z$) and then substitute that into $X(z)$. However, all nonzero $z \in \mathbb{C}$ have $M$ distinct $M$-th roots:

$$ \left\{ r_p,\ r_p e^{\frac{j2\pi}{M}},\ r_p e^{\frac{j2\pi2}{M}},\ \ldots\ ,\ r_p e^{\frac{j2\pi(M-1)}{M}} \right\} $$

$$ = \left\{ r_p,\ r_p W,\ r_p W^{2},\ \ldots\ ,\ r_p W^{M-1} \right\} $$

where $W_k$ is the DFT kernel value $e^{j2\pi k / M}$ referenced in your question, and $r_p$ is what I define to be the principal $M$-th root of the complex value $z$:

$$ r_p = \sqrt[M]{|z|} e^{\frac{j\angle{z}}{M}} $$

That is, $z$'s principal $M$-th root $r_p$ is obtained by converting $z$ to polar form, taking the $M$-th root of $z$'s magnitude (which is a real number), and dividing $z$'s angle by $M$. The resulting values express $r_p$ in polar form.

Why go to all of this trouble? Because, as I noted before, the mapping from $Y_D(z)$'s domain to the domain of $X(z^{1/M})$ is not one-to-one. I'll now begin some handwaving. For any particular value of $z$ that you would like to evaluate $Y_D(z)$ for, there are $M$ corresponding points in $X(z^{1/M})$ that you could map to. Therefore, each of those $M$ points in $X(z^{1/M})$ contribute to the corresponding value of $Y_D(z)$. You then end up with a sum like that shown in the paper:

$$ Y_D(z) = \frac{1}{M} \sum_{k=0}^{M-1} X(r_p(z)W^k) $$

where $r_p(z)$ refers to the principal $M$-th root calculation I showed earlier. In reality, you could pick any of $z$'s $M$-th roots as the principal one; I picked this definition because it's the most straightforward. If you were to properly and rigorously derive this relationship, I believe the factor of $\frac{1}{M}$ comes in because of a derivative of $z^{1/M}$.

In mathematician-speak, I believe this would be referred to as a composition of functions; $Y_D(z) = f(g(z))$, where $f(z) = X(z)$ and $g(z) = z^{1/M}$. In order to unroll the function composition and write $Y_D(z)$ as a function of $z$ only, you would chop the domain of $Y_D(z)$ into chunks that are one-to-one, invert the function over those intervals, and then sum the results with appropriate scaling factors. I've used this technique before to calculate the probability distribution function of a function of a random variable given the original random variable's pdf (e.g. to derive the pdf of $\sqrt{X}$ given $X$'s pdf), but the name of the technique escapes me.