ECE 445S Real-Time Digital Signal Processing Lab - Homework 2 Hints

• Please keep all of the work for each problem together. This makes it much easier for the grader to give you the appropriate full or partial credit due for the solution.

• Problem 2.1.

  Please see the homework #1 hints and solution for problem 1.1.

  The connection between the z-domain and discrete-time frequency domain is z = exp(j w) where w is the discrete-time frequency in rad/sample. If we were to draw the shape of exp(j w) for -pi < w ≤ pi, we would have a circle of radius 1, which is the unit circle. In the z-domain, the discrete-time frequency w is the angle with respect to the positive real axis. Here are the equations and a visualization of this connection from lecture on Feb. 8, 2021.

  For the magnitude response, the angles of poles near but inside the unit circle indicate the passband of the filter, and the angles of the zeros near or on the unit circle indicate the stopband:

  • Lecture 6 slides 6-8 through 6-11.
  • Chapter 5 in Lathi's Linear Systems and Signals book
  • Section 10.5 in Oppenheim & Willsky's Signals and Systems book
  • Chapter 8 in McClellan, Schafer and Yoder's Signal Processing First

  Poles on the unit circle correspond to an impulse response that is pure oscillation. This is from lab #2 and homework problem 0.4 (see homework #0 solution).

  Useful Matlab commands include

  • zplane plots the poles and zeros of a transfer function in the z-domain.
  • freqz plots the magnitude and phase response for a given transfer function.

  The frequency response of a discrete-time linear time-invariant (LTI) system can be found a number of ways:

  • If the region of convergence includes the unit circle, where the unit circle is z = exp(j w) for real-valued frequency w in rad/sample, then one can compute the frequency response from the transfer function in the z-domain:

    Hfreq(w) = H(z) |
                     z = exp(j w)
    

  • Otherwise, one could compute the discrete-time Fourier transform of the impulse response to find the frequency response (see lecture slide 5-15). Using the discrete-time Fourier transform may require you to search for the answer in books, e.g. Linear Systems and Signals by B. P. Lathi (chapter 9) or Signals and Systems by A. V. Oppenheim and A. S. Willsky (section 5.6). Recall that the discrete-time Fourier domain is periodic with period 2 pi.

  Bounded-input bounded-output (BIBO) stability means that whenever any signal that is bounded (finite) in amplitude for all time is input into an LTI system, the output signal will always have bounded amplitude? See lecture 6 slides 6-13 and 6-14, and Handout on BIBO Stability:

  • All FIR filters have BIBO stability.
  • Causal IIR filters must have all poles inside the unit circle.
  • General result for all LTI systems: region of convergence must include the unit circle for LTI system to be BIBO stable

  Part (c). The discrete-time integrator, a.k.a. running summation, has an infinite impulse response. Because the discrete-time integrator is not BIBO stable, the discrete-time integrator is not an infinite impulse response (IIR) filter. Nonetheless, it is quite useful.

  Part (d). When plotting the magnitude response of the bandpass resonator, please use a large value for r such as 0.9 or 0.95 and use a value of w₀ between pi/4 and 3 pi/4.

