Introduction to Systems

Discrete-Time Integrators

A discrete-time integrator implemented through a forward difference and a backward difference technique

An integrator is a very important filter that proves useful in implementation of many blocks of a communication receiver. In continuous-time case, an integrator finds the area under the curve of a signal amplitude. A discrete-time system deals with just the signal samples and hence a discrete-time integrator serves the purpose of collecting a running sum of past samples for an input signal. Looking at an infinitesimally small scale, this is the same as computing the area under the curve of a signal sampled at an extremely high rate. For the following discussion, we assume a normalized sample time T_S =1.

For an input s[n] and output r[n], there are three main methods to implement a discrete-time integrator.

[Forward difference] The forward difference integrator is realized through the following equation.

(1)   \begin{align*}                 r[n] &= \sum _{i=-\infty}^{n} s[i] \nonumber \\                      &= \sum _{i=-\infty}^{n-1} s[i] + s[n] \nonumber \\                      &= r[n-1] + s[n]              \end{align*}

For obvious reasons, the running sum r[n-1] is a common component in all types of integrators; the differentiating factor is the term added to r[n-1] as a replacement for area under the curve. For a forward difference integrator, this additional term is s[n] — the current input. This is drawn in Figure below. Notice that the forward difference integrator computes its output after the arrival of the current sample, i.e., at time n.

A discrete-time integrator implemented through a forward difference and a backward difference technique

The block diagram shown in the figure is used to implement a forward difference integrator in loop filter portion of a Phase Locked Loop (PLL) that is employed in carrier phase and symbol timing synchronization systems of a digital communication system.

[Backward difference] The backward difference integrator is realized through the following equation.

(2)   \begin{equation*}                 r[n] = r[n-1] + s[n-1]             \end{equation*}

Here, the term added to r[n-1] is the previous input s[n-1]. This is also illustrated in the Figure above. In contrast to the forward difference case, the backward difference integrator can compute its output after the last sample, i.e., at time n-1. This minor dissimilarity plays a huge role in analyzing the performance of the actual discrete-time system where the integrator is employed (refer to any DSP text to understand the role of poles and zeros of a system).

The block diagram shown in the figure is used to implement a backward difference integrator in a Numerically Controlled Oscillator (NCO) of a Phase Locked Loop (PLL).

[Average of a backward and forward difference] A third integrator can be implemented as

    \begin{equation*}             r[n] = r[n-1] + \frac{1}{2}\left(s[n-1]+s[n]\right)         \end{equation*}

Why is integrator a lowpass filter?


An integrator maintains a running sum of past samples. In time domain, the summation of a large number of values tends towards a mean value. When some numbers are small, some are large and some are in between, their running sum naturally pulls large variations towards the middle, thus smoothing them out. Large variations represent high frequencies and smoothing out such variations is the function of a lowpass filter.

To verify this fact in frequency domain, first consider that when an impulse is given as input to an integrator in time domain, it again forms a running sum of the impulse, which is a unit step signal. Hence, its impulse response is a unit step signal. The unit step signal is very wide in time domain, so it must be narrow in frequency domain. Consequently, it filters out the high frequency terms and passes only the low frequency content. We can say that it acts as a lowpass filter.

Moving Average Filter

Coefficients of a moving average filter in time domain

The most commonly used filter in DSP applications is a moving average filter. In today’s world with extremely fast clock speeds of the microprocessors, it seems strange that an application would require simple operations. But that is exactly the case with most applications in embedded systems that run on limited battery power and consequently host small microcontrollers. For noise reduction, it can be implemented with a few adders and delay elements. For lowpass filtering, the excellent frequency domain response and substantial suppression of stopband sidelobes are less important than having a basic filtering functionality, which is where a moving average filter is most helpful.

As the name implies, a length-M moving average filter averages M input signal samples to generate one output sample. Mathematically, it can be written as

    \begin{align*}       r[n] &= \frac{1}{M} \sum \limits_{m=0}^{M-1} s[n-m]\\            &= \sum \limits_{m=0}^{M-1} \frac{1}{M} \cdot s[n-m]  \\            &= \sum \limits_{m=0}^{M-1} h[m] s[n-m]     \end{align*}

