EE341 Digital Filtering and Audio Equalizer Solved

Part 1: Digital Filters
To start, create a new Jupyter Notebook in your ee341lab folder and name it to lab4. Import the modules: numpy, pyplot, fftpack, and simpleaudio just as you did in Lab 3. Also import the signal module from SciPy library as follows.

from scipy import signal 
The signal module provides several convenient functions for creating filters and applying them to signals. It will be our main tool for implementing filters in this lab.

Then, download microsoft_stock.txt and music.wav from Canvas and copy them to the same directory where your Jupyter Notebook is stored. We will use the signals from these files to test our filters.  

Background
Every casual LTI digital filter can be described by a linear difference equation that relates the output of the filter 𝑦[𝑛] to the input of the filter 𝑥[𝑛]:
                                                                              𝑁                                                    𝑀

∑ 𝑎𝑘𝑦[𝑛 − 𝑘] = ∑ 𝑏𝑘𝑥[𝑛 − 𝑘]

                                                                              𝑘=0                                              𝑘=0

The time-domain difference equation presented above is equivalent to the frequency-domain transfer function:

 𝐻(𝑧) = 𝑋𝑌((𝑧𝑧)) = ∑∑𝑁𝑀𝑘=0𝑎𝑏𝑘𝑘𝑧𝑧−−𝑘𝑘               𝑧              𝑗ω  𝑤ℎ𝑒𝑟𝑒      = 𝑒

𝑘=0

In the equation, 𝑎𝑘 and 𝑏𝑘 are the filter coefficients. 𝑎𝑘 is an array of length 𝑁 + 1 and 𝑏𝑘 is an array of length 𝑀 + 1. The length of the filter, also called the number of taps, is the maximum of the lengths of the two coefficient arrays, i.e. max(𝑁, 𝑀) + 1. The order of the filter is the length of the filter minus 1, i.e. max(N, M).

By rearranging the equation, we can isolate the current output 𝑦[𝑛] to the left-hand side:
                                                                                           𝑀                                                  𝑁

𝑦  (∑ 𝑏𝑘𝑥[𝑛 − 𝑘] − ∑ 𝑎𝑘𝑦[𝑛 − 𝑘])

𝑎0
                                                                                         𝑘=0                                             𝑘=1

This gives us the formula for calculating the current output of the filter from the current input, past inputs and past outputs of the filter. The leading denominator coefficient 𝑎0 is usually set to 1, since the effect of a non-unit 𝑎0 is equivalent to dividing all other coefficients by 𝑎0.

 In Python, LTI digital filters can be implemented using the lfilter function from the signal module:

y = signal.lfilter(b, a, x) 
In the example code above, x is the input signal, y is the output signal, and a and b are the filter coefficients. All of these 4 are stored as 1D ndarrays.

Types of Filters
Filters can be classified by the shape of their frequency responses, which determine their applications.

•  Low-pass filter removes high-frequency components of the data, allowing only frequencies lower than the cut-off frequency to pass through. It is commonly used for smoothing the data. For audio, applying low-pass filter removes high-pitched sound and upper harmonics, which results in a warm and muddy sound. For image, applying low-pass filter removes sharp details, which results in a blurred image.  

• High-pass filter removes low-frequency components of the data, allowing only frequencies higher than the cut-off frequency to pass through. It is commonly used for removing DC component from the signal and for extracting details from the data. For audio, applying high-pass filter removes the bass component, which results in a sharp and crispy sound.

• Band-pass filter allows only a certain range of frequency to pass through. It is usually used in communication systems, such as radio, to select a channel. Multiple band-pass filters can also form an equalizer which offers detailed control of the level of each band.

•  Band-reject filter removes a certain range of frequency. It is usually used to remove a specific unwanted frequency from the signal, e.g. removing the feedback of a loudspeaker.

• All-pass filter lets all frequency components equally pass through. However, it applies different amount of phase-shift (delay) to different frequency, which is useful for some applications.

FIR Filter
LTI digital filters can also be classified into two types by structure: finite impulse response (FIR) filter, and infinite impulse response (IIR) filter. A FIR filter is a filter with no denominator coefficients 𝑎𝑘, except for 𝑎0, which we set to 1. Therefore, by definition, the output of a FIR filter can be expressed as:

𝑀

𝑦[𝑛] = ∑ 𝑏𝑘𝑥[𝑛 − 𝑘]

𝑘=0

