mp573_hw1_handout.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework Set 1: Signals and Systems, Fourier Transforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## MP 573 HW1\n",
    "## Import some useful stuff\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, Image\n",
    "import matplotlib.image as mpimg\n",
    "from os.path import dirname, join as pjoin\n",
    "import scipy.io as sio\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1: Visualizing various signals\n",
    "\n",
    "## A. Let's start by visualizing complex exponentials \n",
    "\n",
    "Create a complex exponential signalwith frequencies $u_0 = 1/2$ and $v_0 = 1/5$ along x and y, respectively. Try displaying our signal f using the command  imshow. You should get an error - what is the problem? \n",
    "\n",
    "Next, display the real part, imaginary part, magnitude, and phase of  f. Show or sketch these four images in your homework solutions (OK to show your completed jupyter notebook). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate 2D complex sinusoidal signal, and plot its real/imaginary values, and also magnitude/phase\n",
    "N = 501\n",
    "xmax = 4\n",
    "x = np.linspace(-xmax,xmax,N)\n",
    "y = x;\n",
    "X, Y = np.meshgrid(x,y)\n",
    "u0 = 1/2 # This is the frequency along x\n",
    "v0 = 1/5 # This is the frequency along y\n",
    "f = np.exp(1j*(2*np.pi*(u0*X + v0*Y)))\n",
    "\n",
    "\n",
    "## To-Do: Display f directly (should probably give an error)\n",
    "\n",
    "\n",
    "\n",
    "## To-Do: Display the real and imaginary parts, magnitude and phase of f\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## B. Add some initial phase $\\pi/2$\n",
    "\n",
    "\n",
    "Repeat part A, but this time with an initial phase $\\theta = \\pi /2$ (this is the phase at x=0 and y=0). How do the four displayed images (real part, imaginary part, magnitude, and phase) change with respect to part A, if at all? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## To-Do: Add initial phase\n",
    "\n",
    "\n",
    "\n",
    "## To-Do: Display the real and imaginary parts, magnitude and phase \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## C. Try different frequencies\n",
    "Repeat part A, but this time with frequencies  $u_0 = 1/4$ and $v_0 = 3$. How do the four displayed images (real part, imaginary part, magnitude, and phase) change with respect to part A, if at all?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## To-Do: Simulate 2D complex sinusoidal  f\n",
    "\n",
    "\n",
    "\n",
    "## To-Do: Display the real and imaginary parts, magnitude and phase of f\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Support of convolution\n",
    "\n",
    "## 2.2 Practice\n",
    "\n",
    "### A. \n",
    "Verify your result by calculating numerically (eg: in Matlab) the convolution of the following signals: $f_1(x) =\\hbox{rect}(x/2)$ and $f_2(x) =\\hbox{rect}(2x)$. Note that in most computational environments there is no exact notion of a signal $f(x)$, only data arrays. However, you can define an array of x values (x = np.linspace(-5,5,1001)), and then define a discretized rect signal relative to your x values (f1 = 1.0*(np.abs(x)<1), and analogously for $f_2(x)$). Then you can convolve the two arrays to obtain an array $g$, but make sure to keep track of the dimensions of the result and the meaning of these dimensions (eg: consider using the command \\texttt{convolve} with an optional parameter \" mode='same' \".  Plot the result on the appropriate x-scale and include the plot in your homework solutions. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## To-do: Simulate 1D input signal f1(x) and f2(x)\n",
    "N = 1001\n",
    "xmax = 5\n",
    "x = np.linspace(-xmax,xmax,N)\n",
    "dx = x[1]-x[0]\n",
    "f1 = 1.0*(np.abs(x)<1)\n",
    "f2 = 1.0*(np.abs(x)<1/4)\n",
    "\n",
    "## To-Do: calculate the convolution between f1 and f2\n",
    "\n",
    "# To-Do: Plot the two signals f1 and f2, and their convolution g\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### B. \n",
    "Does the support (size of the non-zero region) in your computational result $g$ match your analytical prediction from the previous question (note that these should match, at least up to the discretization level on your arrays)? Does the height of your computational result $g$ (ie: the actual values of $g$) match your analytical prediction (note that this will generally not happen unless you've taken steps to make them match)?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### C. \n",
    "Note that the (discrete) convolution implemented in a package such as Python/scipy is actually a summation rather than an integral: \n",
    "$$g[n] = \\sum_m f_1[m] f_2[n-m]$$ \n",
    "Repeat the operation above (create f1 and f2, and convolve them to obtain g), but this time based on a finer grid of size $4001$ (x = np.linspace(-5,5,4001)) and then based on an even finer grid of size $8001$ (x = np.linspace(-5,5,8001)). Does the height of your convolution result depend on the resolution of your grid? How would you normalize your  calculation so that the height of the resulting convolution height does not depend on the resolution of the grid, and matches the analytical calculation? (Think about approximating the convolution integral defined in class using a summation). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## To-do: Simulate 1D input signal f1(x) and f2(x)\n",
    "N = 1001\n",
    "xmax = 5\n",
    "x = np.linspace(-xmax,xmax,N)\n",
    "Delta_x = x[1]-x[0]\n",
    "f1 = 1.0*(np.abs(x)<1)\n",
    "f2 = 1.0*(np.abs(x)<1/4)\n",
    "\n",
    "## To-Do: calculate the convolution between f1 and f2 \n",
    "g = Delta_x*signal.convolve(f1, f2, mode='same') # Include dx factor to normalize the sum (approximate the integral)\n",
    "\n",
    "# To-Do: Plot the two signals f1 and f2, and their convolution g\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 Separable convolution\n",
    "\n",
    "## 3.3 Practice\n",
    "\n",
    "Use Python (or some other computational tool) to generate a 2D rect function (as defined above) over an (x,y) grid with range (-3,3) along each dimension. \n",
    "\n",
    "Next, calculate the convolution of rect(x,y) with itself, both as a 2D convolution (Python: use the command signal.convolve2d), as well as a product of two 1D convolutions (Python: signal.convolve, which you can apply separately to $f_1$ and $f_2$). Make sure you understand the  optional parameters and outputs of these commands as they handle finite arrays rather than signals on a continuous domain. For instance, consider using a third parameter mode='same' (including the single quotes) when calling the convolve or convolve2d function. Also, note that you cannot just feed 2D arrays to the 1D convolution function convolve and expect that it will know the dimension along which you want to convolve. \n",
    "\n",
    "Show that either approach (single 2D convolution, or product of two 1D convolutions) leads to the same resulting 2D signal. You can show this by calculating the relative difference between the two resulting signals (2D arrays), which we can name g1 and g2. You can calculate this relative difference as RelDiff = norm(g1(:)-g2(:))/norm(g1(:)), which should be a very small number, possibly not exactly zero due to numerical precision limitations. \n",
    "\n",
    "Does your empirically calculated convolution have the same shape as the analytically calculated convolution from the previous question? You can assess this by plotting a profile (say, a horizontal profile through the center of your array, ie: for y=0) and comparing to your resulting expression from the previous question. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate 2D rect signal\n",
    "N = 101\n",
    "xmax = 3\n",
    "x = np.linspace(-xmax,xmax,N)\n",
    "y = x;\n",
    "f1 = 1.0*(np.abs(x)<0.5)\n",
    "f2 = 1.0*(np.abs(y)<0.5)\n",
    "# Make sure to understand how fx and fy are created, \n",
    "# including the logical operators ``abs(x)<0.5'' and ``abs(y)<0.5'' \n",
    "f = np.outer(f1,f2)\n",
    "# Make sure to understand the previus line, \n",
    "# including the outer product of two vectors\n",
    "\n",
    "\n",
    "# To-Do: 2D convolution\n",
    "\n",
    "# To-Do: 1D convolutions, then combine\n",
    "\n",
    "# To-Do: Calculate the difference between g1 and g2\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}