where as usual, s[n] is the input signal, r[n] is the system output and h[n] is the impulse response of the moving average filter given by

    \begin{equation*}         h[n] = \frac{1}{M}, \qquad n = 0,1,\cdots,M-1     \end{equation*}

The above discussion leads us to the following observations.

Impulse response From the above equation, the impulse response of a moving average filter can be seen a rectangular pulse which is drawn in Figure below. In terms of filtering, a constant amplitude implies that it treats all input samples with equal significance. This results in a smoothing operation of a noisy input [1].

Coefficients of a moving average filter in time domain

To see an example of noise reduction in time domain, consider a pulse in Figure below in which random noise is added. Observe that a longer filter performs a better smoothing action on the input signal. However, the longer length translates into wider edges during the transition.

A noisy signal and its smoothed version for two different filter lengths

Magnitude response As a consequence of its flat impulse response, the frequency response of a moving average filter is the same as the frequency response of a rectangular pulse, i.e., a sinc signal. The magnitude response of such a filter is drawn in Figure below. It is far from an ideal brickwall response of an ideal lowpass filter but the magnitude response still has a shape that passes low frequencies with some distortion and attenuates high frequencies to some extent.

Frequency response of a moving average filter

Due to this insufficient suppression of sidelobes, it is probably the worst lowpass filter one can design. On the other hand, as described earlier, its simplicity of implementation makes it useful for many applications requiring high sample rates including wireless communication systems. In such conditions, a cascade of moving average filters is implemented in series such that each subsequent filter inflicts more suppression to higher frequencies of the input signal.

Due to the inverse relationship between time and frequency domains, a longer filter in time domain produces a narrower frequency response. This can be observed in the impulse and frequency response for filter lengths M=5 and M=11.

Phase response As described in the article on FIR filters, a moving average filter generates a delay of 0.5(M-1) samples (for odd M) in the output signal as compared to an input signal. This delay manifests itself in the form of a linear phase in frequency domain.

Convolution Notice from above equations that the filter output is a convolution of the input signal with the impulse response of the moving average filter, a rectangular pulse.

References

[1] The Scientist and Engineer’s Guide to Digital Signal Processing (1st Edition), California Technical Pub., 1997.

Finite Impulse Response (FIR) Filters

Time and frequency response of a lowpass FIR filter designed with Parks-McClellan algorithm for N=33

We learned in the concept of frequency that most signals of practical interest can be considered as a sum of complex sinusoids oscillating at different frequencies. The amplitudes and phases of these sinusoids shape the frequency contents of that signal and are drawn through magnitude response and phase response, respectively. In DSP, a regular goal is to modify these frequency contents of an input signal to obtain a desired representation at the output. This operation is called filtering and it is the most fundamental function in the whole field of DSP.

From this Eq, one can observe that it is possible to design a system, or filter, such that some desired frequency contents of the signal can be passed through by specifying H(k) \approx 1 for those values of k, while other frequency contents can be suppressed by specifying H(k) \approx 0. The terms lowpass, highpass and bandpass then make intuitive sense that low, high and band refer to desired spectral contents that need to be passed while blocking the rest.

Filter Frequency Response


Figure below shows the magnitude response |H(F)| (as a function of continuous frequency) of an ideal lowpass filter. A lowpass filter passes frequencies near 0 while blocks the remaining frequencies.

An ideal lowpass filter and aliases

As explained in the discussion about sampling, in a continuous frequency world, the middle filter is all that exists. However, in a sampled world, the frequency response of the filter — just like a sampled signal — repeats at intervals of F_S (sampling frequency) and there are infinite spectral replicas on both sides shown with dashed lines. The primary zone of the sampled frequency response is from -0.5F_S to +0.5F_S shown within dotted red lines with its spectral contents in solid line. From this Eq on sampling, we know that filter response near 0 is the same as that near F_S and the filter response around +0.5F_S is the same as around -0.5F_S.

