EE314 - Digital Signal Processing
Lab exercises

Session: 1

Related topics:

• Discrete-time sequences

1. Write a python script to generate $\delta [n]$ and plot it using stem method.

2. Write a python script to generate and plot $u[n]$.

3. Write a python script to generate and plot $A{\alpha }^n$ for (i) $\alpha \le -1$, (ii) , (iii) and (iv) $1\le \alpha$. Use subplots.

4. Write a python script to generate and plot $A\cos (\omega n+\varphi )$ for
   $\omega =0,\pi /4,\pi /2,3\pi /4,\pi ,5\pi /4,3\pi /2,7\pi /4,2\pi$. Use subplots.

Session: 2

Related topics:

• time-shifting, upsampling, downsampling, time-inverse

• Convolution

1. Write a python script to generate $x[an+b]$ from a given sequence $x[n]$ and user-provided input $a,b$.

   Verify your answer for $a=2,b=-3$

   $n$	$-3$	$-2$	$0$	$2$	$4$
   $x[n]$	$-1$	$2$	$0.5$	$-2.5$	$5$

2. Write a python script to perform convolution of two sequences in the conventional way, i.e., $x[n]\ast h[n]=\sum _{k}x[k]h[n-k]$

   Verify your answer for the $x[n]$ mentioned in the previous question and the following $h[n]$

   $n$	$0$	$2$	$3$	$4$	$5$
   $h[n]$	$-0.5$	$1$	$2$	$-3$	$1$

3. Write a python script to perform convolution using matrix multiplication. Represent the input sequence $x[n]$ as a column vector, $𝐱$, and multiply with the matrix $𝐇$, obtained from the impulse response $h[n]$, i.e., $𝐲=\mathrm{𝐇𝐱}$, where $𝐲$ is the column vector representation of the output sequence $y[n]$.

Session: 3

Related topics:

• Convolution

• Connection between DTFT and CTFS and DTFT and CTFT

Note: Please use line (continuous) plots while plotting continuous-time signals. Also, unlike discrete-time signals where index $n$ takes only integer values, the $t$ axis of continuous-time signals can have rational and irrational numbers.

1. Write a python script to verify that accumulator and forward difference are inverse of each other by forming matrices of appropriate sizes corresponding to those two systems.

2. Write a python script to verify (a) duality of DTFT and CTFS, (b) connection between DTFT and CTFT.

   For part (a), consider a periodic signal $x(t)=\cos (2\pi 100t)+0.5\cos (2\pi 400t)$ and find its Fourier series coefficients. Form a discrete-time signal from those coefficients while ensuring that the indices of coefficients remain same in the discrete-time sequence. Now calculate the DTFT of the sequence you obtained and visually compare it with the plot of $x(t)$ over a suitable interval and list your observations.

   For part (b), consider a continuous time signal $x(t)=\exp (-|t|)$. First, find its Fourier transform and then find the Fourier transform of $x(t)$ after sampling it with a pulse train. Form a discrete-time signal from those samples and calculate its DTFT. Visually, compare the Fourier transform of the sampled signal and the DTFT, and list your observations.

3. Write a python script to demonstrate Gibb’s phenomenon for discrete-time idea low-pass filter.

   Consider the impulse response of an ideal lowpass filter with cut-off frequency ${\omega }_c$,

   Now to demonstrate Gibb’s Phenomenon, plot $H_{lp,M}(e^{j\omega })$ for different values of $M$. For example, $M=1,3,5,7,13,17,19$, where

   $$H_{lp,M}=\sum _{n=-M}^M{h}_{lp}[n]e^{-j\omega n}.$$

Session: 4

Related topics:

• z-transform

1. Plot $X(z)$ for the following sequences and then extract the values corresponding to $|z|=1$ to plot on a linear scale and then compare with the DTFT of the corresponding $x[n]$

   (i) $x=\{1,2,3,2,1\}$ where $n=\{0,1,2,3,4\}$

   (ii) $x=\{1,2,3,2,1\}$ where $n=\{-3,-2,-1,0,1\}$

   (iii) $x[n]={0.5}^{n}u[n]$

Session: 5

Related topics:

• z-transform

1. Write a python script to perform convolution by multiplying two polynomials in $z^{-1}$. Verify your answer for the sequences mentioned in the Lab sessions: 2 and 3.

