Decimation filter that throws away every $N$th sample, keeping multiple samples: is it valid?

Question

While working on an existing project, I came across a piece of code that does a decimation filter in a way that I find not obvious.

In this code, every 6th sample is being thrown away. The result is being fed to FFT functions, in order to calculate magnitudes of two particular frequencies. This is done to demodulate S-FSK messages.

My question is specifically about validity of such approach when only single $N$th sample is being thrown away. Googling around yielded plenty of beginner-level articles on decimation filtering, which all suggested keeping a single $N$th sample. This makes much more sense to me, since (assuming sufficient oversampling) keeping only a single sample periodically will not change an overall form of a sine wave.

On the other hand, in the implementation that I encountered the end result is a broken wave with phase changing all the time.

Granted, we don' care about the phase in this application, only the magnitude, but still, is this valid?

If yes, what are the drawbacks or potential side effects of such approach?

Answer 1

In this code, every 6th sample is being thrown away.

Bob Adams of Analog Devices did this paper 3 decades ago that addresses exactly this situation. If you run into a pay wall, lemme know, I can get you a copy.

If you are regularly discarding every 6th sample, that means you are keeping 5 out of every 6 samples in a regular manner. Normally the sampling theorem says that you must constrain the bandwidth of the signal being sampled to less than half the sample rate:

$$ B < \frac{f_\mathrm{s}}{2} $$

We call $\frac{f_\mathrm{s}}{2}$ the "Nyquist frequency".

Adams had shown that if you limit your actual bandwidth to

$$B < \frac{N-1}{N} \times \frac{f_\mathrm{s}}{2}$$

or $\frac{N-1}{N}$Nyquist, then you can recover the original signal from the $N-1$ samples that you keep from the original $N$.

Answer 2

My question is specifically about validity of such approach

Every approach is valid if it meets the stated requirements of the application. In other words "if it works, it works" and it all depends on what "works" exactly means.

If yes, what are the drawbacks or potential side effects of such approach?

That all depends on WHY it was done this way and what the alternatives are. It's possible that someone needed a sample rate conversion by a factor of 5/6th "on a budget". This approach has the advantages that is very cheap to implement and it's causal. It has the disadvantage that it creates lots of noise and distortion.

If your application can easily tolerate that amount of noise and distortion, than going for the cheapest implementations is arguable the "best" choice in terms of development & test time and processor resources.

Answer 3

Okay, this is gonna be a long answer that will present and summarize the 1991 paper by Bob Adams of Analog Devices on this topic. I digested this a bit and will present the same concepts in the manner I would if the idea were mine.

---

So here are the prelims and definitions. Functions with brackets (like $x[n]$) are discrete time and the argument can only be an integer ($n \in \mathbb{Z}$).

Discrete-Time Fourier Transform:

$$\begin{align} \mathcal{DTFT}\Big\{x[n]\Big\} &\triangleq X(e^{j \omega}) \\ &= \sum\limits_{n=-\infty}^{\infty} x[n]\,e^{-j \omega n} \\ \end{align}$$

$X(e^{j \omega})$ is always periodic with period $2\pi$:

$$ X\big(e^{j (\omega + 2\pi)}\big) = X(e^{j \omega}) \qquad \forall \omega \in \mathbb{R} $$

and $\omega = \pm\pi$ is the Nyquist frequency expressed in terms of normalized angular frequency.

Unit discrete-time impulse function (the "Kronecker delta"):

$$ \delta[n] \triangleq \begin{cases} 1 \qquad & n = 0 \\ 0 \qquad & n \ne 0 \\ \end{cases} $$

Rectangular function:

$$ \operatorname{rect}(u) \triangleq \begin{cases} 1 \qquad & |u| < \tfrac12 \\ \tfrac12 \qquad & |u| = \tfrac12 \\ 0 \qquad & |u| > \tfrac12 \\ \end{cases} $$

Sinc function:

$$ \operatorname{sinc}(u) \triangleq \begin{cases} 1 \qquad & u = 0 \\ \frac{\sin(\pi u)}{\pi u} \qquad & u \ne 0 \\ \end{cases} $$