On the other hand, a highpass filter passes both positive and negative frequencies near 0.5F_S while blocking everything else. Observe the symmetry around \pm 0.5F_S in Figure below. Similarly, bandpass filters only pass a specific portion of the spectrum.

An ideal highpass filter and aliases

Filter Impulse Response


Since a filter is a system that modifies the spectral contents of a signal, it naturally has an impulse response h[n] just like a system. When this impulse response is finite, the filter is called Finite Impulse Response (FIR). Impulse response h[n] of an FIR filter, i.e., the sequence of its coefficients, can be found by taking the iDFT of the frequency response H[k].

FIR Filter Design


Understanding of FIR filter design starts from the following argument for continuous frequency domain. Since ideal filters have unity passband and zero stopband and a direct transition between these two, this finite frequency domain support produces an impulse response h[n] in time domain which is infinite and hence non-realizable, as explained in the discussion about frequency. When this h[n] is truncated on the two sides, it is equivalent to multiplication of the ideal infinite h[n] by a rectangular window w_{\text{Rect}}[n] as

    \begin{equation*}         h_{\text{Truncated}}[n] = h[n] \cdot w_{\text{Rect}}[n]     \end{equation*}

As mentioned previously, windowing is a process of pruning an infinite sequence to a finite sequence. We know that multiplication in time domain is convolution in frequency domain, and that the spectrum of a rectangular window is a sinc signal. Thus, convolution between the spectrum of the ideal filter and this sinc signal generates the actual frequency response H(F) of the filter.

    \begin{equation*}         H(F) = \text{Ideal brickwall response} * sinc(F)     \end{equation*}

This sinc signal has decaying oscillations that extend from -\infty to \infty. Convolution with these decaying oscillations is the reason why a realizable FIR filter always exhibits ripple in the passband and in the stopband as well as a transition bandwidth between the two. Naturally, this transition bandwidth is equal to the width of the sinc mainlobe, which in turn is a function of the length of the rectangular window in time. Truncation with a rectangular window causes the FIR filters to have an impulse response h[n] of finite duration.

Sidelobes appearing in the stopband of the filter, transition bandwidth and ripple in the passband characterize the performance of a filter. These parameters are usually specified in terms of decibels (or dBs). The filter magnitude in dB is given by

    \begin{equation*}         |H(F)|_{dB} = 20 \log_{10} |H(F)|     \end{equation*}

It should be remembered that the multiplicative factor 20 above is replaced by 10 for power conversions (because power is a function of voltage or current squared). The advantage of converting magnitudes to dBs for plotting purpose is that very small and very large quantities can be displayed clearly on a single plot.

A filter is designed in accordance with the parameters in Table below and shown graphically in the next Figure.

Lowpass filter specifications table

A lowpass filter marked with specs frequencies and their amplitudes

FIR filter design basically requires finding the values of filter taps (or coefficients) that translate into a desired frequency response. Many software routines are available to accomplish this task. A standard method for FIR filter design is the Parks-McClellan algorithm. The Parks-McClellan algorithm is an iterative algorithm for finding the optimal FIR filter such that the maximum error between the desired frequency response and the actual frequency response is minimized. Filters designed this way exhibit an equiripple behavior in their frequency responses and are sometimes called equiripple filters (equiripple implies equal ripple within the passband and within the stopband, but not necessarily equal in both bands).

Example


We design a lowpass FIR filter through Parks-McClellan algorithm in Matlab that meets the following specifications.

    \begin{align*}             &\text{Sample rate}, \quad F_S = 12\: \text{kHz} \\             &\text{Passband edge} \quad F_{\text{pass}} = 2\: \text{kHz} \\             &\text{Stopband edge} \quad F_{\text{stop}} = 3\: \text{kHz} \\             &\text{Passband ripple} \quad \delta_{\text{pass,dB}} = 0.1\: \text{dB} \\             &\text{Stopband ripple} \quad \delta_{\text{stop,dB}} = -60\: \text{dB}           \end{align*}

