{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chebychev design of an FIR filter given a desired $H(\\omega)$\n",
    "\n",
    "A derivative work by Judson Wilson, 5/27/2014.<br>\n",
    "Adapted from the CVX example of the same name, by Almir Mutapcic, 2/2/2006.\n",
    "\n",
    "Topic References:\n",
    "\n",
    "* \"Filter design\" lecture notes (EE364) by S. Boyd\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This program designs an FIR filter, given a desired frequency response $H_\\mbox{des}(\\omega)$.\n",
    "The design is judged by the maximum absolute error (Chebychev norm).\n",
    "This is a convex problem (after sampling it can be formulated as an SOCP),\n",
    "which may be written in the form:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\max |H(\\omega) - H_\\mbox{des}(\\omega)| \n",
    "                             \\quad \\mbox{ for }  0 \\le \\omega \\le \\pi,\n",
    "    \\end{array}\n",
    "where the variable $H$ is the frequency response function, corresponding to an impulse response $h$.\n",
    "\n",
    "## Initialize problem data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cvxpy as cp\n",
    "\n",
    "#********************************************************************\n",
    "# Problem specs.\n",
    "#********************************************************************\n",
    "# Number of FIR coefficients (including the zeroth one).\n",
    "n = 20\n",
    "\n",
    "# Rule-of-thumb frequency discretization (Cheney's Approx. Theory book).\n",
    "m = 15*n\n",
    "w = np.linspace(0,np.pi,m)\n",
    "#********************************************************************\n",
    "# Construct the desired filter.\n",
    "#********************************************************************\n",
    "# Fractional delay.\n",
    "D = 8.25                # Delay value.\n",
    "Hdes = np.exp(-1j*D*w)  # Desired frequency response.\n",
    "\n",
    "# Gaussian filter with linear phase. (Uncomment lines below for this design.)\n",
    "#var = 0.05\n",
    "#Hdes = 1/(np.sqrt(2*np.pi*var)) * np.exp(-np.square(w-np.pi/2)/(2*var))\n",
    "#Hdes = np.multiply(Hdes, np.exp(-1j*n/2*w))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the minimax (Chebychev) design problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem status: optimal\n",
      "final objective value: 0.4999999999999996\n"
     ]
    }
   ],
   "source": [
    "# A is the matrix used to compute the frequency response\n",
    "# from a vector of filter coefficients:\n",
    "#     A[w,:] = [1 exp(-j*w) exp(-j*2*w) ... exp(-j*n*w)]\n",
    "A = np.exp( -1j * np.kron(w.reshape(-1, 1), np.arange(n)))\n",
    "\n",
    "# Presently CVXPY does not do complex-valued math, so the\n",
    "# problem must be formatted into a real-valued representation.\n",
    "\n",
    "# Split Hdes into a real part, and an imaginary part.\n",
    "Hdes_r = np.real(Hdes)\n",
    "Hdes_i = np.imag(Hdes)\n",
    "\n",
    "# Split A into a real part, and an imaginary part.\n",
    "A_R = np.real(A)\n",
    "A_I = np.imag(A)\n",
    "\n",
    "#\n",
    "# Optimal Chebyshev filter formulation.\n",
    "#\n",
    "\n",
    "# h is the (real) FIR coefficient vector, which we are solving for.\n",
    "h = cp.Variable(shape=n)\n",
    "# The objective is:\n",
    "#     minimize max(|A*h-Hdes|)\n",
    "# but modified into an equivelent form:\n",
    "#     minimize max( real(A*h-Hdes)^2 + imag(A*h-Hdes)^2 )\n",
    "# such that all computation is done in real quantities only.\n",
    "obj = cp.Minimize(\n",
    "    cp.max(  cp.square(A_R * h - Hdes_r)     # Real part.\n",
    "           + cp.square(A_I * h - Hdes_i) ) ) # Imaginary part.\n",
    "\n",
    "# Solve problem.\n",
    "prob = cp.Problem(obj)\n",
    "prob.solve()\n",
    "\n",
    "# Check if problem was successfully solved.\n",
    "print('Problem status: {}'.format(prob.status))\n",
    "if prob.status != cp.OPTIMAL:\n",
    "    raise Exception('CVXPY Error')\n",
    "print(\"final objective value: {}\".format(obj.value))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'FIR filter freq. response angle')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 9968 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 18376 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 30344 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plot inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)\n",
    "\n",
    "# Plot the FIR impulse reponse.\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.stem(range(n), h.value)\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('h(n)')\n",
    "plt.title('FIR filter impulse response')\n",
    "\n",
    "# Plot the frequency response.\n",
    "H = np.exp(-1j * np.kron(w.reshape(-1, 1), np.arange(n))).dot(h.value)\n",
    "plt.figure(figsize=(6, 6))\n",
    "# Magnitude\n",
    "plt.plot(w, 20 * np.log10(np.abs(H)),\n",
    "         label='optimized')\n",
    "plt.plot(w, 20 * np.log10(np.abs(Hdes)),'--',\n",
    "         label='desired')\n",
    "plt.xlabel(r'$\\omega$')\n",
    "plt.ylabel(r'$|H(\\omega)|$ in dB')\n",
    "plt.title('FIR filter freq. response magnitude')\n",
    "plt.xlim(0, np.pi)\n",
    "plt.ylim(-30, 10)\n",
    "plt.legend(loc='lower right')\n",
    "# Phase\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(w, np.angle(H))\n",
    "plt.xlim(0, np.pi)\n",
    "plt.ylim(-np.pi, np.pi)\n",
    "plt.xlabel(r'$\\omega$')\n",
    "plt.ylabel(r'$\\angle H(\\omega)$')\n",
    "plt.title('FIR filter freq. response angle')"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