---

Let $0<\alpha\le 1$, then the DTFT of the sinc function (expressed as the impulse response of an ideal, brick wall, low-pass filter) is:

$$ h_\mathrm{LP}[n] = \alpha \operatorname{sinc}(\alpha n) $$

$$\begin{align} \mathcal{DTFT}\Big\{h_\mathrm{LP}[n]\Big\} &\triangleq H_\mathrm{LP}(e^{j \omega}) \\ &= \sum\limits_{m=-\infty}^{\infty} \operatorname{rect}\left(\frac{1}{\alpha} \frac{\omega - 2 \pi m}{2 \pi} \right) \\ \end{align}$$

One can see that if $\alpha = 1$ then $h_\mathrm{LP}[n] = \delta[n]$ and $H_\mathrm{LP}(e^{j \omega})=1$ for all $\omega$ and the LPF is a wire.

This ugly summation is necessary to make the DTFT periodic with period $2\pi$, but if we're willing to consider only the baseband, $|\omega|<\pi$, then we need only the $m=0$ term in the summation and this is sufficiently accurate:

$$\begin{align} \mathcal{DTFT}\Big\{h_\mathrm{LP}[n]\Big\} &\triangleq H_\mathrm{LP}(e^{j \omega}) \\ &= \operatorname{rect}\left(\tfrac{1}{\alpha} \tfrac{\omega}{2 \pi} \right) \qquad \qquad -\pi < \omega < \pi\\ \end{align}$$

So when $|\omega| < \alpha \pi$ then $H_\mathrm{LP}(e^{j \omega})=1$ and when $\alpha \pi < |\omega| < \pi$ then $H_\mathrm{LP}(e^{j \omega})=0$. A brick wall LPF with cutoff frequency at $\alpha \pi$ or $\alpha \times$Nyquist.

---

So the idea that Bob cooked up was a filter that had, for it's impulse response, coefficients that are zero for every $N$th sample. $N$ must be a finite and positive integer:

$$ N \in \mathbb{Z}, \qquad N > 0 $$

We start out with the low-pass filter above in which

$$ \alpha = \frac{N-1}{N} $$

$$ h_\mathrm{LP}[n] = \tfrac{N-1}{N} \operatorname{sinc}\big(\tfrac{N-1}{N} n \big) $$

One can see that at every multiple of $N$, except at the $0$th multiple, that the impulse response of this filter is zero.

$$\begin{align} h_\mathrm{LP}[mN] &= \tfrac{N-1}{N} \operatorname{sinc}\big(\tfrac{N-1}{N} mN \big) \qquad & m \in \mathbb{Z} \\ &= \tfrac{N-1}{N} \operatorname{sinc}\big((N-1) m \big) \\ &= \begin{cases} \tfrac{N-1}{N} \qquad & m=0 \\ 0 \qquad & m \ne 0 \\ \end{cases} \\ \end{align}$$

This has frequency response of

$$ H_\mathrm{LP}(e^{j \omega}) = \operatorname{rect}\left(\tfrac{N}{N-1} \tfrac{\omega}{2 \pi} \right) \qquad \qquad -\pi < \omega < \pi $$

So when $|\omega| < \frac{N-1}{N} \pi$ then $H_\mathrm{LP}(e^{j \omega})=1$ and when $\frac{N-1}{N} \pi < |\omega| < \pi$ then $H_\mathrm{LP}(e^{j \omega})=0$. A brick wall LPF with cutoff frequency at $\frac{N-1}{N} \pi$ or $\frac{N-1}{N} \times$Nyquist.

Now that impulse response is zero every $N$th sample except for the $0$th sample (when $n=0$). That single exception can be dealt with by subtracting a correctly scaled Kronecker impulse from it:

$$ h_\mathrm{HP}[n] = \tfrac{N-1}{N} \left( \operatorname{sinc}\big(\tfrac{N-1}{N} n \big) - \delta[n] \right) $$