First, noting from above Figure that the maximum amplitude in the passband of frequency domain is 1+\delta_{\text{pass}}, we can convert dB units into linear terms as

    \begin{align*}       20 \log_{10} \big(1+\delta_{\text{pass}}\big) &= \delta_{\text{pass,dB}} \\       1+\delta_{\text{pass}} &= 10^{\delta_{\text{pass,dB}}/20} \\         \delta_{\text{pass}} &= 0.012     \end{align*}

For the stopband,

    \begin{align*}       20 \log_{10} \delta_{\text{stop}} &= \delta_{\text{stop,dB}} \\       \delta_{\text{stop}} &= 10^{-60/20} = 0.001     \end{align*}

A Matlab function firpm() returns the coefficients of a length N+1 linear phase (which implies real and symmetric coefficients, see the concept of phase) FIR filter. For finding an approximate N, the Matlab command firpmord() — which stands for FIR Parks-McClellan order — can be used.

    \begin{equation*}         N = \text{firpmord} ([F_{\text{pass}}\quad F_{\text{stop}}],[1 \quad 0],[\delta_{\text{pass}} \quad \delta_{\text{stop}}],F_S)     \end{equation*}

where the vector [1\quad 0] specifies the desired amplitude pointing towards passband and stopband, respectively. For our example,

    \begin{equation*}         N = \text{firpmord} ([2000\quad 3000],[1 \quad 0],[0.012 \quad 0.001],12000) =29     \end{equation*}

A word of caution: it can raise confusion in the sense that firpmord() uses the band widths for input amplitudes while firpm() uses the band edges. For example, firpmord() uses the desired amplitude vector [1 \quad 0] to denote passband and stopband. On the other hand, firpm() uses the amplitude vector [1\quad  1\quad 0 \quad 0] to denote band edges.

For a trial value of N=29,

    \begin{align*}       h &= \text{firpm}(N-1,[0 \quad F_{\text{pass}} \quad F_{\text{stop}} \quad 0.5F_S]/(0.5F_S),[1\quad 1\quad 0\quad 0]) \\         &= \text{firpm}(28,[0 \quad 2000 \quad 3000 \quad 6000]/6000,[1\quad 1\quad 0\quad 0])     \end{align*}

where the vector [1\quad 1 \quad 0 \quad 0] specifies the desired amplitude at band edges of passband and stopband, respectively. The normalization by 0.5 F_S is of no significance as it is just a requirement of this function. Furthermore, the filter is designed with N+1 coefficients by Matlab and that is why we used N-1 as its argument.

If you try this little code, you will see that the result thus obtained had stopband levels of only around -45 dB and a passband ripple of 0.05 dB. The sidelobes do not meet the 60 dB attenuation while the passband ripple exceeds the requirement of 0.1 dB. Therefore, N can be increased until 60 dB attenuation is achieved for N = 41. To avoid excessive filter length increase which increases the processor workload, different weighting factors can be employed for passband and stopband ripples as well. For example, a weighting factor of 1:10 for passband to stopband ([1 \quad 10] in vector form) yields,

    \begin{equation*}       h = \text{firpm}(28,[0 \quad 2000 \quad 3000 \quad 6000]/6000,[1\quad 1\quad 0\quad 0],[1 \quad 10])     \end{equation*}

In this case, the sidelobe level reached -55 dB, a 10 dB improvement over the same N=29 when there were no band weights. Next, the filter order is gradually increased until it meets the desired specifications at N = 33, a difference of 8 taps from N=41 and no weighting factors. Such trial and error approach is often necessary to design a filter with given specifications.

The impulse response h[n] and frequency response 20 \log_{10} H(F) of the final filter are drawn in Figure below. Note the stopband attenuation of 60 dB and passband ripple within 0.1 dB.

Time and frequency response of a lowpass FIR filter designed with Parks-McClellan algorithm for N=33

Filtering as Convolution


We also learned above that the higher the number of taps, the more closely the filter follows an expected response (e.g., lowpass, highpass, bandpass, etc.) at a cost of increased computational complexity. Therefore, the defining feature of an FIR filter is the number of coefficients that determines the response length, the number of multipliers and the delay in samples required to compute each output. An FIR filter is built of multipliers and adders. A delay element, which is just a clocked register, is used between coefficients. A single delay element is represented with a D_1 symbol. Figure below illustrates a 5-tap FIR filter.