3. Write a python script to plot poles and zeros of a rational transfer function on the $z$-plane. Show your plots for the following transfer functions:

   1. $$H(z)=\frac{z^{-1}+z^{-2}}{1-\frac{1}{6}{z}^{-1}-\frac{1}{6}{z}^{-2}}.$$

   2. $$H(z)=\frac{3-7z^{-1}+5z^{-2}}{1-\frac{5}{2}{z}^{-1}+z^{-2}}.$$

   3. $$H(z)=\frac{1}{{(1+z^{-1}/2)}^2(1-2z^{-1})(1-3z^{-1})}.$$

Session: 6

Related topics:

• Block diagram representation

• DFT

1. Generate the $N\times N$ DFT matrix discussed in the class. Multiply the sequence $x[n]$ given in Lab session: 2 to find its DFT. Compare your output with the numpy’s fft function’s output.

3. Simulate the following systems on Simulink by creating their block diagram representation

   1. Moving average filter:

      $$y[n]=\frac{1}{N}\sum _{k=0}^{N-1}x[n-k],n,k\in \mathbb{Z}.$$

   2. Auto-regression filter:

      $$y[n]=\sum _{i=1}^M{\alpha }_{i}y[n-i]+\beta x[n],n,i\in \mathbb{Z};{\alpha }_i,\beta \in ℝ.$$

   3. Auto-regression moving average (ARMA):

      $$y[n]=\sum _{i=1}^M{\alpha }_{k}y[n-i]+\sum _{k=0}^{N-1}{\beta }_{k}x[n-k],n,i,k\in \mathbb{Z};{\alpha }_k,{\beta }_k\in ℝ.$$

Session: 7

Related topics:

• DFT, IDFT

• Circular & Linear convolution

• Circulant matrices

1. Write a python script that performs $N$ point circular convolution of two sequences. Note that $N$ must be equal to or larger than the length of both sequences. Use this script to perform linear convolution of two given sequences. Consider the sequences given in Lab sessions 2 and 3.

3. Perform linear convolution by first calculating the DFT of two given sequences and then inverse DFT of the element-wise product of resultant sequences, i.e.,

   $$x[n]\ast h[n]=\mathrm{IDFT}\{\mathrm{DFT}\{x[n]\}.\mathrm{DFT}\{h[n]\}\}.$$

   Verify your answer by comparing it with the results obtained from other methods you have learned in the earlier sessions.

5. Verify that the inverse DFT matrix is indeed the eigenvector matrix of a circulant matrix and the eigenvalues are the DFT of the sequence from which the circulant matrix is formed. Generate a random sequence of length $N=5,10,20$ and the corresponding inverse DFT matrices for this task.

Session: 8

Related topics:

• FFT

• CTFT, DTFT

1. Write a python script that finds $N$ point DFT of a given sequence using the decimation-in-time algorithm and find the DFT of the following sequences

   $x[n]=$ $\{3,-1,-7,-4,2\}$ (1)

   $y[n]=$ $\{6,-1,2,-3,-7,5,2,-9,4,1,0,-3\}$ (2)

   Note that you need to pad zeros to both sequences before applying $8$ and $16$ point FFT algorithms, respectively.

   Cross-check your answers against the output of numpy.fft.fft(). Find the norm of the differences of the sequences returned by your algorithm and the numpy’s fft algorithm. It should be a very small quantity, in the range of ${10}^{-14}$ or lower, if both the sequences are identical.

3. Repeat the above tasks for decimation-in-frequency algorithm.

5. Redo the Task 2(b) of Session#3. Find the frequency that encompasses 95% of the spectrum of the given signal and use that to find the Nyquist rate.

Session:9

Related topics:

• Sampling

• Discrete-time processing of continuous-time signals

• Continuous-time processing of discrete-time signals

1. Consider the signal mentioned in Tasks 2(b) of Session#3, i.e., $x(t)=\exp (-|t|)$. In the previous lab session, you have identified the Nyquist rate for this signal (albeit after truncating it to a band-limited signal). Now, reconstruct the signal from its samples taken at the Nyquist rate by adding appropriately time-shifted and amplitude-scaled $\mathrm{sinc}$ functions, i.e.,

   $$x_r(t)=\sum _{n}x[n]\frac{\sin (\pi (t-nT)/T)}{\pi (t-nT)/T},$$

   where $T$ is the sampling period chosen according to the Nyquist-Shannon theorem to avoid aliasing.