• Problem 2.2:

  Please download the audio file twosignals.wav. If reading the wave file gives an error, then there is a Matlab function called wavchunksizefix to fix it.

  Part (a). Please see the homework hints and solution for problem 1.2 concerning spectrogram plots. Here's a set of parameters for the spectrogram that will reveal the two signals, including the frequencies in which they overlap:

  [twosignals, fs] = audioread('twosignals.wav');   
  nwin = 512;       % divide signal into blocks of nwin samples with each
  noverlap = 511;   % block overlapping by noverlap samples with previous
  nfft = 512;       % number of frequency points to calc Fourier transform
  figure;
  spectrogram( twosignals, nwin, noverlap, nfft, fs, 'yaxis' );
  

  Parts (b) and (c). For LTI systems, frequencies on the output had to be present on input, a.k.a. LTI systems do not create new frequencies (see lecture slides 3-9, 3-21, and 5-15). In the discrete-time frequency domain, Y(w) = H(w) X(w) where X(w) and Y(w) are the input and output discrete-time Fourier transforms, respectively and H(w) is the frequency response of the system. You'll be designing filters to produce a desired frequency response H(w) to pass certain frequencies and attenuate others.

  One way to design a finite impulse response (FIR) filter using the Parks-McClellan algorithm (Remez exchange algorithm) is to use the MATLAB function firpm, which is introduced in JSK on pages 47 and 48. (Octave users should see the footnote on page viii of JSK.) Here is an example template of code to design a lowpass filter:

  %% Set the Nyquist frequency to be half of the sampling rate fs.
  %% The sampling rate fs is set when reading the twosignals file 
  %% using the audioread command, and will be used by the sound 
  %% command to send audio samples to the sound card for playback.
  fnyquist = fs/2;
  %% Define the passband frequency fpass in Hz
  fpass = ??;
  %% Define the stopband frequency fstop in Hz
  fstop = ??;
  ctfrequencies = [0 fpass fstop fnyquist];
  idealAmplitudes = [1 1 0 0 ];
  pmfrequencies = ctfrequencies / fnyquist;
  %% Number of coefficients is filter order plus one
  filterOrder = ??;
  filterCoeffs = firpm( filterOrder, pmfrequencies, idealAmplitudes );
  

  There are many possible values for the passband frequency, stopband frequency and filter order. Please fill in your choices in the above code template.

  For Matlab versions 2016a and 2016b, the implementation of the Parks-McClellan FIR filter design algorithm was failing to converge for filter orders of 50 and higher. There are a couple of workarounds. First, try to see if there is a filter order lower than 50 that is able to filter out the bird chirping and keep the gong sound. Second, you could design the FIR filter using the filter design and analysis tool (filterDesigner) and export the FIR coefficients for use in a Matlab script.

  Once the FIR filter coefficients have been computed, you could filter the signal by using the MATLAB command conv or filter. The conv function would only work for FIR filtering, whereas the filter function can work for both FIR and IIR filters. Please see the handout Four Ways to Filter a Signal for more information.

  You could analyze the resulting signal using the plotspec function from the MATLAB files that accompany the Software Receiver Design book.

  Part (d). Downsampling is the removal of samples in a regular fashion. In downsampling by 2, the first sample and every other sample thereafter is kept, and hence, the second sample and every other sample thereafter is removed. Downsampling by 2 would halve the number of samples.

  Play back the downsampled filtered signal at the same sampling rate as the original gong file. The Matlab code to downsample a vector vec by 2 follows:

  vecDownsampledBy2 = vec(1:2:length(vec));
  

  The MATLAB code 1:2:N generates the vector [1 3 5 7 .... ].

  Part (e). Upsampling is the insertion of samples in a regular fashion. After every sample in a signal, upsampling by 2 would insert a sample of value zero. Upsampling by 2 would double the number of samples, which would either double the sampling rate or double all of the frequencies if the sampling rate were kept the same, depending on how the upsampled signal is converted to continuous time. Please play the upsampled audio at twice the sampling rate of the twosignals audio signal. Please note that upsampling might lead to a form of aliasing.

  After midterm #1, we'll make heavy use of upsampling in a communications transmitter and downsampling in a communications receiver.

  Here's Matlab code for upampling by 2 of a signal whose amplitude samples are stored in a row vector vec:

  upsampledLength = 2*length(vec);
  vecUpsampledBy2 = zeros(1,upsampledLength);
  vecUpsampledBy2(1:2:upsampledLength) = vec;
  

  You could define a signal to try out the above code, e.g.

  vec = cumsum( ones(1,10) );          %%% [1 2 3 4 ... 10 ]
  