Shown by the equation, the output of a FIR filter is simply a weighted average of the current and past inputs, with the weights being the nominator coefficients 𝑏𝑘. For this reason, people also call it as the moving average filter. Also notice that the equation exactly describes the convolution operation: 𝑦[𝑛] = 𝑏[𝑛] ⊛ 𝑥[𝑛], with the impulse response being the nominator coefficients 𝑏𝑘. Because there can only be a finite number of nominator coefficients, the length of the impulse response is also finite; hence the name “finite impulse response”.

 Since the impulse response has finite length, you can implement the filter by directly computing the convolution instead of using signal.lfilter. To calculate convolution, use signal.convolve as follows:

y = signal.convolve(x, h) 
where h is the impulse response of the filter, stored as a 1D ndarray. The advantage of convolve is that it can leverage the circular convolution theorem to use FFT to speed-up the calculation.

Task 1: Implement a Moving Average Filter
Load the data in signal_data.txt to an ndarray using the following statement:

data = np.genfromtxt('signal_data.txt') 
The data is a 100-second sinusoidal signal with random noise. Plot the data. Notice the quick changes in the data which correspond to the high frequency component. People tend to remove the noise by using a moving average filter. Consider the following 30-second moving average filter:

𝑦 

which corresponds to the impulse response:

ℎ[𝑛] =   (𝑢[𝑛] − 𝑢[𝑛 − 30])

Apply this filter to the 100-second signal using both lfilter and convolve. For each the two implementations, plot the original data and the filtered data in the same plot. Does the filtered result look smoother? Plot the magnitude of the frequency response of the filter to verify that this system corresponds to a low-pass filter. 

In your notebook, include the plot for the lfilter implementation, the plot for the convolve implementation and the plot of the system frequency response. Be sure to label the axes and give each plot a meaningful title. Explain the differences in the results of the two implementations. How would you adjust the implementations so that they give identical results? 

IIR Filter
We already know that FIR filters don’t have any denominator coefficients 𝑎𝑘 for 𝑘 ≥ 1. IIR filters, in contrast, have some non-trivial denominator coefficients. Therefore, by definition, the output of an IIR filter can be expressed as:

                                                                                        𝑀                                                  𝑁

𝑦[𝑛] = ∑ 𝑏𝑘𝑥[𝑛 − 𝑘] − ∑ 𝑎𝑘𝑦[𝑛 − 𝑘]

                                                                                      𝑘=0                                             𝑘=1

As shown by the equation, the output of an IIR filter is the weighted average of not only the current and past inputs, but also past outputs. Since the current value of 𝑦 depends on the previous values of 𝑦, IIR filter is also known as the recursive filter. Because of the recursion, the impulse response of a stable IIR filter decays over time but never dies off, making it infinite length, and hence the name “infinite impulse response”.

Apart from the length of the impulse response, FIR and IIR have many other different properties as listed below.   Order: FIR filters require a large filter order to meet the requirements of passband, stopband, ripple, and roll-off. In contrast, IIR filters can accomplish them with a lower order. Therefore, FIR filters are computationally more expensive compared to IIR filters.  

•   Frequency Response: The frequency response of a FIR filter has ripples in both passband and stopband and has a slow roll-off at the cut-off frequency. In contrast, the frequency response of an IIR filter is much smoother and has sharper roll-off at the cut-off frequency.

•   Stableness: The poles of a FIR filter can only be at the origin, so it must be stable. In contrast, the poles of an IIR filter can be anywhere, so it can be either stable or unstable.

•   Phase: FIR filters can be designed to have a linear phase which preserves the shape of the waveform. In contrast, IIR filters have non-linear components in the phase response that can distort the shape of the waveform.

Using Classic Filter Designs
Instead of designing filters by ourselves, we can use some classic filter designs provided by the SciPy package.

For IIR filter, we can use the Butterworth design, which is a design that has no ripple in the frequency response.

To create a Butterworth filter, use the signal.butter function:[1]

b, a = signal.butter(order, W, type) 
where b and a are the filter coefficients, order is the order of the filter, freq is the cutoff frequency and type can be either ’lowpass’ or ’highpass’. The unit of W is the normalized frequency, which has the range 0 ≤ 𝑊 ≤ 1 where 1 corresponds to the Nyquist frequency (half of the sampling rate). This range also maps to the DTFT or DFT range: 0 ≤ ω ≤ π or 0 ≤ 𝑓 ≤ 0.5 we used in Lab 3.

For FIR filter, we can use the window-method design, which approximates the impulse response of an ideal filter. To create a window-method filter, use the signal.firwin function:

# Lowpass: 

b = signal.firwin(num_taps, W) 
# Highpass: 

