Digital Signal Processing in One Lesson

9 min readJul 26, 2017

Some years ago, I worked at a position where we regularly had new students join our team as research assistants. I found that often the students needed a refresher on digital signal processing. Thus, I made a MATLAB script to illustrate key concepts in digital signal processing. I’ve turned the script into a publication here on medium, but I have embedded the snippets of code needed to build the figures and videos shown.

This code demonstrates the use of the Discrete Fourier Transform (DFT) via the Fast Fourier Transform (FFT). Note that the FFT only operates on signals with a total number of samples equal to an integer power of 2. If a signal with a different value of samples is input, MATLAB takes the nearest power of 2 and then interpolates the result to give the spectral signal of the input’s length.

Sampling Periodic Functions and FFTs

Here’s an example of sampling a continuous signal to produce a digital signal. Note that “continuous signals” can’t really be created in MATLAB, so we model a continuously sampled signal with a very finely sampled digital signal.

So now let’s look at the Fourier transform of this signal.

The Fourier transform of y is meaningless unless we have a way to interpret the output! We know what frequencies the k values correspond to so we can take advantage of this fact. We also know that the amplitude of signals is multiplied by N in the discrete frequency domain. Therefore, we can scale down by N.

We observe that there is no sample at exactly 1 Hz, so we get a tapering effect. The nearest value to 1 Hz has the largest value and its value isn’t quite 0.5 (the amplitude of y). So can we get closer to the spectrum we expect to see? The absolute value of the Fourier Transform of a sinusoid with frequency fc is 0.5 ×dirac delta function center on ±fc. (Therefore, the spectrum we expect to see has an impulse with height 0.5 at ± 1 Hz.) First we will address the peak at 9 Hz.

Aliasing

Why does it look like there is a signal at 9 Hz? The reason for this is aliasing, and can be predicted by using the Nyquist Rate. The Nyquist rate says that the sampling frequency must be greater than or equal to the bandwidth of the signal to avoid aliasing. i.e. Fs ≥ 2×B. The bandwidth is essentially the highest frequency component of a signal (for us 1 Hz). Below a demonstration of aliasing is given. p = “prime”.

In the first figure we can the plot that we expect. In the second plot, we show that a 9 Hz signal with a sampling frequency of 10 samp/sec looks just like a 1 Hz signal with this sampling frequency! Moreoever a 10 Hz signal looks like a DC signal! Thus the word aliasing: when not sampled at a high eough sampling frequency (undersampling) these signals disguise (alias) themselves as other signals. e.g. the 10 Hz signals looks DC. So firstly, is there a way to predict what a signal will alias itself as given a sampling frequency? Yes, and it is given by the following equation.

|fₐ| = |f — mFs|, m is in the set of all integers. fₐ is the frequency that the aliased signal will look like. f is the frequency of the real signal Fs is the sampling frequency. Using this equation can be a little tricky. You want to find the first value of m that will make f ₐ ≤ Fs/2. E.g. fₛ = 10, f = 9. The first value that makes fₐ ≤ Fs/2 = 5 Hz is m = 1. Therefore, |fₐ| = |f — mFs| = |9–1×10| = |-1| = 1. Notice that this is exactly what we observed! The 9 Hz signals appears to be a 1 Hz signal at a sampling frequency of 10 samp/sec.

This phenomenon is exactly what the Fourier transform is showing when it displays a peak a 1 Hz and 9 Hz. It is saying that there could be a signal at 1 Hz or 9 Hz. Therefore, it must be known apriori what the bandwidth of the signal you are concerned with is. One may now observe that the FFT will be periodic with a period of Fs. So there will be another peak at 11 Hz, 19 Hz, 21 Hz, 29 Hz, and so on.

This observation is also true on the negative frequency axis. Just like the Fourier Transform, F(-f) = F*(f), where F is the FFT of a signal. So the absolute value of the FFT of a signal is symetric about the y-axis. So there will also be peaks at -1Hz, -9Hz, -11Hz, etc. This brings us to the fftshift function.

fftshift — Viewing Only the Unalaised Portion of a Spectrum

The introduction to the fftshift is given by the text at the end of the last cell. The purpose of the fftshift, is it shows the portion of the spectrum of a signal that is not aliased: -Fs/2 ≤ f ≤ Fs/2. An example is given below.

This new plot has the exact same spectral information as the other plots but without displaying aliased signals.

Getting signal amplitude and Phase Exactly Right

To get the FFT to provide the same transform as the Fourier Transform, the frequency of the signal must be exactly an integer multiple of df. This requires that N is an integer multiple of Fs.

