Least Squares FIR

least_squares_fir.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "from scipy.linalg import solve_toeplitz\n",
    "\n",
    "######\n",
    "##\n",
    "## Algorithm that assumes \n",
    "## a fixed length impulse response (length imp_length)\n",
    "## and calculates the least-squares solution of the \n",
    "## convolution equation y = convolve(x, w). With optional l2-regularization.\n",
    "##\n",
    "##\n",
    "## Arguments:\n",
    "## reg - l2 regulariztion strength (optional).\n",
    "## x - input signal length N\n",
    "## y - output signal length N\n",
    "## imp_length - desired impulse response length.\n",
    "##\n",
    "## Returns:\n",
    "## w - Learned impulse response (length imp_length)\n",
    "##\n",
    "## Remember to zero-pad x/y if normal convolution (not circular) is desired.\n",
    "## This is typically the case.  \n",
    "    \n",
    "def least_squares_fir(x,y,imp_length,reg=0.0):\n",
    "\n",
    "    n=len(x)\n",
    "    X = np.fft.rfft(x)\n",
    "    Y = np.fft.rfft(y)\n",
    "\n",
    "    b = np.fft.irfft(np.conj(X)*Y,n)[0:imp_length] # Normal equation vec b\n",
    "    \n",
    "    tmp_a = np.fft.irfft(np.conj(X)*X,n) # Normal Equation matrix A is first (imp_length x imp_length) elements of the circulant matrix made from this vector.\n",
    "    \n",
    "    \n",
    "    \n",
    "    c = tmp_a[0:imp_length]\n",
    "    r = tmp_a[::-1][0:imp_length-1]\n",
    "    r = np.concatenate((tmp_a[[0]],r))\n",
    "\n",
    "    r[0]+=reg \n",
    "    c[0]+=reg\n",
    "    \n",
    "    return solve_toeplitz( (c,r), b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have:\n",
    "\n",
    "\n",
    "Wet Signal: $\\vec{y}$ (length n)\n",
    "\n",
    "Dry Signal: $\\vec{x}$ (length n)\n",
    "\n",
    "Impulse response $\\vec{w}$ (length m < n)\n",
    "\n",
    "Convolution can be described as a Linear problem given by:\n",
    "$\\vec{y} = C I_{n,m} \\vec{w}$.\n",
    "\n",
    "$I_{n,m}$ is a linear operator that zero-pads $\\vec{w}$ to length n.\n",
    "$C$ is the circulant matrix of $\\vec{x}$.\n",
    "\n",
    "$C = W^H X W$, where X is the complex diagonal matrix with the fourier-components of x, and W is the DFT(n) matrix.\n",
    "\n",
    "Normal equation for the system can be written as:\n",
    "\n",
    "$A^H\\vec{y} = I_{n,m}^T W^H X^H W\\vec{y}$.\n",
    "\n",
    "$A^HA = I_{n,m}^T (W^HX^HXW) I_{n,m}$\n",
    "\n",
    "Algorithmically:\n",
    "\n",
    "$A^H\\vec{y}$ is first M elements of IFFT(conj(FFT($\\vec{x}$)).*FFT($\\vec{y}$))\n",
    "\n",
    "$A^HA$: is the first (m x m) elements of the circulant matrix made from IFFT(conj(FFT($\\vec{x}$)).*FFT($\\vec{x}$)))\n",
    "\n",
    "The resulting matrix is itself not circulant, but topelitz, and the system can be solved in $O(m^2)$ using Levinson recursion method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}