This is sorta a high-pass filter and the impulse response is zero for every integer multiple of $N$.

$$ h_\mathrm{HP}[mN] = 0 \qquad \qquad \forall m \in \mathbb{Z} $$

We see that the frequency response is

$$\begin{align} H_\mathrm{HP}(e^{j \omega}) &= \operatorname{rect}\left(\tfrac{N}{N-1} \tfrac{\omega}{2 \pi} \right) - \tfrac{N-1}{N} \qquad \qquad -\pi < \omega < \pi \\ &= \begin{cases} \tfrac{1}{N} \qquad & |\omega| < \tfrac{N-1}{N}\pi \\ -\tfrac{N-1}{N} \qquad & \tfrac{N-1}{N}\pi < |\omega| < \pi \\ \end{cases} \\ \end{align}$$

So, now for the low frequencies $|\omega| < \tfrac{N-1}{N}\pi$, we want to scale it so that the gain is $1$ (like a wire), and this is our final reconstruction filter we will be using:

$$ h[n] = (N-1) \left( \operatorname{sinc}\big(\tfrac{N-1}{N} n \big) - \delta[n] \right) $$

and

$$\begin{align} H(e^{j \omega}) &= N \operatorname{rect}\left(\tfrac{N}{N-1} \tfrac{\omega}{2 \pi} \right) - (N-1) \qquad \qquad -\pi < \omega < \pi \\ &= \begin{cases} 1 \qquad & |\omega| < \tfrac{N-1}{N}\pi \\ -(N-1) \qquad & \tfrac{N-1}{N}\pi < |\omega| < \pi \\ \end{cases} \\ \end{align}$$

It is still the case that every $N$th sample of the impulse response $h[n]$ is zero.

$$ h[mN] = 0 \qquad \qquad \forall m \in \mathbb{Z} $$

So, when the missing samples line up with these zero coefficients, they don't matter.

---

So above we defined a discrete-time filter with impulse response $h[n]$ and frequency response $H(e^{j\omega})$ that has the following properties:

1. The impulse response is zero at every $N$th sample: $$h[mN]=0 \qquad \forall m \in \mathbb{Z}$$
2. The frequency response is flat and equal to $1$ for all frequencies below $\frac{N-1}{N}\times$Nyquist. $$ H(e^{j\omega}) = 1 \qquad |\omega|<\tfrac{N-1}{N}\pi $$ That means all frequency components of the input are unaffected below that threshold frequency.

Let the input to the filter be $x[n]$ having DTFT of $X(e^{j\omega})$ and the output of the filter be $y[n]$ having DTFT of $Y(e^{j\omega})$ which is:

$$\begin{align} y[n] &= h[n]*x[n] \\ & = \sum\limits_{i=-\infty}^{\infty} x[i] h[n-i] \\ & = \sum\limits_{i=-\infty}^{\infty} h[i] x[n-i] \\ \end{align}$$

$$ Y(e^{j\omega}) = H(e^{j\omega}) X(e^{j\omega}) $$

Now, here is the biggie condition: Suppose that $x[n]$ is bandlimited to $\tfrac{N-1}{N}\times$Nyquist. That is:

$$ X(e^{j\omega}) = 0 \qquad \text{for } \ \tfrac{N-1}{N}\pi \le |\omega| \le \pi $$

Then, no matter what finite value that the frequency response $H(e^{j\omega})$ takes on for those frequencies $\tfrac{N-1}{N}\pi \le |\omega| \le \pi$, the output $Y(e^{j\omega})$ will also be zero at those frequencies.

$$\begin{align} Y(e^{j\omega}) &= H(e^{j\omega}) \cdot X(e^{j\omega}) \\ &= H(e^{j\omega}) \cdot 0 \qquad \qquad \tfrac{N-1}{N}\pi \le |\omega| \le \pi \\ &= 0 \\ &= X(e^{j\omega}) \\ \end{align}$$

But below that threshold $|\omega| < \tfrac{N-1}{N}\pi$, the frequency response is $1$ so then all of the input gets passed to the output unchanged.

