Useful links
Goal
Construct a variety of filters using zeros and poles. Apply notch filter to audio sample to remove a single-frequency hum.
Instructions
See Lab #1 for a longer narrative about lab mechanics. The 1-sentence version:
Complete the tasks below, submit your task files on-line before the deadline, and complete your check-off interview within a week of the file submission deadline.
As always, help is available in 32-083 (see the Lab Hours web page).
Introduction
You may have noticed that ocassionally inexpensive audio equipment generates a humming sound. This humming is produced by the 60Hz AC electrical power "leaking" into the audio signal processing. This hum is at such low frequency that usually it is not annoying. On most aircraft, the electrical power is at 400Hz and when 400Hz "leaks" into audio signals, the result is very annoying. In this problem you will try to eliminate the 400Hz hum from a digital signal, but of course you will try to accomplish this without eliminating much of the desired signal.
An example of a signal with and without the 400Hz hum are in the files with_hum.wav and without_hum.wav, respectively. These are sampled versions of the audio signal, where the audio was sampled 8000 times a second. You can load the files into numpy arrays by using
with_hum = lab6.read_sound('with_hum.wav')
without_hum = lab6.read_sound('without_hum.wav')
To hear the sound (use headphones if you at a public terminal!!!), you can use lab6.play_sound, passing it a numpy array of audio samples. For example
lab6.play_sound(with_hum)
Note: The appropriate Python sound libraries are installed in the course Athena locker. If you're working on your own computer, you will need to download and install the pyaudio module using the link above.
In this lab you will be developing functions that will allow you to design a notch filter, and you can then apply the notch filter to eliminate the annoying hum from the audio signal.
Task 1: Determine magnitude of frequency response (1 point)
In this task we'll compute magnitude of frequency response for our two channels. To compute the frequency response of a channel we'll first compute the unit-sample response for the channel. You can use your code for Lab #2 Task #3 assuming it's written to accept an arbitrary channel as input -- if it's not, please write a version of unit_sample_response that does take an arbitrary channel function as input and use it for all the tasks here in Lab #6.
Now that we have the h's for the unit-sample response, compute the magnitude of the frequency response using in the formula
|H(ejΩ)| = |Σm h[m]e-jΩm| for m=0 to M
where M is the number of samples in the unit-sample response. Write a function freq_res_usr starting from our template in lab6_1.py. The function should take a numpy array with the unit-sample response and return two numpy arrays giving the test angles and the magnitude of the frequency response at those angles. In our template, we used 401 frequency points in the range 0 ≤ Ω &le π.
One can use complex numbers in Python just as one uses ints and floats. To create a complex constant with real part r and imagininary part i, use complex(r,i) or (r + i*1j). The built-in arithmetic operators and functions will accept complex operands, as will numpy's array operators. For example, if you have an array u with complex elements, the python statement
v = numpy.exp(u)
will create a numpy array v whose length is the same as the length of u, containing the values v[0]=exp(u[0]), v[1]=exp(u[1]), etc. Note also that numpy.abs computes the magnitude of a complex number. lab6_1.py is the template file for this task:
# template file for Lab #6, Task #1
import numpy
import matplotlib.pyplot as p
import lab6
from lab2_3 import unit_sample_response
def freq_res_usr(usr, num_freqs=401):
# 0 <= omega <= pi
omega_list = numpy.linspace(0,numpy.pi,num=num_freqs)
# create and fill in a numpy array hejw containing
# the complex frequency response at each omega in
# omega_list
pass # your code here...
return omega_list,numpy.abs(hejw)
def freq_res(channel, num_freqs=401):
usr = unit_sample_response(channel) # Get channel's u.s.r.
return freq_res_usr(usr,num_freqs)
if __name__ == '__main__':
# generate plots of freq_2_usr for both channels
# using both of the methods
omega_1, mag_hejw_chan1 = freq_res(lab6.channel1)
lab6.plot_freq_response(omega_1, mag_hejw_chan1,'channel 1')
omega_2, mag_hejw_chan2 = freq_res(lab6.channel2)
lab6.plot_freq_response(omega_2, mag_hejw_chan2,'channel 2')
p.show()
Task 2: Derive unit-sample response for a filter using zeros (2 points)
Write a function zfreqs_to_usr that takes a list of frequencies, and computes a unit-sample response for a filter whose frequency response at each of the list of frequencies is zero. You will find this program much easier to write if you create filters that zero each individual frequency and convolve together their unit-sample responses using numpy.convolve.
You should scale the final unit-sample response so that the maximum amplitude of the frequency response is 1. Given a unit-sample response h, the maximum amplitude of the frequency response can be determined empirically: max(freq_res_usr(h)[1]).
We have given you a template with a simple test case in the file lab6_2.py.
# template file for Lab #6, Task #2
import numpy
import matplotlib.pyplot as p
import lab6
from lab6_1 import freq_res_usr
def zfreqs_to_usr(omega_list):
pass # your code here...
if __name__ == '__main__':
# A quick test for zfreqs_to_usr
omega_list = [numpy.pi/20.0, numpy.pi/5.0, numpy.pi/2.5]
usr = zfreqs_to_usr(omega_list)
# Plot the frequency response of the result
omega, mag_hejw = freq_res_usr(usr)
lab6.plot_freq_response(omega, mag_hejw,'three-zero filter')
p.show()
Task 3: Filter design (1 point)
Using your code from Task #2, write a Python function low_pass that takes a cutoff frequency Ωco and returns the appropriately scaled unit-sample response for a three-zero filter where the magnitude of the frequency response is
Similarly, write a function high_pass that does the same thing except the magnitude of the frequency response should be
The region where the magnitude is supposed to be zero is called the stop band of the filter; the region where the magnitude is one is called the pass band. We'll measure the quality of your filters using the following ratio:
(minimum magnitude in pass band) (maximum magnitude in stop band)
When computing the min and max values we'll omit points within ±0.1π of the cut-off frequency to leave room for a transition region in your frequency response. Design your functions to place the zeros to maximize this ratio.
lab6_3.py calls your functions, plotting their frequency response and printing out the ratio above.
# template file for Lab #6, Task #3
import numpy
import matplotlib.pyplot as p
import lab6
from lab6_1 import freq_res_usr
from lab6_2 import zfreqs_to_usr
# return appropriately scaled unit-sample response
# for a low-pass filter with specified cutoff frequency
def low_pass(cutoff_freq):
pass # your code here
# return appropriately scaled unit-sample response
# for a high-pass filter with specified cutoff frequency
def high_pass(cutoff_freq):
pass # your code here
if __name__ == '__main__':
# A low pass filter, cut-off of 0.6pi
lp_cutoff = 0.6 * numpy.pi
usr_lp = low_pass(lp_cutoff)
# Plot the response of the result
omega_lp, mag_hejw_lp = freq_res_usr(usr_lp)
lab6.plot_freq_response(omega_lp,mag_hejw_lp,
'low-pass filter @ $\Omega=%3.2g$' % lp_cutoff)
# Determine the performance ratios for lowpass
quality_lp = \
lab6.filter_quality((lp_cutoff+0.1*numpy.pi,numpy.pi),
(0,lp_cutoff-0.1*numpy.pi),
omega_lp,mag_hejw_lp)
print 'quality ratio for low-pass filter=',quality_lp
# A high pass filter, cut-off of 0.4pi
hp_cutoff = 0.4 * numpy.pi
usr_hp = high_pass(hp_cutoff)
# Plot the response of the result
omega_hp, mag_hejw_hp = freq_res_usr(usr_hp)
lab6.plot_freq_response(omega_hp,mag_hejw_hp,
'high-pass filter @ $\Omega=%3.2g$' % hp_cutoff)
quality_hp = \
lab6.filter_quality((0,hp_cutoff-0.1*numpy.pi),
(hp_cutoff+0.1*numpy.pi,numpy.pi),
omega_hp,mag_hejw_hp)
print 'quality ratio for high-pass filter=',quality_hp
p.show()
Task 4: Bandpass and notch filters (1 point)
Write a Python function bandpass that takes two frequencies Ωlo and Ωhi and returns the appropriately scaled unit-sample response for a filter where the magnitude of the frequency response is
A series combination of appropriately-chosen low-pass and high-pass filters will do the job. Use your low_pass and high_pass functions from Task #3 and then convolve their responses to get the response for the series combination.
Write a Python function notch that takes two frequencies Ωlo and Ωhi and returns the appropriately scaled unit-sample response for a filter where the magnitude of the frequency response is
A parallel combination of appropriately-chosen low-pass and high-pass filters will do the job. Use your low_pass and high_pass functions from Task #3 and then sum their responses to get the response for the parallel combination.
lab6_4.py is a template file that plots the response of your two filters.
# template file for Lab #6, Task #4
import numpy
import matplotlib.pyplot as p
import lab6
from lab6_1 import freq_res_usr
from lab6_3 import low_pass,high_pass
# return appropriately scaled unit-sample response
# for a band-pass filter with specified cutoff frequencies
def bandpass(lo,hi):
pass # your code here
# return appropriately scaled unit-sample response
# for a notch filter with specified cutoff frequencies
def notch(lo,hi):
pass # your code here
if __name__ == '__main__':
omega_lo = 0.4 * numpy.pi
omega_hi = 0.6 * numpy.pi
# A band-pass filter
usr_band = bandpass(omega_lo,omega_hi)
# Plot the response of the result
omega_band, mag_hejw_band = freq_res_usr(usr_band)
lab6.plot_freq_response(omega_band,mag_hejw_band,
'band-pass filter @ $\Omega=(%3.2g,%3.2g)$' % (omega_lo,omega_hi))
# A notch filter
usr_notch = notch(omega_lo,omega_hi)
# Plot the response of the result
omega_notch, mag_hejw_notch = freq_res_usr(usr_notch)
lab6.plot_freq_response(omega_notch,mag_hejw_notch,
'notch filter @ $\Omega=(%3.2g,%3.2g)$' % (omega_lo,omega_hi))
p.show()
Task 5: Single-frequency notch filter (1 point)
Write a Python function notch_one_freq that returns the unit-sample response for a single-frequency notch filter using a single complex conjugate pair of zeros. The function you wrote in Task #2 should be helpful.
Using notch_one_freq write a Python function eliminate_hum that accepts a numpy array of audio samples, passes the samples through your single-frequency notch filter and returns the result. You'll have pick the appropriate frequency for the notch filter. In a discrete-time system, Hz and Ω are related by the formula
&Omega = 2π*freq_in_Hz/samples_per_second
lab6_5.py is template file for this task that plots the frequency response of your filter, reads in the wav file with the hum, processes the samples with your function and then plays the result:
# template file for Lab #6, Task #5
import numpy
import matplotlib.pyplot as p
import lab6
from lab6_1 import freq_res_usr
from lab6_2 import zfreqs_to_usr
# return unit-sample response for a single-frequency notch
# filter at the specified frequency
def notch_one_freq(freq):
pass # your code here
# Using notch_one_freq, create a notch filter at the frequency
# of the hum (you'll have to convert 400Hz to the appropriate
# omega), and then apply the filter to a numpy array of audio
# samples.
def eliminate_hum(audio_samples):
pass # your code here
if __name__ == '__main__':
# plot freq response of notch_one_freq @ pi/2
omega_notch, mag_hejw_notch = \
freq_res_usr(notch_one_freq(numpy.pi/2))
lab6.plot_freq_response(omega_notch,mag_hejw_notch,
'notch filter @ $\Omega=\pi/2$')
print "Please close plot window to start audio playback."
p.show()
# Read in the good sound and play
without_hum = lab6.read_sound('without_hum.wav')
lab6.play_sound(without_hum)
# Read in the corrupted sound and play
with_hum = lab6.read_sound('with_hum.wav')
lab6.play_sound(with_hum)
# Filter the sound and play
hum_removed = eliminate_hum(with_hum)
lab6.play_sound(hum_removed)
Task 6: Derive unit-sample response for a filter with one pole (2 points)
Write a Python function pole_to_coeffs that takes the pole frequency and a scale factor r (r < 1). Your function should create a complex conjugate pair of poles and and then use the poles to compute the coefficients (a0, a1, a2) for the difference equation
a0*y[n] + a1*y[n-1] + a2*y[n-2] = x[n]
so that the frequency response of the resulting linear system has the specified poles. Please refer to course annotated notes for a correction to the definition of a pole. The polynomial associated with the difference equation above is
a0 + a1*z + a2*z2 = (z - 1/p0)*(z - 1/p1)
= (z - 1/(r*ejΩ))*(z - 1/(r*e-jΩ))
= (1/r2) - (2/r)*cos(&Omega)*z + z2
where Ω is the pole frequency. You should scale the values for a0, a1, and a2 so that the maximum magnitude of the frequency response at the pole location is 1. The appropriate scaling factor can be determined by evaluting frequency response of your single-pole filter at the pole frequency (clearly the maximum value of the response for this single-pole filter) and taking its magnitude.
lab6_6.py includes a function pole_to_usr that computes the unit-sample response of your difference equation. The testing code creates a filter with a pole at π/2 and plots its frequency response.
import numpy
import matplotlib.pyplot as p
import lab6
from lab2_3 import unit_sample_response
from lab6_1 import freq_res_usr
# compute coefficients for a difference equation with
# a pair of complex conjugate poles at the specified
# frequency and with a scale factor of r.
def pole_to_coeffs(pole_freq,r):
pass # your code here
# use the differential equation solver to convert
# diffeq coefficients into a unit-sample response
def pole_to_usr(pole_freq,r):
ycoeffs = pole_to_coeffs(pole_freq,r)
xcoeffs = [1.0] + [0.0]*(len(ycoeffs) - 1)
channel = lab6.de_solver(ycoeffs,xcoeffs)
return unit_sample_response(channel)
if __name__ == '__main__':
# get unit-sample response for filter with a pole at pi/2
pole_freq = numpy.pi/2
r = 0.98
usr_pole = pole_to_usr(pole_freq,r)
# plot the frequency response
omega_pole, mag_hejw_pole = freq_res_usr(usr_pole)
lab6.plot_freq_response(omega_pole,mag_hejw_pole,
'pole filter @ $\Omega=%3.2g$, r=%3.2g' % (pole_freq,r))
p.show()
Task 7: Improved hum elimination (2 points)
Using your function from Task #6, add a single pole to your implementation of notch_one_freq in Task #5 -- if you choose a good location for the pole, the frequency response should be much more like an ideal notch (see lecture notes). You can add the pole filter in series with the filter from Task #5, convolving their unit-sample responses and rescaling to get the unit-sample response for the series combination.
lab6_7.py is a template file for this task:
# template file for Lab #6, Task #7
import numpy
import matplotlib.pyplot as p
import lab6
from lab6_1 import freq_res_usr
from lab6_2 import zfreqs_to_usr
from lab6_6 import pole_to_usr
# NEW! Now improved by the addition of a pole!
# return unit-sample response for a single-frequency notch
# filter at the specified frequency
def notch_one_freq(freq):
pass # your code here
# Using notch_one_freq, create a notch filter at the frequency
# of the hum (you'll have to convert 400Hz to the appropriate
# omega), and then apply the filter to a numpy array of audio
# samples.
def eliminate_hum(audio_samples):
pass # your code here
if __name__ == '__main__':
# plot freq response of notch_one_freq @ pi/2
omega_notch, mag_hejw_notch = \
freq_res_usr(notch_one_freq(numpy.pi/2))
lab6.plot_freq_response(omega_notch,mag_hejw_notch,
'notch filter @ $\Omega=\pi/2$')
print "Please close plot window to start audio playback."
p.show()
# Read in the good sound and play
without_hum = lab6.read_sound('without_hum.wav')
lab6.play_sound(without_hum)
# Read in the corrupted sound and play
with_hum = lab6.read_sound('with_hum.wav')
lab6.play_sound(with_hum)
# Filter the sound and play
hum_removed = eliminate_hum(with_hum)
lab6.play_sound(hum_removed)