Instructions | Deliverables | Walkthrough | References

Detailed Instructions (Recommended Steps in Python)
  1. Every plotting step in this lab will plot exactly the same five things: a signal, the magnitude of its DTFT, the phase of its DTFT, a magnitude model, and a phase model. Each of these should be plotted on appropriate axes, with an appropriate title, and with a zero-line to show you when the thing being plotted is negative. In order to save time, create a function that does those plotting steps:
    import cmath
    import matplotlib.pyplot as plt
    def plot_signal_and_model(signal,omega,dtft,ampmodel,phasemodel):
    zero_line = np.zeros(len(omega))
    plt.subplot(5,1,1)
    plt.stem(signal)
    plt.title('Signal')
    plt.subplot(5,1,2)
    magnitude_response = [ abs(x) for x in dtft ]
    plt.plot(omega,magnitude_response)
    plt.title('Magnitude of DTFT')
    plt.subplot(5,1,3)
    phase_response = [ cmath.phase(x) for x in dtft ]
    plt.plot(omega,phase_response,omega,zero_line)
    plt.title('Phase of DTFT')
    plt.subplot(5,1,4)
    plt.plot(omega,ampmodel,omega,zero_line)
    plt.title('Amplitude Model')
    plt.subplot(5,1,5)
    plt.plot(omega,phasemodel,omega,zero_line)
    plt.title('Phase Model')
  2. The first signal we're going to work with is a length-11 rectangle. First, define the length of the analysis window, L, and the length of the rectangle, N. We want the analysis window to be about ten times the rectangle length, so that we have lots of extra frequency samples in the DTFT; this trick is called zero-padding. Try something like:
    L = 100
    N = 11
     In python, you can create a length-N list by first creating a length-one list, then multiplying the list object by N. Similarly, you can concatenate two lists by using the + symbol. Use these tricks to create a rectangle signal that contains N ones, followed by (L-N) zeros:
    x_signal = [1]*N + [0]*(L-N)
     Take its DTFT:
    import scipy.fftpack
    X_dtft = scipy.fftpack.fft(x_signal)
     The amplitude of X_dtft is abs(X_dtft). Similarly, the phase of each term is computed by [ cmath.phase(X) for X in X_dtft ]. The thing is, we know in advance what these should be, because we learned them in class. The phase should be -omega*N/2. The amplitude should be sin(omega*N/2)/sin(omega/2), except at omega=0, where this is undefined; at omega=0, it should be N. Create a function that computes an omega axis, an amplitude model, and a phase model:
    import math
    def dtftmodel(N,L):
    w_axis = [ 2*math.pi*k/L for k in range(0,L) ]
    a = [N]+[ math.sin(w*N/2)/math.sin(w/2) for w in w_axis[1:L] ]
    t = [ -w*(N-1)/2 for w in w_axis]
    return (w_axis, a, t)     Now plot your signal and its model, and save the result:
    (omega_axis, ampmodelX, phasemodelX) = dtftmodel(N,L)
    plt.figure(1)
    plot_signal_and_model(x_signal,omega_axis,X_dtft,ampmodelX,phasemodelX)
    plt.savefig('lab7fig1.png')
  3. Notice that the phase model is perfectly straight, while the true phase is jagged. There are two reasons for this, which we'll now explore one after the other. First, the cmath.phase function returns only the principal value of the phase, that is, it always returns a number between -pi and pi. Create a function wrap_phase that wraps the phase, so that it is between -pi and pi:
    def wrap_phase(phi):
    theta=phi.copy()
    for k in range(0,len(theta)):
    while theta[k] < -math.pi:
    theta[k]+=2*math.pi
    while theta[k] > math.pi:
    theta[k] -= 2*math.pi
    return theta
     The phi.copy() step is important, because there will be several times when we'll want to keep both the wrapped and unwrapped phase models. For now, try applying this function to your phase model, and plot the result in lab7fig2.png:
    phasemodel2 = wrap_phase(phasemodel)
    plt.figure(2)
    plot_signal_and_model(x_signal,omega_axis,X_dtft,ampmodelX,phasemodel2)
    plt.savefig('lab7fig2.png')
  4. OK, wrapping the phase helped, but the true phase still has twice as many discontinuities as the model. On a related note: the amplitude model is signed, whereas the true amplitude of the DTFT is always positive. The reason for this oddity: (-1) is equal to exp(j*pi). Every time the amplitude model goes negative, it has the effect of adding pi to the phase model. Change your wrap_phase function to reflect this:
    def wrap_phase2(A,phi):
    theta=phi.copy()
    B = A.copy()
    for k in range(0,len(theta)):
    if B[k] < 0:
    theta[k] += math.pi
    B[k] = abs(B[k])
    while theta[k] < -math.pi:
    theta[k]+=2*math.pi
    while theta[k] > math.pi:
    theta[k] -= 2*math.pi
    return theta
     Now try normalizing the phase and amplitude models again, and plot the result as lab7fig3.png:
    (ampmodel3,phasemodel3) = wrap_phase2(ampmodelX,phasemodelX)
    plt.figure(3)
    plot_signal_and_model(x_signal,omega_axis,X_dtft,ampmodel3,phasemodel3)
    plt.savefig('lab7fig3.png')
  5. Now let's create a different signal --- let's call it h_signal --- which is a rectangle of length 5. Find its amplitude model and phase model, wrap the phase, plot, and save the result as lab7fig4.png:
    M = 5
    h_signal = [1]*M + [0]*(L-M)
    H_dtft = scipy.fftpack.fft(h_signal)
    (omega_axis, ampmodelH, phasemodelH) = dtftmodel(M,L)
    (ampmodel4,phasemodel4)=wrap_phase2(ampmodelH,phasemodelH)
    plt.figure(4)
    plot_signal_and_model(h_signal,omega_axis,H_dtft,ampmodelH,phasemodelH)
    plt.savefig('lab7fig4.png')
     You should find that the spectrum and its model are identical. If not, take a moment to debug.
  6. Remember the main reason the DTFT was invented? It's because convolution of two signals is the same as multiplication of their DTFTs. Convolve x and h, truncate the result to a length of just L samples, and then take its DTFT:
    y_signal = np.convolve(x_signal,h_signal)
    y_signal = y_signal[0:L]
    Y_dtft = scipy.fftpack.fft(y_signal)
     Now create the model. The amplitude of Y should be H times X; the phase of Y should be phase(H)+phase(X):
    phasemodelY=[ phasemodelH[k]+phasemodelX[k] for k in range(0,L) ]
    ampmodelY=[ ampmodelH[k]*ampmodelX[k] for k in range(0,L) ]
    (ampmodel5,phasemodel5)=wrap_phase2(ampmodelY,phasemodelY)
    plt.figure(5)
    plot_signal_and_model(y_signal,omega_axis,Y_dtft,ampmodel5,phasemodel5)
    plt.savefig('lab7fig5.png')

Deliverables (Required)

By 3/14/2017 23:59, upload to compass a zip file named MYNAME_LAB5.ZIP containing the following things:

  1. Five images, from the five parts listed above.
  2. A program that generates your figures.

Walkthrough

Here's a video walkthrough of the lab.

References