3. Write a python code to perform low-pass filtering of a continuous-time signal using a discrete-time low-pass filter. Either use the signal mentioned in the Task 1 or a sum sinusoids with three frequencies $1000,1500,3000$.

   At first, find the cut-off frequency of the continuous-time low-pass filter. Then, calculate the cut-off frequency of the discrete-time low-pass filter. Convolve the samples with the impulse response of the discrete-time low-pass filter. Reconstruct the continuous-time signal from the samples obtained after convolution. Compare the continuous-time signals (input and output) to verify whether the filtering is done correctly.

5. Change the cut-off frequency of the (effective) continuous-time low-pass filter by changing the sampling interval and apply on the sum of sinusoids mentioned in the previous task.

Session: 10

Related topics:

• Sampling, Resampling

• Audio signal processing

Note: You need to install sox to execute the commands mentioned in the tasks listed below. The installation file (Windows) is included in the directory where codes are located.

1. Reading & plotting of a .wav file

   $\star$ Python: Download the file “readPlayPlotAudio.py” from the shared cloud directory and modify it to the read audio files in the “WavFiles” cloud directory.

   $\star$ MATLAB: Download the file “readPlayPlotAudio.m” from the shared cloud directory and modify it to the read audio files in the “WavFiles” cloud directory.

   $\star$ Sox: Read the info of the file using

   >sox.exe --info <filename>

   (Sox doesn’t have any plotting or playing functionality. You may use sounder.exe, included in the same directory that contains MATLAB, Python codes, to play an audio file. The command to play an audio file is sounder.exe <file_name>).

3. Synthesizing an audio file – single and multitone signals

   Generate synthetic audio files and combine them to create a song. Use the following frequencies to generate synthetic audio files:

   Notation	Frequency (Hz)	Note
   $f_a$	220	a
   $f_b$	246.94	b
   $f_c$	523.25	c
   $f_d$	585.93	d
   $f_e$	657.68	e

   Table 1: Table of frequencies for different notes.

   Then combine them in the following order: [a, b, c, d, e]. Repeat the Task 1 for the generated audio file.

   $\star$ Python: Download the file “syntheticAudioFileGeneration.py” from the shared cloud directory and modify it to accomplish the task mentioned above.

   $\star$ MATLAB: Download the file “syntheticAudioGen.m” from the shared cloud directory and modify it to accomplish the task mentioned above.

   $\star$ Sox: Use the following commands to synthesize and combile audio files:

   sox -n -b 16 output.wav synth 2.25 sine 300 // generates a 2.25 second sine wave at 300Hz with 16-bit precision

   sox file1.wav file2.wav combined.wav // concatenate two files to create a single file

5. Processing audio files – filtering, resampling

   Write Python/MATLAB codes to resample, filter the shared audio files and repeat the Task 1 to observe changes in the processed file.

   $\star$ Use the following commands to resample and filter audio files:

   sox input.wav output.wav lowpass 100.0 // Apply a 100 Hz low pass filter on the input file to generate the output file.

   sox input.wav -r 22050 -b 8 output.wav // Change the sampling rate to 22050 and the precision to 8-bit

Session: 11

Related topics:

• Window sequences

• Spectral analysis using DFT

• Short-time Fourier transform

1. Plot the following windows in the time domain and their frequency spectrum in $-\pi$ to $\pi$ range using MATLAB’s wvtool, and document your observations

   (a) rectangular, (b) Bartlett, (c) Hamming, (d) Hann, (e) Blackman, (f) Kaiser

3. Spectrum analysis using DFT

   Write a python script that estimates the frequency spectrum of a continuous time signal from its samples by first splitting them into smaller blocks by windowing and then calculating their DFT.

   Consider different types of windows mentioned in the Task 1 and observe their effects on the estimated spectrum. Pad zeros to blocks you obtained earlier and find the DFT of the zero padded blocks. Does zero padding help in getting a better estimation?

   Test your script on the following signals:

   (a) $x(t)=e^{-|t|}$ (refer to the Task 1 of Session: 9)

   (b) $x(t)=\cos (600\pi t)+\sin (3000\pi t)$

5. Time-dependent Fourier transform or Short-time Fourier transform (STFT)

   Write a python script to find the STFT of the signals mentioned in the Task 2. You have to compute the DFT of the shifted windowed signal to find the samples of the DTFT and plot them in a 3D plot.

   In addition to the two signals mentioned earlier, find the STFT of $\cos (2\pi t^2)$ and plot it with respect to time index ($n$) sampled frequency ($2\pi k/N$).