FIR filter implementation

Taking into account the structure in Figure above, this 5-tap filter output can be mathematically written as

    \begin{equation*}         r[n] = h_0  s[n] + h_1  s[n-1] + h_2  s[n-2] + h[3]  s[n-3] + h[4]  s[n-4]     \end{equation*}

For a general FIR filter of length N,

(1)   \begin{equation*}         r[n] = \sum \limits _{m=0} ^{N-1} h[m] s[n-m]     \end{equation*}

If this equation sounds familiar to you, you guessed it right. Compare it with convolution Eq, and this is nothing but a convolution between the input signal and the filter coefficients. Just like convolution, at each clock cycle, samples from input data and filter taps are cross-multiplied and the outputs of all multipliers are added to form a single output for that clock cycle. On the next clock cycle, the input data samples are shifted right by 1 relative to filter taps, and the process continues.

Filter Delay


Due to convolution, the length of the output signal r[n] is given as in the article on convolution, which is repeated below.

    \begin{equation*}         \textmd{Length}\{r[n]\} = \textmd{Length}\{s[n]\} + \textmd{Length}\{h[n]\} - 1     \end{equation*}

which is longer than the input by \textmd{length}\{h[n]\} - 1 samples. We can say that every FIR filter causes the output signal to be delayed by a number of samples given by

(2)   \begin{equation*}         D = \frac{\textmd{Length}\{h[n]\} - 1}{2}     \end{equation*}

As a result, D samples of the output at the start — when the filter response is moving into the input sequence — can be discarded. Similarly, the last D samples — when the filter is moving out of the input sequence — can be discarded as well.

Filter Phase Response


As discussed in the concept of phase, the phase response of filters is also important because the shape of the signal in time domain is determined by phase alignment of its constituent complex sinusoids. A filter with unity magnitude across the spectrum does not affect the amplitude of those sinusoids but if its phase response is not linear, it can completely change the shape of the filtered signal due to different sinusoids adding up with different phases at the output. A linear phase on the other hand just delays the output without distorting it. In the post on convolution, we discussed that the output response of an LTI system to a complex sinusoid at its input is nothing but the same complex sinusoid at the same frequency but with different amplitude and phase. If this phase of the filter is a linear function of frequency k, then rearranging Eq on time and phase shifts,

    \begin{flalign*}         \begin{aligned}               \textmd{Phase shift} \quad - \textmd{constant} \cdot k \cdot \textmd{delay} ~& \xrightarrow{\hspace*{1cm}}~ s[n-\textmd{delay}] \quad \textmd{Time shift}         \end{aligned}     \end{flalign*}

for all complex sinusoids in its spectrum. Fortunately, maintaining a linear phase response in FIR filters is straightforward by making the coefficients symmetric, see this Figure as an example. This is one of the main reasons for widespread use of FIR filters in DSP applications.

What is a System?

A system with input and output

In wireless communications and other applications of digital signal processing, we often want to modify a generated or acquired signal. A device or algorithm that performs some prescribed operations on an input signal to generate an output signal is called a system.

In a previous article about transforming a signal, we saw how a signal can be scaled and time shifted, or added and multiplied with another signal. These are all examples of a system. Our main focus in these articles will be on a particular class of systems which are linear and time-invariant.

Additive White Gaussian Noise (AWGN)

Computing noise power within a specified bandwidth

The performance of a digital communication system is quantified by the probability of bit detection errors in the presence of thermal noise. In the context of wireless communications, the main source of thermal noise is addition of random signals arising from the vibration of atoms in the receiver electronics.

It is called additive white Gaussian noise (AWGN) due to the following reasons:

[Additive] The noise is additive, i.e., the received signal is equal to the transmitted signal plus noise. This gives the most widely used equality in communication systems.

(1)   \begin{equation*}                 r = s + w             \end{equation*}