$$\begin{align} Y(e^{j\omega}) &= H(e^{j\omega}) \cdot X(e^{j\omega}) \\ &= 1 \cdot X(e^{j\omega}) \qquad \qquad |\omega| < \tfrac{N-1}{N}\pi \\ &= X(e^{j\omega}) \\ \end{align}$$

So it doesn't matter whether the frequency is above or below the bandlimit threshold, $\tfrac{N-1}{N}\pi$, as long as the input $x[n]$ is bandlimited to $\tfrac{N-1}{N}\times$Nyquist, the output of this filter is the same as the input:

$$ Y(e^{j\omega}) = X(e^{j\omega}) $$

or

$$ y[n] = x[n] $$

---

Now consider that we contaminate the input $x[n]$ with some crap (we shall call "$\epsilon[\cdot]$"), every $N$th sample. This contaminated signal is

$$ \hat{x}[n] \triangleq \begin{cases} x[n] \qquad \qquad & n \ne mN \quad m\in\mathbb{Z} \\ x[n] + \epsilon[m] \qquad \qquad & n = mN \\ \end{cases} $$

If $0<|\epsilon[m]|<\infty$ and $\epsilon[m]$ is unknown, this error of unknown value added to $x[mN]$ essentially makes $x[mN]$ a "missing sample".

Now consider the discrete convolution for the output $y[n]$:

$$\begin{align} y[n] &= h[n]*x[n] \\ &= \sum\limits_{i=-\infty}^{\infty} x[i] h[n-i] \\ &= \sum\limits_{k=-\infty}^{\infty} \sum\limits_{i=0}^{N-1} x[kN+i] h[n-(kN+i)] \\ &= \sum\limits_{k=-\infty}^{\infty} \left( x[kN] h[n-kN] + \sum\limits_{i=1}^{N-1} x[kN+i] h[n-(kN+i)] \right) \\ \end{align}$$

Now we have already established that when the uncontaminated (but bandlimited to $\frac{N-1}{N}\pi$) $x[n]$ goes in, exactly the same $x[n]$ comes out. What comes out when the contaminated $\hat{x}[n]$ goes into this filter?

