Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It only takes a minute to sign up.
Sign up to join this community
This question stems from the discussion we had on a previous question of mine.
The point of contention is whether the downsampling operator is time invariant or not. Which gives us a new condition to check whether such systems are LTI independent of the input is bandlimited or not. Essentially, if there's a criteria other than bandwidth which can allow for the us to check if multirate systems are LTI or not.
Directly, no - but in an important way, yes: downsampling is LTI within recovery, or "rLTI" for short. Downsampling is always linear - "summed inputs $\Leftrightarrow$ summed outputs"; time-invariance is "delayed input $\Leftrightarrow$ delayed output" : $x(t - t_0) \Leftrightarrow y(t - t_0)$. To be precise, in this post I'll use "subsampling" to refer to downsampling that only skips samples, no other steps (lowpass, etc). Post isn't basic reading, but "Basic answer" is further below.
Bandlimited subsampling is rLTI because upsampling of the subsampling is LTI. If $y[n]$ results from subsampling of $x[n]$ that's then upsampled, and $z[n]$ results from subsampling of $x[n - n_0]$ that's then upsampled, then $z[n]$ will indeed be $y[n]$ shifted by $n_0$ - i.e. $z[n] = y[n - n_0]$.
To avoid trivial-sounding statements, consider OP's original context - convolution. Let's amend the question to "Is downsampled convolution LTI for bandlimited inputs and filters?", and let $y[n]$ be the non-subsampled output of convolution of $x[n]$ with $h[n]$ What above is really saying is, if we subsample some $x[n - n_0]$, then convolve with (subsampled) $h[n]$, and both are bandlimited, then the upsampling of this result correctly recovers the original samples between the missing samples (within constant scaling, explained below) - or samples as they would be if we never subsampled.
I say "trivial-sounding" because it really is trivial: the perfect upsampling of subsampling of $x[n - n_0]$ is $x[n - n_0]$. Bandlimit is what enables this perfection. Put another way, in general, upsampling of subsampling of $x[n - n_0]$ won't recover $x[n - n_0]$ - of course it won't! Though, we're not quite done - note my amended question reads, "bandlimited inputs and filters". What if only one is? It's true that "conv in time $\Leftrightarrow$ mult in freq", and that in frequency, convolution expressed with two variables, $X\cdot H$, is same as a general signal expressed with one variable, $Z$ - and if either $X(\omega)$ or $H(\omega)$ is zero past some $\omega_\text{bw}$, then so is $Z(\omega)$. So, all good?
The catch is, it's not multiplication of spectra of $x$ and $h$, but of their subsamplings: we first subsample and then convolve. For this, we use Subsampling in time $\Leftrightarrow$ Folding in Fourier to realize, that if $x$ is bandlimited but $h$ isn't, then the spectrum of subsampled $x$ is identical but shorter (just missing zeroes), while that of $h$ changes. Same for vice versa.
Are non-bandlimited inputs ever rLTI? Yes, plenty. Two categories:
rLTI clarification: rLTI does not mean "upsampling of shifted subsampling matches shift of original" - that's a stronger result, let's call it "roLTI" for "LTI within recovery of original". Bandlimited subsampling is roLTI, and we can have rLTI without roLTI but not vice versa - full comparison below in "LTI vs other LTI".
The generalized insight is, without the caveat in point 2 (rather, accounting for it), operators are rLTI if their invertibility is unaffected by domain translations of input. This presumes an encoder-decoder pair, or simply, the existence of an (not necessarily single) inverse operation. Below sections clarify these points.
Subsampling of
$$ a = [0, 1, 2, 3, 4, 5, 6, 7], $$
by e.g. 2,
$$ a_{0, 2} = [0, 2, 4, 6] $$
isn't time-invariant, because $a$ shifted by -1,
$$ a_{-1} = [1, 2, 3, 4, 5, 6, 7, 0], $$
subsampled,
$$ a_{-1, 2} = [1, 3, 5, 7;], $$
isn't the un-shifted result, $a_{0, 2}$, shifted by -1.
The shifts here are circular, but it's same for linear, just more complicated in notation. The linearity (not same "linear" as in previous sentence) of subsampling is obvious and not discussed further.
Let $x$ be a length $N$ sequence, and $X[k]$ its spectrum (DFT via $e^{-j2\pi k n/N}$, no $1/N$). Let $x_{s, d}$ be the subsampling of $x$ shifted by $s$, $x_s[n] = x[n - s]$, denoted by $\texttt{sub}\{x_s, d\}[n]$.
We know, that $X_s[k] = e^{-j2\pi s k/N}X[k]$, i.e. $X_s[k] = X_0[k] e^{-j2\pi s k/N}$. Then,
$$ \begin{aligned} X_{s, d}[k] = \texttt{DFT}\{\texttt{sub}\{x_s\}\}[k] &= \frac{1}{d}\sum_{m=0}^{d-1} \exp{\left(-j2\pi \frac{s}{N} \left(k + m\frac{N}{d}\right) \right)} X_0\left[k + m\frac{N}{d}\right] \\ &= e^{-j2\pi s (k/N)} \frac{1}{d} \sum_{m=0}^{d - 1} e^{-m\cdot j2\pi(s/d)} X_0\left[k + m\frac{N}{d}\right] \\ &= e^{-j2\pi s (k/N)} \frac{1}{d} \sum_{m=0}^{d - 1} X_0\left[k + m\frac{N}{d}\right]\ \ \ (s/d \in \mathbb{Z}) \\ &= e^{-j2\pi s (k/N)} X_{0, d}[k] \end{aligned} $$
Thus, $\boxed{X_{s,d}[k] = e^{-j2\pi s k/N}X_{0,d}[k]}$. Note, this "$k$" is different from original - it spans $0$ to $N/d - 1$. Let $N_d = N/d$ - to determine how much we shift in subsampled space, rewrite as $e^{-j2\pi (s)k/N} = e^{-j2\pi (s/d) k/N_d}$. So, the subsampled sequence is shifted by $s/d$ - which is same as shifting the original sequence by $s$ (see "Basic answer"). This proves $s/d \in \mathbb{Z}$ is LTI for all $x$.
Note, this reasoning can't be recycled for convolution! Convolution and subsampling don't commute if there's aliasing, and there's also a constant scaling difference - handled next.
Let $x, h$ be length $N$ sequences, and $y[n] = (x \circledast h)[n]$ their circular convolution. Let also the $s$ in $Y$ refer only to shift in $x$ - so, $Y_s[k] = X_s[k]H[k]$ (but it also works for $XH_s$). All else is same as before.
First and foremost, we acknowledge
$$ \texttt{sub}\{x \circledast h, d\} \neq \texttt{sub}\{x, d\} \circledast \texttt{sub}\{h, d\} $$
This can be seen from the relation used in the previous section: if they were equal, then the following (their DFTs) would be equal
$$ \begin{aligned} & \left(\frac{1}{d}\sum_{m=0}^{d-1} X\left[k + m\frac{N}{d}\right]\right) \left(\frac{1}{d}\sum_{m=0}^{d-1} H\left[k + m\frac{N}{d}\right]\right) \stackrel{?}{=} \\ & \frac{1}{d}\sum_{m=0}^{d-1} X\left[k + m\frac{N}{d}\right]H\left[k + m\frac{N}{d}\right] \end{aligned} $$
yet we immediately see, one has $1/d^2$ and the other $1/d$. This aside, the proof that this isn't an equality, is that product of sums $\neq$ sum of products, and both $X$ and $H$ are fully general. So, we can't directly plug $Y=XH$ into the previous proof (i.e. $X_0\rightarrow X_0H$) and reason that a convolution is equivalently a general signal and hence the conclusion carries.
The trick is, we realize, that the "outputs" here for which LTI-ness must hold, are always subsampled - we're comparing unshifted subsampling with shifted subsampling. Let
$$ X_{0, d}[k] = \left(\frac{1}{d}\sum_{m=0}^{d-1} X\left[k + m\frac{N}{d}\right]\right) \\ H_{0, d}[k] = \left(\frac{1}{d}\sum_{m=0}^{d-1} H\left[k + m\frac{N}{d}\right]\right) $$
What we seek to show is,
$$ Y_{s,d}[k] = e^{-j2\pi s k/N}Y_{0,d}[k] $$
meaning
$$ X_{s, d}[k] H_{0, d}[k] = e^{-j2\pi s k/N} X_{0, d}[k] H_{0, d}[k] $$
yet, this is exactly same as the result of the first proof: the $H_{0, d}[k]$ are on both sides, and cancel. So, with corrected reasoning, the signal proof is recycled for convolution.
The $1/d^2$ vs $1/d$ remark is non-applicable here (terms drop, including all $1/d$s) - but it is applicable for roLTI, discussed next.
This section isn't necessary - it suffices to follow the invertibility arguments. But the section also proves another, related case, which cannot be (at least easily) deduced from same arguments. By proving for convolution, of course the case of a single signal is also proven.
Instead of proving with upsampling directly, we prove that the subsampling-convolution inequality in the previous section is equality (within $d$) only if both inputs are bandlimited, and then infer equivalence from the definition of "roLTI". To do this, we first prove that it's not an equality if only one input is bandlimited.
The identity is
$$ \sum_{k \in K} a_k b_k = \left(\sum_{k \in K} a_k\right) \left(\sum_{k' \in K} b_{k'}\right) - \sum_{k, k' \in K; k \neq k'} a_k b_{k'}, $$
and we seek to show that the corrective term is non-zero. The corrective term is
$$ \begin{aligned} & \frac{1}{d} \sum_{m,m'\in [0:d-1];\\m\neq m'} X\left[k + m\frac{N}{d}\right] H\left[k + m'\frac{N}{d}\right] \\ \end{aligned} $$
First, acknowledge again that this "$k$" differs from original, and sweeps $0$ to $N/d - 1$. For $d=2$:
So, only outer of $H$ and inner of $X$ cross (or in signal terms, low frequencies in $H$ and high in $X$), and the two $m, m'$ don't contribute to each other (to any same $k$). Since there are non-zero terms, the sum is non-zero. If, however, $X$ is also bandlimited, then the "inner of $X$" is zero, hence there is no overlap and the corrective sum drops, yielding
$$ \boxed{ \sum_{m=0}^{d-1}X\left[k + m\frac{N}{d}\right]H\left[k + m\frac{N}{d}\right] = \left(\sum_{m=0}^{d-1} X\left[k + m\frac{N}{d}\right]\right) \left(\sum_{m=0}^{d-1} H\left[k + m\frac{N}{d}\right]\right) \\ \Leftrightarrow \\ \texttt{sub}\{x \circledast h, d\} = \texttt{sub}\{x, d\} \circledast \texttt{sub}\{h, d\} \cdot d,\\ H[k']=0, X[k']=0\ |\ N/4\leq k' < N - N/4 } $$
The argument generalizes for any other $d$. This also explains my earlier "within constant scaling". Note the $\cdot d$ isn't needed for the signal-only (no convolution) case.
Refer to "LTI vs other LTI"; we've shown that $s / d \in \mathbb{Z}$ is LTI, i.e. both signal and convolution subsampling is fLTI - now we seek to show that this also guarantees they're rLTI under the same (integer ratio) constraint, i.e. $\text{fLTI} \Rightarrow \text{frLTI}$. So, we seek to show
$$ \texttt{sub}\{\texttt{shift}\{x, s\}, d\} = \texttt{shift}\{\texttt{sub}\{x, d\}, s\},\ s / d \in \mathbb{Z} \\ \Rightarrow \\ \texttt{up}\{\texttt{sub}\{\texttt{shift}\{x, s\}, d\}, d\} = \texttt{shift}\{\texttt{up}\{\texttt{sub}\{x, d\}, d\}, s\},\ s / d \in \mathbb{Z} $$
It suffices to show, that the second line is same as
$$ \texttt{up}\{\texttt{sub}\{\texttt{shift}\{x, s\}, d\}, d\} = \texttt{up}\{\texttt{shift}\{\texttt{sub}\{x, d\}, s\}, d\},\ s / d \in \mathbb{Z} $$
which is the first line wrapped with $\texttt{up}\{\cdot\}$, and if the inputs to $\texttt{up}$ are equal, then so are its outputs (one-to-one). Hence, we must show
$$ \texttt{shift}\{\texttt{up}\{\texttt{sub}\{x, d\}, d\}, s\} = \texttt{up}\{\texttt{shift}\{\texttt{sub}\{x, d\}, s\}, d\} $$
Immediately, we realize that this must be true, since both $\texttt{up}$ and $\texttt{shift}$ utilize the same information in either case - one upsamples what's already shifted, another shifts what's upsampled without the upsampler injecting information. To formalize, we consider the differences in operation:
Hence the operations are entirely commutative, the only difference is in whether multiplying or zero-stuffing happens first or second.
If subsampling is the system, then shifting its input by $s$ shifts its output by $s/d$, which in general is fractional and hence not LTI, but when it's integer, it maps directly to subsampling of a shift.
The part after "when it's integer" is true, but "shifts its output by $s/d$" needs clarification; note the proof above requires $s/d \in \mathbb{Z}$. If output was indeed straightforwardly shifted by a fraction, then fractionally shifting the output of an unshifted input would yield the shifted output, yet the proof rules this out.
Let $z[n]$ be the subsampling of $x[n - s]$. What's true is, $z[0]=x[n - s]$ and $z[1] = x[n - (s + d)]$, and $z[s_0] = x[n - (s + s_0\cdot d)]$ - or, shift in subsampled space by $s_0$ is shift in original space by $s_0\cdot d$. So, the converse is true - but the original is trickier.
Intuitively, that's because in general (non-bandlimited), subsampling of $x[n]$ and of $x[n - 1]$ yield completely different sequences in information content. For $x[n - 1]$ and even $d$, the two sequences share no points in common - but this isn't necessary. To "fractionally shift" means to utilize knowledge of points in-between - i.e., points in original space - yet the original space is lost.
It directly follows from previous paragraph, that output is shifted by $s/d$ if input is bandlimited.
Now, I also didn't say it's "false". That's because it depends on what we mean by "fractional shift" - if it means applying the Fourier shift property, then it's false. But just because we can't shift fractionally, doesn't mean it can't be done, or that in-between samples in the original sequence don't exist.
It's the "Falsy proof" but with second clause omitted ("which in general is fractional and hence not LTI"). Consider an integer $s/d$ for $a$ under "Basic answer":
$$ \begin{aligned} a_{-2} &= [2, 3, 4, 5, 6, 7, 0, 1] \\ \Rightarrow a_{-2, 2} &= [2, 4, 6, 0] \end{aligned} $$
Then, $\texttt{op}\{\texttt{shift}\{a, -2\ \text{sec}\}, 2\} = \texttt{shift}\{\texttt{op}\{a, 2\}, -2\ \text{sec}\}$, which is LTI - where the operation $\texttt{op}$ is $\texttt{sub}$, and I used "sec" for "seconds" to refer to the underlying continuous time (and assumed 1 sample = 1 sec). That's because "TI" is defined in terms of shifts with respect to the same domain, and subsampling is equivalently domain dilation.
I didn't present just this proof instead of the other, more complicated one, because frankly I didn't realize this one at first - but the other one has the benefit of proving the convolution and bandlimited cases, and revealing non-interchangeability of easy-to-confuse concepts. Otherwise, depending on what we mean by "fractional shift", since lossy-aliasing != lossy, it's possible to make it so that $\texttt{sub}\{\texttt{shift}\{a, s\ \text{sec}\}, d\} = \texttt{shift}\{\texttt{sub}\{a, d\}, s\ \text{sec}\}$ is true for all $s, d$ - that is, replace "bandlimited" with "recoverable".
Let $\texttt{up}\{(\cdot), d\}$ be DFT-upsampling of $(\cdot)$ by $d$, and $s$ have a dimensionality (see $\text{sec}$ remark under "Simple proof") but be integer in original sequence.
Now for distinguishing examples:
Code unambiguously passes for certain combinations of $s$, $d$, and bandlimit status, that this post doesn't cover, nor that I can explain, but they're pretty random-looking and probably not useful anyway. My guess is, some cisoid magic happens in the "subsampling <=> folding" relation. These cases are printed in code at the end of testing - currently:
Number of optional (unexpected) successes:
Case 0: 0 ( 0%)
Case 1: 8532 (73.8%)
Case 2: 8532 (73.8%)
Case 3: 4800 (41.5%)
Case 4: 7704 (66.7%)
Example optional successes for case 2:
real, x_bandlimited, h_bandlimited, x_only, N, d, s
(True, True, True, True, 16, 2, 1)
(True, True, True, True, 16, 2, 3)
(True, True, True, True, 16, 2, 5)
(True, True, True, True, 16, 2, 7)
(True, True, True, True, 16, 2, 9)
It's stuff I discovered while working through unnecessarily complicated math, that I found interesting, so I leave it here instead of nowhere - and they aren't totally unrelated. In summary, for integer-frequency cisoids:
Explanations follow.
(1): Subsampling an integer frequency cisoid by a factor of $d$ is same as trimming a cisoid of twice the frequency by a factor of $d$, for even-length cisoids. Conversely, trimming a cisoid by a factor of $d$ is same as subsampling a cisoid of half the frequency by a factor of $d$, for even-length cisoids, but the subsampling must be of an integer frequency, meaning the trimming be of an even frequency. I've only verified this for even lengths and even subsampling factors - unsure if it holds otherwise.
(2):
$$ \sum_{m,m'\in [0:d-1];\\m\neq m'} e^{m\cdot j2\pi f_0} = -\sum_{m=0}^{d - 1} e^{m\cdot j2\pi f_0},\quad f_0 \in \mathbb{Z} $$
It says, the sum with two non-overlapping indices equals the sum with one index. It exploits the unit circle, or zero-sum of an integer-frequency cisoid.
It arose as the corrective term in the identity for "sum of products <=> product of sums", when I tried to solve the original problem by considering the shift term alone. The non-corrective term drops (below omits one of the sums since the other is zero):
$$ \begin{aligned} \frac{1}{d}\sum_{m=0}^{d - 1} \exp{\left(-j2\pi \frac{s}{N} \left(k + m\frac{N}{d}\right) \right)} &= \frac{1}{d} e^{-j2\pi s (k/N)} \sum_{m=0}^{d-1} e^{-m\cdot j2\pi (s/d)} \\ &= \frac{1}{d} e^{-j2\pi s (k/N)} \frac{1-e^{-j2\pi s}}{1-e^{-j2\pi s/d}} \\ &= 0\quad (s\in\mathbb{Z}\Rightarrow 1 - e^{-j2\pi s} = 0) \end{aligned} $$
The corrective term is
$$ \begin{aligned} & \frac{1}{d} e^{-j2\pi s (k/N)} \sum_{m,m'\in [0:d-1];\\m\neq m'} e^{-m\cdot j2\pi (s/d)} Y_0\left[k + m'\frac{N}{d}\right] \\ &= -\frac{1}{d} e^{-j2\pi s (k/N)} \sum_{m=0}^{d - 1} e^{-m\cdot j2\pi(s/d)} Y_0\left[k + m\frac{N}{d}\right] \\ \end{aligned} $$
where $k, s, d \in \mathbb{Z}$ are of course assumed. Note that the corrective term is subtracted, so it equals the original sum, bringing us full circle.
To see the equality, first rewrite $e^{-j2\pi s (m/d)}$ - now it's clear, the iteration in $m$ is along a cisoid with integer frequency $s$. Then, for any fixed $m'$, we have
$$ Y_0\left[k + m'\frac{N}{d}\right] \sum_{m\in [0:d-1]; m \neq m'} e^{-j2\pi s (m/d)} \\ = - Y_0\left[k + m'\frac{N}{d}\right] e^{-j2\pi s (m' / d)} $$
because,
$$ \sum_{m=0}^{d - 1} e^{-j2\pi s(m/d)} = 0 $$
i.e.,
$$ e^{-j2\pi s (m' / d)} + \sum_{m\in [0:d-1]; m \neq m'} e^{-j2\pi s (m/d)} \triangleq \sum_{m=0}^{d - 1} e^{-j2\pi s(m/d)} = 0 $$
so
$$ \sum_{m\in [0:d-1]; m \neq m'} e^{-j2\pi s (m/d)} = -e^{-j2\pi s (m' / d)} $$
i.e., the sum that excludes iterating $m'$ equals the negative of iteration $m'$. This rewrites the original sum from $m'=0$ to $m' = d -1$, from there just change the variables.
It's clearer explicitly: for a fixed $m'$, we're doing some $Yx_0 + Yx_1 + ... = Y[x_0 + x_1 + ...] = Yb$, where $x_0, x_1, ...$ are values of an integer frequency complex sinusoid. We know that $a + b = 0$, where $a$ is the skipped iteration, since the full iteration sums to zero. Using $Yb = Y(a + b) - Ya$, since $a + b = 0$ --> $Yb = -Ya$.
Off-topic 1: (code may not generalize to other combos but I've done it)
import numpy as np
N, d = 14, 2
t0 = np.arange(N)/N
t1 = np.hstack([t0[:int(np.ceil(N/(2*d)))], t0[-(N//(2*d)):]])
a = np.cos(2*np.pi * fg * t0)[::d]
b = np.cos(2*np.pi * fg*d * t1)
assert np.allclose(a, b)
LTI testing: available at Github. Except, it doesn't test rLTI & roLTI, because I'm not willing to share fft_upsample at the moment. The interested reader can either find it, implement it, or "Follow" this answer until I do share it. For now, take my word for it, all assertions pass. Another option is to manually inspect the spectra using what's learned from this answer. If anyone finds an open-source fft_upsample that works in the code, feel free to comment and I'll add it in.
Is downsampling LTI for bandlimited inputs?
No.
Its linear but it's definitely does NOT meet the formal criterion for time invariance. Assuming $y[n] = T\{x[n]\}$ time invariance means
$$ T\{x[n-k]\} = y[n-k]$$
If you down-sample by $M$ and shift the input by $k$ the output shifts by $k/M$. A fractional shift is difficult to interpret in terms of a discrete system but it clearly doesn't meet the TI criterion