and is shown in Figure below. Remember that the above equation is highly simplified due to neglecting every single imperfection a Tx signal encounters, except the noise itself.

A received signal as an addition of transmitted signal and noise

Moreover, this noise is statistically independent of the signal.

[White] Just like the white color which is composed of all frequencies in the visible spectrum, white noise refers to the idea that it has uniform power across the whole frequency band. As a consequence, the spectral density of white noise is ideally flat for all frequencies ranging from -\infty to +\infty.

Constant spectral density of white noise

[Gaussian] The probability distribution of the noise samples is Gaussian, i.e., in time domain, the samples can acquire both positive and negative values and in addition, the values close to zero have a higher chance of occurrence while the values far away from zero are less likely to appear. This is shown in Figure below.

Gaussian distribution of noise

As a result, the time domain average of a large number of noise samples is zero. We will use this attribute over and over again throughout the system design.


Nyquist investigated the properties of thermal noise and showed that its power spectral density is equal to k \times T, where k is a constant and T is the temperature in Kelvin. As a consequence, the noise power is directly proportional to the equivalent temperature and hence the name thermal noise. Historically, this constant value indicated in noise spectral density Figure is denoted as N_0/2 Watts/Hz.

When we view the constant spectral density as a rectangular sequence (we do not discuss random sequences here, so this discussion is just for a general understanding), its iDFT must be a unit impulse. Furthermore, in the article on correlation, we saw that the iDFT of the spectral density is the auto-correlation function of the signal. Combining these two facts, an implication of a constant spectral density is that the autocorrelation of the noise in time domain is a unit impulse, i.e., it is zero for all non-zero time shifts.

In words, each noise sample in a sequence is uncorrelated with every other noise sample in the same sequence.

In reality, the ideal flat spectrum from -\infty to +\infty is true for frequencies of interest in wireless communications (a few kHz to hundreds of GHz) but not for higher frequencies. Nevertheless, every wireless communication system involves filtering that removes most of the noise energy outside the spectral band occupied by our desired signal. Consequently after filtering, it is not possible to distinguish whether the spectrum was ideally flat or partially flat outside the band of interest. To help in mathematical analysis of the underlying waveforms resulting in closed-form expressions — a holy grail of communication theory — it can be assumed flat before filtering as shown in noise spectral density Figure.

For a discrete signal with sampling rate F_S, the sampling theorem dictates that the bandwidth is constrained by an above mentioned lowpass filter within the range \pm F_S/2 to avoid aliasing. This filter is an ideal lowpass filter with

    \begin{equation*}         H(F) =             \begin{cases}             1 &  -F_S/2 < F < +F_S/2 \\             0 &  \textmd{elsewhere}             \end{cases}     \end{equation*}

The resulting in-band power is shown in red in Figure below, while the rest is filtered out.

Computing noise power within a specified bandwidth

As with all densities, the value N_0 is the amount of noise power P_w per unit bandwidth B.

(2)   \begin{equation*}         N_0 = \frac{P_w}{B}     \end{equation*}

Plugging B=F_S/2 in the above equation, the noise power in a sampled bandlimited system is given as

(3)   \begin{equation*}          P_w = N_0\cdot \frac{F_S}{2}     \end{equation*}

Thus, the noise power is directly proportional to the system bandwidth at the sampling stage.

Fooling the Randomness


As we will see in our discussions about communication systems, a communication signal has a lot of structure in its construction. However, due to Gaussian nature of noise acquiring both positive and negative values, the result of adding a large number of noise-only samples tends to zero.

    \begin{equation*}             \sum _n w[n] \rightarrow 0 \qquad \text{for large} ~ n         \end{equation*}

Therefore, when a noisy signal is received, the Rx exploits this structure through correlating with what is expected. In the process, it estimates various unknown parameters and detects the actual message while averaging over a large number of observations. This correlation brings order out and averaging mitigates the harsh effects of noise.

An AWGN channel is the most basic model of a communication system. Even simple practical systems suffer from various kinds of imperfections in addition to AWGN. Examples of systems operating largely in AWGN conditions are space communications with highly directional antennas and some point-to-point microwave links.