$$\begin{align} \hat{y}[n] &= h[n]*\hat{x}[n] \\ &= \sum\limits_{i=-\infty}^{\infty} \hat{x}[i] h[n-i] \\ &= \sum\limits_{k=-\infty}^{\infty} \sum\limits_{i=0}^{N-1} \hat{x}[kN+i] h[n-(kN+i)] \\ &= \sum\limits_{k=-\infty}^{\infty} \left( \hat{x}[kN] h[n-kN] + \sum\limits_{i=1}^{N-1} \hat{x}[kN+i] h[n-(kN+i)] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( (x[kN]+\epsilon[k]) h[n-kN] + \sum\limits_{i=1}^{N-1} x[kN+i] h[n-(kN+i)] \right) \\ \end{align}$$

Note that we got rid of the hat for $x[n]$ in the last line, but now the error $\epsilon[\cdot]$ is in there. We cannot expect that $\hat{y}[n]=y[n]$ (which is $x[n]$) for every $n$, but consider the samples that are at indices that are multiples of $N$, that is $n=mN$.

$$\begin{align} y[mN] &= \sum\limits_{k=-\infty}^{\infty} \left( x[kN] h[mN-kN] + \sum\limits_{i=1}^{N-1} x[kN+i] h[mN-(kN+i)] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( x[kN] h[(m-k)N] + \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( x[kN] \underbrace{h[(m-k)N]}_{=0} + \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \right) \\ \\ \\ \hat{y}[mN] &= \sum\limits_{k=-\infty}^{\infty} \left( (x[kN]+\epsilon[k]) h[mN-kN] + \sum\limits_{i=1}^{N-1} x[kN+i] h[mN-(kN+i)] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( (x[kN]+\epsilon[k]) h[(m-k)N] + \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( (x[kN]+\epsilon[k]) \underbrace{h[(m-k)N]}_{=0} + \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \left( \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \right) \\ \\ &= y[mN] \\ &= x[mN] \\ \end{align}$$

So we have recovered our missing sample $x[mN]$ with this explicit formula:

$$\begin{align} x[mN] &= \sum\limits_{k=-\infty}^{\infty} \sum\limits_{i=1}^{N-1} x[kN+i] h[(m-k)N-i] \\ &= \sum\limits_{k=-\infty}^{\infty} \sum\limits_{i=1}^{N-1} x[kN+i] (N-1) \left( \operatorname{sinc}\big(\tfrac{N-1}{N} ((m-k)N-i) \big) - \delta[(m-k)N-i] \right) \\ &= \sum\limits_{k=-\infty}^{\infty} \sum\limits_{i=1}^{N-1} x[kN+i] (N-1) \Big( \operatorname{sinc}\big(\tfrac{N-1}{N} ((m-k)N-i) \big) - \underbrace{\delta[(m-k)N-i]}_{=0} \Big) \\ &= (N-1) \sum\limits_{k=-\infty}^{\infty} \sum\limits_{i=1}^{N-1} x[kN+i] \operatorname{sinc}\big( (N-1)(m-k-\tfrac{i}{N}) \big) \\ \end{align}$$

Answer 4

If the "S" in "S-FSK" means "spread", as in it's a very wide frequency difference compared to the data rate, and if the signal to noise ratio is high in the received signal, then just swallowing every 6$^{th}$ bit could make sense. (So, basically, what Hilmar said).

Probably what happened is that someone, at some point in the process, was in a bind for processing power and realized that taking five samples out of six would be good enough, and that's what got realized in the final software.

As a necessary digression: while the usual statement of the Nyquist theorem assumes uniform sampling, that's not necessary.

The usual statement of the Nyquist theorem is that you're taking uniform samples of the signal value at a rate of $\frac{1}{T_s}$, and that for exact reconstruction the bandwidth must be $B < \frac{1}{2 T_s}$. That's correct (and fraught with possibilities for error, because as you approach $T_s \to 2B$ your system design gets exponentially more complicated and expensive -- that seems to get glossed over in introductory DSP courses).

That's not the full Nyquist theorem. To satisfy the full Nyquist criterion you just need to take sufficiently independent samples sufficiently fast. Things that can work are sample trains of the signal and it's derivative (this is more theory than practice), or samples that come in groups of the signal filtered in multiple different ways (inphase and quadrature, for instance, or OFDM), or sample trains with hiccups (like your 5-of-6 sampling).

So, if you need to sample a 60Hz signal, then in theory you could sample the first 120 derivatives of that signal once every second, and be just fine. In practice, this would be utter insanity -- if nothing else, getting those 120 derivatives with physical hardware would be tremendously sensitive to noise, and getting them all gain matched would be impossible in practical terms.

On the other hand, quadrature sampling and OFDM are cases where you can get vectors of m independent samples at a rate of $\frac{2}{m\ B}$ and have things work just fine.

In your case you're taking 5 independent samples in every interval $6 T_s$. That means that if the original signal, before sampling, were bandlimited to $\frac{5}{6} \frac{1}{2T_s}$, you could reconstruct it. But you would do it with a polyphase filter*, not by just concatenating groups of five samples into a continuous vector without addressing the gaps.

The first enormous drawback of this approach (decimating -- or heximating -- without accounting for the gap) is that it's going to do weird things to the signal. Signals at some frequencies are just going to get translated in frequency by 6/5, some are going to get their amplitude changed, and I am pretty sure that some would get nulled out. As Hilmar pointed out, it works, so maybe it's OK.

Another, hopefully less serious, drawback of this approach is that if you were sampling the original signal then heximating**, then doing the absolute best reconstruction of it, then the result would be more noise sensitive than doing a proper reconstruction using polyphase filtering. In the case of heximation the difference won't be that much (the noise floor will raise by a factor of 6/5 in power, or about 0.8dB), but as a general technique

* You can look that up -- just search on "polyphase filter" and "sample rate conversion".

** Sorry, can't resist.