How is the DFT involved in downsampling?

Question

I recently had my first proper class of DSP and upon completing the review of Signals and Systems my professor started working on the concepts of upsampling and downsampling. Of course as is their style, they first introduce the topic and explain it in its entirety and then tell us what the topic is.

Now from what I understood during class in order to get $ y[n] = x[nM]$ (Downsampling) the intermediary step is to get:

$$y_1[n] = \frac 1M \sum_{k=0}^{M-1} x[n] e^{j\frac{2\pi}{M}kn}$$

which is basically convolving the signal in time-domain with an impulse train that has impulses at $Mn$ only. This way we can only take the $Mth$ samples. Then from here once we have zeroes between our every $M$ samples we can easily get rid of them $\left(Y(e^{j\omega}) = Y_1(e^{j\frac{\omega}{M}})\right)$ and get our desired downsampled signal. However, the expression for $y[n]$ looks an awful lot like the expression for the Discrete Fourier Transform (which I have not studied yet). So I guess my question is that what is the purpose of the DFT and is this the only way to downsample?

Answer 1

(My original post included an incorrect introduction stating that the OP's formula was similar to a moving average which caused some understandable confusion. I have since removed this while retaining the original answer to the OP's question regarding the DFT and its use in downsampling)

The OP's expression will result in $y_1[n]=0$ for all values of $n$ that are not an integer multiple of $M$, and $y_1[n]=x[n]$ otherwise, resulting in zeros between every $M$ samples. This is not the DFT and just simply down-sampling by selecting every $Mth$ samples (and needlessly zeroing the other samples that are discarded regardless).

The details below do not imply any other result for $y_1[n]$ other than what I already stated above, it's trivial to see that the samples are zero for all $n$ other than $y[n]= x[Mn]$. The OP stated this as well so I don't think that was even a question. The OP specifically asked about the similarity and difference of the formula given to the DFT, as well as how the DFT can be used for decimation applications. This is what I go on to explain below.

That said, I show the similarity and difference to the actual DFT (or here the inverse DFT), and then I will show interestingly how the DFT could actually be used in a decimation application (which combines necessary filtering together with such downsampling) as a channelizer, and how performance can be improved by combining the DFT with polyphase filtering as a high performance decimating channelizer.

The expression is almost that for the inverse DFT, but $x[n]$ in that case would refer to a frequency variable typically given as $X[k]$ where $k$ is the frequency index:

Formula as given in OP:

$$y_1[n] = \frac 1M \sum_{k=0}^{M-1} x[n] e^{j\frac{2\pi}{M}kn}$$

Formula for inverse DFT as typically written:

$$x[n] = \frac 1M \sum_{k=0}^{M-1} X[k] e^{j\frac{2\pi}{M}kn }$$

$k$ is typically used to represent the frequency index, and $n$ as the time index. Lower case $x[n]$ represents a function in the time domain, and upper case of the same letter $X[k]$ represents the same function in the frequency domain. An important difference here to note that in the actual inverse DFT (second form above), we would index through each different value in $X[k]$ as we do the summation, and rotating by the phasor $e^{j2\pi k n/M}$ which is also rotating by a fixed rotation on each increment given by whatever $n$ is for that output value: when $n=0$ for the very first sample in time, there is no rotation and we simply sum all the values in $X[k]$ and divide by $M$- an average of the frequency domain values. When $n=1$ we rotate in angular steps of $2\pi k/M$ as $k$ increments from $0$ to $M-1$, so once around the unit circle. A rotation in the frequency domain of one cycle is a shift in the time domain of one sample, so we are taking the average of the frequency domain values after the time domain waveform has been effectively shifted one sample: this is perfect for computing what $x[1]$ would be! And as we increment higher for $k=2$ we rotate twice (shift two samples in time), etc. and the result is the time domain reconstruction. That is the inverse DFT and not the formula provided by the OP:

So knowing this we can see that the equation given is NOT an inverse DFT for the key reason that for any given $n$, the value in the summation $x[n]$ is not changing. The result is simply zeroing out all but every M samples, the non-zero samples which are selected in down-sampling (select every Mth sample). When we combine the use of an anti-alias filter prior to down-sampling, the result is called "decimation".

With that, to the OP's question, here are some additional important related points for down-sampling (and decimation) and how we can actually use the DFT for a channelizer (simultaneously decimate multiple separate channels), as well as including the important "anti-alias" filter as a critical part of decimation:

To lower the sampling rate from a higher rate to an integer sub-multiple $M$, we can do this by selecting every $M$th sample, and discard the rest- this operation is called down-sampling. However, it is very important to low pass filter the higher sampled signal first, as the down-sampling operation results in aliasing of any signal or noise that is in the higher frequency regions. This is depicted by the graphic below showing the portions of spectrum that would alias into our final output if we were to (for example) down-sample with $D=9$. Here $9f_{s2}$ represents the higher sampling rate of the input, and $f_{s2}$ represents the final sampling rate of the output.

To further understand how this aliasing occurs, study the aliasing that happens in A/D Conversion, as it is the same process (resampling in A/D conversion is resampling from an infinite sampling rate basically). This post also explains the aliasing effects further.

