Downsampling a signal with decimate

Question

I'm experimenting with decimating a signal, in this case a unit impulse.

I'm using Python, with pylab. First, I create a unit impulse, and decimate it by 5.

x = r_[zeros(0), 1, zeros(100)]
N = 2 ** 14
q = 5

y = decimate(x, q, ftype="fir")
subplot(211)
title("Original")
stem(range(len(x)), x)
subplot(212)
title("Decimated - FIR")
stem(range(len(y)), y)

figure()
subplot(211)
semilogx(log(abs(fft(x, N))))
subplot(212)
y = decimate(x, q, ftype="fir")
semilogx(log(abs(fft(y, N))))

This results with the following plots

I then add a few samples of delay before the impulse, by changing x to:

x = r_[zeros(3), 1, zeros(100)]

This results in the following plots

In the second set of plots, the resulting decimated signal is no longer a single sample, but has been distorted.

If I delay the signal with 5 - and any multiple of q - samples, I get the first set of plots again.

The source code for the decimate function is, https://github.com/scipy/scipy/blob/master/scipy/signal/signaltools.py#L1570

def decimate(x, q, n=None, ftype='iir', axis=-1):
    if not isinstance(q, int):
        raise TypeError("q must be an integer")

    if n is None:
        if ftype == 'fir':
            n = 30
        else:
            n = 8

    if ftype == 'fir':
        b = firwin(n + 1, 1. / q, window='hamming')
        a = 1.
    else:
        b, a = cheby1(n, 0.05, 0.8 / q)

    y = lfilter(b, a, x, axis=axis)

    sl = [slice(None)] * y.ndim
    sl[axis] = slice(None, None, q)
    return y[sl]

I'm using a fir low pass filter before decimating, the impulse response of the filter is

This explains why the impulse is distorted when there is a delay, the decimation is selecting parts of the impulse response, when the delay is a multiple of the decimation, it only selects the zero's of the impulse response, and one non-zero sample at the peak.

Is there a way to decimate a unit sample with an arbitrary delay, which results in a scaled unit sample output?

Answer 1

You seem to understand correctly what is going on. However, I'm not sure what you expect to get. Take your initial example. Let your input signal be $x[n]$:

$$ x[n] = \delta[n] $$

The first step in the decimation process is that the input signal is convolved with the anti-aliasing filter's impulse response $h[n]$:

$$ \begin{align} x_f[n] &= x[n] * h[n] \\ &= \delta[n] * h[n] \\ &= h[n] \end{align} $$

Next, the filtered signal is downsampled by a factor of $q$ ($5$ in your example).

$$ x_d[n] = x_f[qn] = h[qn] $$

As you noted, for FIR filters whose orders are a multiple of $q$ (actually, since the filter is linear-phase, the order need only be a multiple of $\frac{q}{2}$), the taps at time delays $qn$ are zero for all $q \ne 0$. Therefore, $x_d[n]$ is nonzero only at $n=0$, as you found.

When you apply a time delay to the input impulse $x[n]$, then the filtered output is just delayed by the same amount, since the filter is linear and time invariant (LTI):

$$ x_f[n] = h[n-D] $$

$$ x_d[n] = x_f[qn] = h[qn-D] $$

Again, as you noted, this has the effect of plucking out a different set of taps from the filter's response, such that the decimated output signal is no longer zero for all but one sample (i.e. it doesn't look like an impulse any more). This is to be expected. Why?

Remember that the discrete-time Fourier transform (DTFT) of a discrete impulse is flat (in magnitude) in the frequency domain. If $x_d[n]$ is to be equivalent to a delayed impulse, it must also have flat magnitude in the frequency domain. However, its DTFT just a scaled copy of the filter's frequency response:

$$ x_d[n] = h[qn-D] \Leftrightarrow e^{-j\omega D} H\left(\frac{\omega}{q}\right) $$

where $H(\omega)$ is the filter's frequency response. In order for the decimated output $x_d[n]$ to be equal to a delayed impulse, the filter must have a brick-wall response that is perfectly flat across its passband (up to the post-decimation Nyquist frequency) and zero everywhere else (so no aliasing leaks back in after downsampling, making the result non-flat). This isn't realizable unless you have an infinite amount of time and resources.

Since the filter is the source of the "distortion" that you don't want, you might consider trying the process again without a filter. But, consider what you would get then:

$$ x_f[n] = x[n] = \delta[n-D] $$

$$ x_d[n] = x_f[qn] = \delta[qn-D] $$

If $q$ is not a multiple of $D$, then $x_d[n] = 0\ \forall\ n$, which is probably not what you would want either.