{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "56bf4898-66b4-4e46-8bee-e23bf3fc66bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.family\"] = \"serif\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e0ebd31-d7be-44db-9e15-30aea93fdd55",
   "metadata": {},
   "source": [
    "Make a transfer function (lambda)\n",
    "\n",
    "Make a weighting function\n",
    "\n",
    "Make a delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "359ac0ff-873a-49f0-9db5-948839bf461f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 37748 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example frequency range (logarithmic scale)\n",
    "frequencies = np.logspace(-1, 2, 10000)  # Frequencies from 10 to 1000 Hz\n",
    "\n",
    "# Example transfer function H(s) = 1 / (s + 1) evaluated at jw (jω)\n",
    "# H(jω) = 1 / (jω + 1)\n",
    "# s = jw where w = 2 * pi * f\n",
    "s = 1j * 2 * np.pi * frequencies  # jω = j * 2π * frequency\n",
    "\n",
    "# Calculate the transfer function response\n",
    "P = -s*(1/10)/((1/100)*s+1)/((1/10)*s+1)\n",
    "W_2 = 0.21*s/(0.1*s + 1)\n",
    "delta = 1/np.cos(s/1e3)\n",
    "P_tilde = (1+delta*W_2)*P\n",
    "P_env = (1+W_2)*P \n",
    "\n",
    "# Calculate magnitude and phase\n",
    "magnitude = np.abs(P)\n",
    "phase = np.angle(P)\n",
    "\n",
    "magnitude_2 = np.abs(P_tilde)\n",
    "phase_2 = np.angle(P_tilde)\n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, 20 * np.log10(magnitude))  # Plot magnitude in dB\n",
    "ax1.set_title('Bode Plot - Continuous')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylabel('Magnitude (dB)')\n",
    "ax1.set_ylim([-60,60])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, np.degrees(phase))  # Plot phase in degrees\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylabel('Phase (Degrees)')\n",
    "ax2.set_ylim([-180,180])\n",
    "ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('nom_cont.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bdc7fd3d-053a-4b56-974b-b74a4f69ff8b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22360 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example frequency range (logarithmic scale)\n",
    "frequencies = np.logspace(-1, 2, 20)  # Frequencies from 10 to 1000 Hz\n",
    "\n",
    "# Example transfer function H(s) = 1 / (s + 1) evaluated at jw (jω)\n",
    "# H(jω) = 1 / (jω + 1)\n",
    "# s = jw where w = 2 * pi * f\n",
    "s = 1j * 2 * np.pi * frequencies  # jω = j * 2π * frequency\n",
    "\n",
    "# Calculate the transfer function response\n",
    "P = -s*(1/10)/((1/100)*s+1)/((1/10)*s+1)\n",
    "W_2 = 0.21*s/(0.1*s + 1)\n",
    "delta = 1/np.cos(s/1e3)\n",
    "P_tilde = (1+delta*W_2)*P\n",
    "P_env = (1+W_2)*P \n",
    "\n",
    "# Calculate magnitude and phase\n",
    "magnitude = np.abs(P)\n",
    "phase = np.angle(P)\n",
    "\n",
    "magnitude_2 = np.abs(P_tilde)\n",
    "phase_2 = np.angle(P_tilde)\n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, magnitude, 'x')  # Plot magnitude in dB\n",
    "ax1.set_ylabel('Gain')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylim([-0.5,1.5])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, phase/2, 'x')  # Plot phase in degrees\n",
    "ax2.set_ylabel('Phase')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylim([-1.5, 1.5])\n",
    "#ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('nom_dis.png')\n",
    "ax1.set_title('Bode Plot - Discretized')\n",
    "fig.savefig('nom_dis_labeled.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3de5e17b-6a3e-4ef2-b24a-30373135286d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.06270696 0.0900099  0.12891192 0.18379078 0.25967539 0.36060425\n",
      " 0.48559287 0.62301098 0.75012556 0.84460884 0.89698469 0.90811321\n",
      " 0.87900145 0.80767598 0.69664593 0.56179549 0.42761927 0.31261995\n",
      " 0.22310444 0.15715682] [-0.8199144  -0.83498485 -0.85653951 -0.88718468 -0.93023678 -0.98933846\n",
      " -1.06715042 -1.16291424 -1.27062149 -1.38134313 -1.48910635  1.54760516\n",
      "  1.44179513  1.33236964  1.2219375   1.11848039  1.03025951  0.96092209\n",
      "  0.90937199  0.87227083]\n"
     ]
    }
   ],
   "source": [
    "print(magnitude, phase/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c92adaf4-7e95-4732-9b92-6554a11abc40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22704 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate magnitude and phase\n",
    "magnitude = np.abs(P)\n",
    "phase = np.angle(P)\n",
    "\n",
    "beta = 0.5\n",
    "noise = np.random.normal(scale = beta, size = 20)\n",
    "magnitude_2 = np.abs(np.abs(P) + np.abs(P)*np.random.normal(scale = beta, size = 20))\n",
    "phase_2 = np.angle(P)*np.random.normal(scale = beta, size = 20)+np.angle(P)\n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, magnitude_2, 'x')  # Plot magnitude in dB\n",
    "ax1.set_ylabel('Gain')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylim([-0.5,1.5])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, phase_2/2, 'x')  # Plot phase in degrees\n",
    "ax2.set_ylabel('Phase')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylim([-1.5, 1.5])\n",
    "#ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('nom_dis_1.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "548d13c9-d8dd-4bae-bd99-1ec39c7a4d50",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 19908 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "beta = 0.25\n",
    "noise = np.random.normal(scale = beta, size = 20)\n",
    "magnitude_3 = noise \n",
    "phase_3 = noise \n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, magnitude_3, 'x')  # Plot magnitude in dB\n",
    "ax1.set_ylabel('Gain')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylim([-1.5,1.5])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, phase_3/2, 'x')  # Plot phase in degrees\n",
    "ax2.set_ylabel('Phase')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylim([-1.5, 1.5])\n",
    "#ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('nom_dis_2.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cbb8a5ae-d769-4876-9f2c-b00c503c9bc2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36152 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example frequency range (logarithmic scale)\n",
    "frequencies = np.logspace(-1, 2, 10000)  # Frequencies from 10 to 1000 Hz\n",
    "\n",
    "# Example transfer function H(s) = 1 / (s + 1) evaluated at jw (jω)\n",
    "# H(jω) = 1 / (jω + 1)\n",
    "# s = jw where w = 2 * pi * f\n",
    "s = 1j * 2 * np.pi * frequencies  # jω = j * 2π * frequency\n",
    "\n",
    "# Calculate the transfer function response\n",
    "P = -s*(1/10)/((1/100)*s+1)/((1/10)*s+1)\n",
    "W_2 = 0.21*s/(0.1*s + 1)\n",
    "delta = 1/np.cos(s/1e3)\n",
    "P_tilde = (1+delta*W_2)*P\n",
    "P_env = (1+W_2)*P \n",
    "\n",
    "# Calculate magnitude and phase\n",
    "magnitude = np.abs(P)\n",
    "phase = np.angle(P)\n",
    "\n",
    "magnitude_2 = np.abs(P_tilde)\n",
    "phase_2 = np.angle(P_tilde)\n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, 20 * np.log10(magnitude_2), '-r')  # Plot magnitude in dB\n",
    "ax1.set_title('Bode Plot - Continuous')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylabel('Magnitude (dB)')\n",
    "ax1.set_ylim([-60,60])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, np.degrees(phase_2), '-r')  # Plot phase in degrees\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylabel('Phase (Degrees)')\n",
    "ax2.set_ylim([-180,180])\n",
    "ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('per_cont.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fa23e33f-0891-47ee-ac11-a39c43f73f62",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20844 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example frequency range (logarithmic scale)\n",
    "frequencies = np.logspace(-1, 2, 20)  # Frequencies from 10 to 1000 Hz\n",
    "\n",
    "# Example transfer function H(s) = 1 / (s + 1) evaluated at jw (jω)\n",
    "# H(jω) = 1 / (jω + 1)\n",
    "# s = jw where w = 2 * pi * f\n",
    "s = 1j * 2 * np.pi * frequencies  # jω = j * 2π * frequency\n",
    "\n",
    "# Calculate the transfer function response\n",
    "P = -s*(1/10)/((1/100)*s+1)/((1/10)*s+1)\n",
    "W_2 = 0.21*s/(0.1*s + 1)\n",
    "delta = 1/np.cos(s/1e3)\n",
    "P_tilde = (1+delta*W_2)*P\n",
    "P_env = (1+W_2)*P \n",
    "\n",
    "# Calculate magnitude and phase\n",
    "magnitude = np.abs(P_tilde)\n",
    "phase = np.angle(P_tilde)\n",
    "\n",
    "magnitude_2 = np.abs(P_tilde)\n",
    "phase_2 = np.angle(P_tilde)\n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, magnitude, 'xr')  # Plot magnitude in dB\n",
    "ax1.set_ylabel('Gain')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylim([-0.5,3.5])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, phase/2, 'xr')  # Plot phase in degrees\n",
    "ax2.set_ylabel('Phase')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylim([-1.5, 1.5])\n",
    "#ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('per_dis.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c95b8aea-9859-400a-9f23-f164b53e6e6e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21852 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate magnitude and phase\n",
    "magnitude = np.abs(P)\n",
    "phase = np.angle(P)\n",
    "\n",
    "beta = 0.5\n",
    "noise = np.random.normal(scale = beta, size = 20)\n",
    "magnitude_2 = np.abs(np.abs(P_tilde) + np.abs(P_tilde)*np.random.normal(scale = beta, size = 20))\n",
    "phase_2 = np.angle(P_tilde)*np.random.normal(scale = beta, size = 20)+np.angle(P_tilde)\n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, magnitude_2, 'xr')  # Plot magnitude in dB\n",
    "ax1.set_ylabel('Gain')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylim([-0.5,1.5])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, phase_2/2, 'xr')  # Plot phase in degrees\n",
    "ax2.set_ylabel('Phase')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylim([-1.5, 1.5])\n",
    "#ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('per_dis_1.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8d06b583-59aa-4d12-b128-fc4283192e25",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21540 characters]",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "beta = 0.25\n",
    "noise = np.random.normal(scale = beta, size = 20)\n",
    "magnitude_3 = noise\n",
    "phase_3 = noise \n",
    "\n",
    "# Create subplots\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 4))\n",
    "\n",
    "# Magnitude plot\n",
    "ax1.semilogx(frequencies, magnitude_3, 'rx')  # Plot magnitude in dB\n",
    "ax1.set_ylabel('Magnitude')\n",
    "ax1.set_xlabel('Frequency (Hz)')\n",
    "ax1.set_ylim([-1.5,1.5])\n",
    "ax1.grid(True)\n",
    "\n",
    "# Phase plot\n",
    "ax2.semilogx(frequencies, phase_3, 'rx')  # Plot phase in degrees\n",
    "ax2.set_ylabel('Phase')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylim([-1.5, 1.5])\n",
    "#ax2.set_yticks(np.arange(180,-181,-60))\n",
    "ax2.grid(True)\n",
    "\n",
    "# Show the plots\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig('per_dis_2.png')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