A very simple anti-alias filter is a moving average filter, since it can be done with an FIR filter with unity gain coefficients (no multipliers) and ultimately for a decimator a very efficient moving average and down-sampler is combined to be a "CIC Filter." At this point note how the first bin ($k=0$) of a DFT is simply the sum of $M$ samples, so identical an average of those samples scaled by $M$; here below is the formula for the DFT for $k=0$:

$$X[k=0]= \sum_{n=0}^{M-1}x[n]e^{-j2\pi no/M}= \sum_{n=0}^{M-1}x[n]$$

Thus it is very clear that if we were to take a block by block DFT for every $M$ samples, the first bin is the decimated version of the input. Similarly if we were to move sample by sample through the input data, taking an M Sample DFT after each one sample shift, the first bin is identical to a moving average of the input (scaled by $M$: thus if we select every $M$th sample the mathematical operation is identical to a CIC filter as a moving average filter followed by a downsampler. [NOTE: I never said the moving average filter was a "great" filter, just a simple one, and in many applications sufficient so used due to its simplicity. For more complete details and example of the anti-aliasing provided, please see the post linked at the bottom of this post]

What is interesting is that we can instead use a bandpass filter to select any alias instead of the primary low frequency signal, which then combines frequency translation and decimation, and with the DFT we get this as well resulting in a channelizer!

As I demonstrated, the first bin of the DFT ($k=0$) operates identical to a moving average filter, if we for every output sample, take a new $M$ point DFT as we slide through a longer stream of data one sample at a time. Similarly each of the higher bins is a bandpass filter. See this post for further details. Thus we could implement a channelizer (M outputs, one for each band) by stepping the block $M$ samples at a time and then computing a new DFT, resulting in $M$ decimated outputs. Since we are going to throw away all but every $M$th sample at the output in the down-sampling process, then we only need to take the DFT after every $M$ samples and thus can use block by block DFT processing, resulting in the functional diagram below. From this we see how the DFT can be used in downsampling:

I show specifically what I mean by aliasing suppression in the plot below, which is identical to that in a first order CIC filter. The thin blue traces similar to Sinc functions (the Dirichlet Kernel specifically which is an "aliased Sinc") are the frequency response curves for each bin, due to the DFT. We see how the response is maximized on each bin center and then provides maximum rejection at the center of all other bins. As we move off center the rejection decreases. All decimation approaches will consist of an anti-alias filter that "prevents" aliasing to the degree needed for that application. All filters have a finite rejection so there will of course be aliasing in any decimation implementation, the point is to reject it sufficiently to meet SNR and other distortion requirements. I will also note that the rejection and effects of aliasing at DC are identical to every other bin, and the contribution of each bin continues to decrease as the number of bins is increased.

When the anti-alias performance of a simple moving average is not sufficient, aliasing rejection is further improved by combining polyphase partitioned filters with an IDFT; for complete details of those and other implementations using the DFT and IDFT for downsampling, I recommend fred harris' book "Multirate Signal Processing for Communications Systems" as well as V.V.T's answer at this post as well as this related post demonstrating both approaches:

Use of DFT for Decimating Channelizers

Answer 2

The DFT's downsampling relation is described by Subsampling in time <=> Folding in Fourier. Combined with frequency domain multiplication, it enables efficient lowpass or bandpas filtering, among other things, as described in that post. That's all there is to say about the title of your post. Your equation, however, has little to do with DFT.

the intermediary step

No. This step does nothing for $y[n]$ that $y[n] = x[nM]$ doesn't.

which is basically convolving the signal in time-domain with an impulse train that has impulses at $Mn$ only

No, it's multiplying. Notice that $x[n]$ can be moved outside the summation; the result is a convenient mathematical expression for unit impulses spaced by $M$ samples, which then multiply $x[n]$. In the DFT or iDFT, the summation isn't factorable - that makes all the difference.

we can easily get rid of them $\left(Y(e^{j\omega}) = Y_1(e^{j\frac{\omega}{M}})\right)$

No, this confuses time subsampling with frequency subsampling, and the expression isn't reducible like this in frequency.

Re: Dan's answer

The resulting expression when only the decimated output values are considered (n=0,M,2M…) is identical to a block by block average over M samples as a moving average filter

No it's not. Any part of the answer that explains $y[n]$ as such is hence invalid.

channelizer

This is a fancy way of saying STFT with maximum hop size and rectangular windowing, and it has nothing to do with downsampling or preventing aliasing. Even if Dan's description was amended to reduce aliasing, it'd still have nothing to do with 1D-to-1D downsampling of the original signal.

Dan's proposal does opposite of what's advertised: every output channel suffers extreme aliasing. The channel at DC suffers maybe acceptable aliasing. If the intended "relevance" is that we can get $M$ separate "views" into original, said views are completely useless for most purposes.

Dan's updates after my post make his answer dangerously misleading. Unless one already knows what Dan is referring to, personally I recommend not reading it at all.

Proving aliasing in "channelizer"

Here's the frequency response only for non-negative channels:

This is with $N=256$ and $M=8$. Here's just one channel, without subsampling:

Now imagine taking the right half and adding it to left half. That's the effect of subsampling by 2, as shown in the top linked post. Now do it again, and again, to obtain subsampling by 8. This is introducing more aliasing than the original signal could ever dream of.

There's also phase, and Nyquist-related aliasing, that this doesn't cover. For a fuller analysis, refer to this post.

Plots are reproduced with this code, that's a modification of my code to this article (which also explains this subsection).