Session: 12

Related topics:

• Frequency response of discrete-time systems

1. Plot pole-zero plots, Magnitude and phase response of the following systems using MATLAB^™’s freqz(), zplane() functions and comment on what kind of filter (lowpass, highpass, bandpass, bandstop etc.) they represent.

   1. Difference equation:

      $$y[n]=x[n]-x[n-N]$$

      Use $N=1,5,10$

   2. Moving average filter:

      $$y[n]=\frac{1}{N}\sum _{k=0}^{N-1}x[n-k]$$

      Use $N=3,5,10$

   3. Auto-regression filter:

      $$y[n]=\sum _{i=1}^M{\alpha }_{i}y[n-i]+\beta x[n]$$

      Consider the following values for the parameters

      $M$	${\alpha }_i$	$\beta$
      1	0.5	1
      1	0.1	5
      3	(1/3, 1/3, 1/3)	1
      3	(1/2, 1/4, 1/8)	1

   4. Auto-regression moving average (ARMA):

      $$y[n]=\sum _{i=1}^M{\alpha }_{i}y[n-i]+\sum _{k=0}^{N-1}{\beta }_{k}x[n-k]$$

      Consider the following values for the parameters

      $(M,N)$	${\alpha }_i$	${\beta }_k$
      (1,3)	0.5	(1,1,1)
      (3,3)	(1/3, 1/3, 1/3)	(1,1/2,1/4)
      (3,3)	(1/3, 1/3, 1/3)	(1/3, 1/3, 1/3)
      (3,3)	(1/2, 1/4, 1/8)	(1,1/2,1/4)

Session: 13

Related topics:

• Frequency response of LTI system

• Filter design

1. Plot the normalised frequency response of the following systems (review the “firstorder” and “secondorder” MATLAB files and make necessary changes )

   (a) A single pole/zero located at $0.3\exp (j\pi /6)$

   (b) A pair of complex-conjugate poles/zeros located at $0.3\exp (j\pi /6)$ and $0.3\exp (-j\pi /6)$

   (c) An allpass system with complex-conjugate poles $0.8\exp (j\pi /4)$ and $0.8\exp (-j\pi /4)$

3. Minimum, mixed and Maximum phase systems

   Review the MATLAB file with name “minimumphase” and modify it to observe the changes in the phase response due to the positioning of the two zeros of the transfer function. What happens when (i) both the zeros are inside the unit circle, (ii) both are outside the unit circle, (iii) one is inside and the other one is outside. Note that the provided MATLAB file has codes corresponding to a single zero system.

5. Filter design

   (a) Explore MATLAB’s filterDesigner app and design a Butterworth, a Chebyshev-I, a Chebyshev-II, a Elliptic, an FIR lowpass filter with the following specification

   $0.89125\le |H(e^{j\omega })|\le 1,$ $0\le \omega \le 0.2\pi$

   $|H(e^{j\omega })|\le 0.17783,$ $0.3\pi \le \omega \le \pi$

   Note down your observations about the magnitude, phase response, filter orders of the filters designed with the app.

   (b) Use the MATLAB functions butter, buttord, cheby1, cheby1ord, cheby2, cheby2ord, ellip, ellipord, fir1, fir2 to design a filter with the specifications mentioned in the part (a).

References

• Alan V. Oppenheim, Ronald W. Schafer, John R. Buck, “Discrete-time signal processing” 2nd Ed, Prentice Hall

• John G. Proakis, Dimitris G. Manolakis., “Digital Signal Processing – Principles, Algorithms, and Applications” Fourth Edition, Pearson 2007.

• Sounder https://elifulkerson.com/projects/commandline-wav-player.php

• https://realpython.com/python-wav-files/

• https://docs.python.org/3/library/wave.html#module-wave

• https://docs.scipy.org/doc/scipy/reference/generated/scipy.io.wavfile.read.html

• https://linguistics.berkeley.edu/plab/guestwiki/index.php?title=Sox_in_phonetic_research

• http://class.ece.iastate.edu/mmina/ee186/labs/audio.htm

• https://ccrma.stanford.edu/~jos/JFB/JFB_2up.pdf Slide nos. 3 and 4, Source: CCRMA,Stanford University

• https://www.cs.cornell.edu/courses/cs6220/2017fa/CS6220_Lecture3.pdf Lecture notes, Cornell University

• https://in.mathworks.com/help/signal/ref/convmtx.html MATLAB convolution matrix