b = signal.firwin(num_taps, W, pass_zero = False) 
where b is the numerator coefficient or impulse response, num_taps is the length of the filter and W is the cutoff frequency which has the unit of the normalized frequency.

Task 2: Analyze the Characteristics of Filters
For this task, you will use the analyze function provided in Appendix A to analyze the characteristics of 4 different filters. In particular, we will look at how FIR and IIR filters behave differently in terms of frequency response, impulse response and poles / zeros.

To start, copy the analyze function in Appendix A to a new code cell and execute it. Then, implement the 4

following filters. Each filter should have an order of 10 and a cutoff frequency of ω = 0.1𝜋.

•       Low-pass IIR filter using Butterworth design.

•       Low-Pass FIR filter using window-method design.

•       High-pass IIR filter using Butterworth design.

•       High-pass IIR filter using window-method design.

Then for each filter, apply the analyze function, which will produce plots of the frequency response, impulse response and poles / zeros.

Finally, apply each filter to the each of the following signals:

• The data in microsoft_stock.txt, which is the daily stock price of Microsoft over 4 years.

•   A pulse of length 20:

𝑥[𝑛] = 𝑢[𝑛] − 𝑢[𝑛 − 20] where 𝑥[𝑛] is of total length 60 (so append 40 zeros to the end).

In your notebook, include the 4 plots produced by the analyze function, and for each input signal, plot the original signal and the output of each filter in the same set of axes. Comment on the differences in the frequency response of the two filters (magnitude and phase) and how this impacts the outputs. When commenting on frequency responses, consider how close the filter matches an ideal filter.

Part II: Audio Equalizer
In this section, you will implement an audio equalizer using three different filters. Equalizer works by first breaking down the audio signal into multiple frequency bands using filters. Each band will then receive a different gain. Finally, all bands are summed together to produce the output signal. Audio mixing boards do exactly this, though with many more bands than the three bands we have here. The term “equalizer” comes from the fact that people often want to boost the gain of low-energy frequency ranges of the sound or to reduce the gain of high-energy frequency ranges of the sound. By doing so, each frequency range of the sound will have about equal loudness, allowing them to mix together nicely.

Shown below is the signal flow diagram of the 3-band equalizer you will implement. 𝐺1,𝐺2 and 𝐺3 are correspondingly the gains of the bass, mid and treble frequency ranges. The coefficients of each filter are provided in Table 1. However, you need to figure out by yourself which filter does each set of coefficients correspond to. The magnitude response of each filter is shown in Figure 2.  

Figure 1: Signal flow of the equalizer 

Table 1: Filter coefficients of the equalizer 

Filter Coefficients
Which filter? (LP, HP or BP?)
b1 = [11.713, 0, -23.4257, 0, 11.713] a1 = [1800, -6656.7, 9341, -5896.1, 1413.4] 
 
b2 = [0.2982, 0.89471, 0.89471, 0.2982] a2 = [1800, -4989.8, 4625.2, -1433] 
 
b3 = [688.1, -2752.5, 4128.71, -2752.47, 688.12] a3 = [1000, -3256.6, 4033.77, -2246, 473.51] 
 
Task 3: Implement the 3-Band Equalizer
First, identify which filter does each set of coefficients in Table 1 correspond to. Then, implement the equalizer by writing a function that takes a sound signal and the gain terms 𝐺1,𝐺2 and 𝐺3 as inputs, applies the filters, multiplies the filter outputs by the gain terms, and sums the results. The gain terms supplied to the function should be in decibels (dB).

After you implemented your equalizer, load the sound file music.wav using the function provided in Appendix B (the same one you used in Lab 3). Apply the equalizer to the music with 𝐺1 = 𝐺2 = 𝐺3 = 0 𝑑𝐵 and play the output. Verify that it sounds the same as the original input. Then, experiment with different sets of gains (for example 𝐺1 = 𝐺2 = 0 𝑑𝐵 and 𝐺3 = −40 𝑑𝐵, or 𝐺1 = −40 𝑑𝐵 and 𝐺2 = 𝐺3 = 0 𝑑𝐵). Discuss how the filtered sound sounds different in at least two different cases.  

To play sound, use the function below.

def play(samples, sample_rate): 
  sa.play_buffer((np.clip(samples, -1, 1) * 32767).astype('int16'),                  1, 2, sample_rate).wait_done() 
In your notebook, include your answers to the blanks in Table 1, the code of your equalizer function and the code of playing the outputs of the equalizer. Discuss how the filtered sound sounds different under different sets of gains.
.