• Problem 2.3:

  The design of a bandstop filter to be applied to streaming music and music files to lessen tinnitus loudness was the subject of Problem 3 of Midterm #1 in Fall 2016.

  Passband Ripple. In the Matlab filter design and analysis tool (filterDesigner), a passband ripple of 1 dB would mean that the passband magnitude response in deciBels (dB) would be over -0.5 dB to 0.5 dB for FIR filters and over -1 dB to 0 dB for IIR filters. The Matlab filterDesigner will plot these interpretations if you select "Filter Specifications" under the "Analysis" menu. It is also available under the Matlab Introduction to Digital Filter Design help page.

  Validating Designs. Be sure to validate that the FIR filter designed by the filter design and analysis tool (filterDesigner) meets the specification. First, you will need to display the magnitude response in filterDesigner by choosing "Magnitude Response" under the "Analysis" menu or the appropriate icon shortcut. Then, force the magnitude response plot to show all values by selecting "Full View" under the "View" menu. Finally, you will need to zoom into the passbands and stopband by using the zoom tool (available as the icon of the magnifying glass with a plus sign in the middle of it). You can hit CONTROL-Z to undo the most recent zoom.

  Here's a handout on verifying FIR filter specifications by Dan Jacobellis (TA).

  For additional validation of the design, you could export a filter design from the filter design and analysis tool (filterDesigner) to the Matlab workspace by selecting "Export" under the "File" menu in filterDesigner. Then, you can then plot the magnitude and phase responses by using freqz(impulseResponse) where impulseResponse is the Matlab variable name that holds the filter coefficients (the default variable name is Num).

  Part (a). For a demonstration of using the filter design and analysis tool (filterDesigner) in Matlab to design FIR filters, please start viewing the spring 2014 video for lecture 6, part 5, at the 7:45 mark.

  The FIR Equiripple Design Method is also known as the Remez Exchange or Parks-McClellan filter design algorithms. The Parks-McClellan algorithm is an iterative algorithm that attempts to fit a Chebyshev polynomial into the desired magnitude filter design specification to determine the FIR filter coefficients.

  Please select the minimum order design.

  For the FIR Equiripple Design Method, the estimation of the minimum filter order has an accuracy of +/- 10%. So, finding the lowest filter order requires designing an FIR filter using the equiripple method for many different filter orders. As a user of filterDesigner, we can enter filter specifications that are different than the design specifications and then evaluate the designed filter with respect to the original design specifications. We can force the minimum filter order design algorithm to use a lower filter order by reducing the stopband attenuation. For example, try a stopband attenuation of 78 dB.

  For FIR least squares design, you might consider using the Matlab command firls to obtain the filter coefficients instead of using filterDesigner. The usage is very similar to firpm (remez) you used in problem 2.2. Once you've designed the filter, you can use freqz to plot the the magnitude response. The advantage of this method is that you can tweak the filter order and re-run simulations until the filter specification is met. The FIR Least Squares method commonly misses the stopband specification.

  In certain versions of filterDesigner, the specification for the Kaiser window method uses cutoff frequencies instead of passband and stopband frequencies. The cutoff frequency is the value of the magnitude response that half of the passband magnitude response, or equivalently 6 dB less than the passband magnitude response. You can estimate the cutoff frequency as the average frequency in transition band.

  For the Kaiser window method in filterDesigner, you could use the baseline filter order obtained from the Parks-McClellan algorithm as a reference. Instead of playing with the filter order itself, I would recommend changing the filter parameters entered in the design tool to see if a lower filter order can be obtained. Instead of making the filter parameters more difficult to meet, one could decrease the stopband attenuation and increase the passband ripple.

  You can validate a filter designed in filterDesigner by zooming into the magnitude response to check the magnitude response in dB at the stopband and passband frequencies, as well as at DC and pi.

  Part (b). In filterDesigner, you can plot the impulse response by choosing the appropriate option under the Analysis menu. You can also export the impulse response by using the Export option under the File menu to export the impulse response to the Matlab workspace in the variable Num (by default). In Matlab, you could plot the impulse response using the stem command.

  FIR filters whose impulse response is either even or odd symmetric about its midpoint have linear phase. Realizable IIR filters cannot have linear phase. See lecture slides 5-19, 5-20, and 5-21.

  Part (d). For an FIR filter of N coefficients with input signal x[n] and output signal y[n] is defined in the discrete-time domain as

  y[n] = h0 x[n] + h1 x[n-1] + ... + hN-1 x[n - (N-1)] 
  

  for n ≥ 0.
  1. Computational complexity. To compute one output sample, we can manually count the number of multiplications needed. There are N terms being added together, and each term requires one multiplication. Computing one output would need N multiplications and N-1 additions.
  2. Memory storage. Storage includes the N-length vector of current and N-1 previous input values, an N-length vector of coefficients, and y[n], which totals 2N+1 words. A word means the size of the data type needed to store a value. On a 32-bit processor, the size of a single-precision floating point number is 32 bits (4 bytes), which consists of an 8-bit exponent (with offset) and a 24-bit mantissa (with a sign bit). The DSP processor that we used to use in the lab component had a 32-bit single-precision floating point format (float data type in C) and a 64-bit double-precision floating-point format (double data type in C).

  For FIR filters, circular buffers are far more efficient at managing the buffer of the current input and N-1 previous input values than using a linear array. Updating a circular buffer when a new input value arrives takes 2 reads and 2 writes, whereas updating a linear array takes N reads and N writes.

bevans@ece.utexas.edu