Disambiguate Signals that are Close in Frequency

How can we tell the difference between two sinusoids that are close in frequency? This is disambiguation. The key here is to set df small enough to tell these signals apart. The only way to get df smaller is to observe a signal for a longer period of time. df = 1/(total observation time). Observe this fact below. It is unclear what signals are present in this spectrum.

Observe the change in the frequency spectrum as the total observation time is increased by increasing N.

Zero-Padding

Often we need to perform some operation in the frequency domain on two signals of different lengths. How can we do this? The answer is add zeros to the shorter of the two signals until they are the same length (padding), but what effect does zero-padding have on the spectrum of a signal? A demonstration is given below.

The simplest interpretation for the effect of zero padding is that it adds a kind of ripple to the spectrum of the signal. Usually not too much harm is done by zero-padding. A very common example where zero padding is useful is for calculating circular convolution. Another more trivial use is to provide more frequency resolution by decreasing df (since N is larger). An example of using zero-padding for increased frequency resolution is given below.

In the spectrum of the zero-padded signal a peak near 1.5 Hz is clearly visible. This is not true of the non-padded signal’s spectrum, which is roughly constant in that region.

Convolution with FFTs

One of the most useful features of the FFT results from its relationship with convolution. The relationship between the Fourier Transform and convolution is given as: f∗g = F^(-1)(F(f)∗F(g)), where F is the Fourier Transform. There is a similar relationship between the FFT and circular convolution. This is given as: f∘g = ifft(fft(f)×fft(g)). The difference between convolution and circular convolution is demonstrated below.

Circular convolution of f and g. Note that MATLAB’s cconv function provides more utility than the mathematical definition allows. By the mathematical definition of cconv, the output should be the length of the longest input (f in this case). MATLAB’s cconv can automatically zero-pad to make cconv equivalent to conv.

The Fast Fourier Transform may be used to perform circular convolution very quickly and efficiently. This functions due to the Discrete Fourier Transform property that convolution time domain is the same multiplication in the frequency domain. This is shown as follows.

Matched Filter

In many systems, signals are attenuated by Additive White Guassian Noise (AWGN). When searching for a specific signal in a noisy signal, the best means of finding that signal in the noise is via a matched filter. This is illustrated in this section.

The received signal is still observable in this case. This figure simply illustrates what is happening to the received signal when attenuated by the noise. The following figure illustrates the more realistic case, where the received signal cannot be visually seen.

We would like to eliminate the noise in the signal. To do this, we take advantage of the fact that the mean of the AWGN is 0. Therefore, if the noise could be summed, the mean should be nearly 0.

Consider what would happen if the noisey signal were convoled with a flipped version of the signal we are searching for (called the reference signal). Since convolution performs and flip and shift, the reference flipped reference signal is simply unflipped by convolution. Therefore, by flipping the reference signal before convolution, the process becomes merely a shifting of the reference signal over the noisy signal. This is referred to as a matched filter.

The principle idea of the matched filter being that the noise should be approximately summed to zero, but where the two signals overlap a maximum value in the convolution should be produced. Note that when implemented, the convolution is performed as circular convolution. This is visualized below.

The peak of the matched filter output should be located at nearly the amount that the reference signal needs to be shifted by for it to align with the signal buried in the noisy signal.

The matched filter is efficiently performed using the FFT as shown below. Also note that in general the complex conjugate of the yref must be taken; however, this is only necessary if the reference signal is complex.

Another desire arising from the matched filter’s output is to be able to shift the reference signal so that it aligns with the signal buried in the noisy signal. This is performed by using an Discrete Fourier Transform property known as the Shift Property (or sometimes Shift Theorem). This is shown in the next section.

Shifting 2D Signals with FFTs

This section demonstrates how signals may be shifted using the FFT. The shift theorem states that multiplication by a complex exponential in frequency domain results in a shift in the time domain.

Shifting 3D Signals

Shifting 3D functions is more difficult and requires several subtle concepts. One of the trickiest parts about shifting a 3D function in 2 dimensions is that structure of the matrix does not match the visualization. Observe in the plot that the point (0,0) on the xy-axis is in the lower left.

This is critical because the fx=0 point and fy=0 point in the shift must match the (0,0) point on the xy-axis. Therefore, the shift matrix must be constructed as follows. Furthermore note that .’ must be used to perform the transpose as opposed to ‘. In MATLAB .’ is a regular transpose. The ‘ operator is actual the conjugate transpose!