{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "348dd2c6",
   "metadata": {},
   "source": [
    "# Noise suppression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b114e7b",
   "metadata": {},
   "source": [
    "Very simplified noise suppression algorithm to show how to use the CMSIS-DSP Python wrapper and how it can help translating a float algorithm to a fixed point one ready to be implemented on a Cortex-M or Cortex-A processor."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba37a7e0",
   "metadata": {},
   "source": [
    "The CMSIS-DSP package must first be installed. This command will take some time to execute since the CMSIS-DSP will be built."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19680a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install cmsisdsp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "934cf8b0",
   "metadata": {},
   "source": [
    "The packages required for this example must be loaded:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bfb671b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerical packages\n",
    "import cmsisdsp as dsp \n",
    "import cmsisdsp.fixedpoint as fix\n",
    "import numpy as np\n",
    "from numpy.lib.stride_tricks import sliding_window_view\n",
    "from scipy.signal.windows import hann\n",
    "\n",
    "# Package for plotting\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Package to display Audio widget in the notebook and to upload sound files\n",
    "import ipywidgets\n",
    "from ipywidgets import FileUpload\n",
    "from IPython.display import display,Audio\n",
    "\n",
    "# To convert sound file to NumPy array\n",
    "import io\n",
    "import soundfile as sf\n",
    "\n",
    "# To load test pattern from Arm VHT Echo Canceller demo\n",
    "from urllib.request import urlopen\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae39ac8d",
   "metadata": {},
   "source": [
    "## Loading the audio file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "153b71bb",
   "metadata": {},
   "source": [
    "You can either upload a .wav file or download the test pattern used in our echo canceller demo running on Arm Virtual Hardware."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8130c184",
   "metadata": {},
   "source": [
    "### Test pattern from Arm Virtual Hardware demo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70e33e66",
   "metadata": {},
   "source": [
    "This is based on the Yes/No pattern from the TensorFlow Lite Microspeech demo but with added noise and a sequence of Yes/No which is different. It is working with the Q15 implementation (which is not the case of all speech pattern you may test)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3c298f96",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_pattern_url=\"https://github.com/ARM-software/VHT-SystemModeling/blob/main/EchoCanceller/sounds/yesno.wav?raw=true\"\n",
    "f = urlopen(test_pattern_url)\n",
    "filedata = f.read()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fa31b19",
   "metadata": {},
   "source": [
    "### Wav file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "918a12b0",
   "metadata": {},
   "source": [
    "Select a .wav file to upload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 451,
   "id": "0f4d0cff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d582a97e382f4858a5e07321af8b1243",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FileUpload(value={}, accept='wav', description='Upload')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "uploader=FileUpload(accept='wav',multiple=False)\n",
    "display(uploader)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0853e7a9",
   "metadata": {},
   "source": [
    "The two next cell is loading the audio data from the uploaded file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 436,
   "id": "ad14cf43",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename=list(uploader.value.keys())[0]\n",
    "filedata=uploader.value[filename][\"content\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95948301",
   "metadata": {},
   "source": [
    "### Audio player"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66ab2923",
   "metadata": {},
   "source": [
    "You can now listen to the sound to check it is available in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "46f7169b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480952 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audio=Audio(data=filedata,autoplay=False)\n",
    "audio"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2eaeae6",
   "metadata": {},
   "source": [
    "A NumPy array is created from the audio and a Q15 and Q31 versions are created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9b7c4b14",
   "metadata": {},
   "outputs": [],
   "source": [
    "data, samplerate = sf.read(io.BytesIO(filedata))\n",
    "if len(data.shape)>1:\n",
    "    data=data[:,0]\n",
    "data = data + np.random.randn(len(data))*0.00\n",
    "data=data/np.max(np.abs(data))\n",
    "dataQ15 = fix.toQ15(data)\n",
    "dataQ31 = fix.toQ31(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6dfbffa",
   "metadata": {},
   "source": [
    "The audio waveform is plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1f446dfc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13868 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "677cbf18",
   "metadata": {},
   "source": [
    "The audio is split into overlapping windows. Those windows will be used in the notebook to implement the processing.\n",
    "The signal will be rebuilt at the end with an overlap and add.\n",
    "\n",
    "If the window durations are changed, the thresholds for the VAD will also have to be changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4427ae73",
   "metadata": {},
   "outputs": [],
   "source": [
    "winDuration=30e-3\n",
    "winOverlap=15e-3\n",
    "\n",
    "# Decrease the duration so that Q15 is more likely to work with sampling freq of 48 kHz\n",
    "# (to avoid too long FFT length which may saturate in Q15)\n",
    "winDuration = winDuration/6 \n",
    "winOverlap = winOverlap/6\n",
    "\n",
    "winLength=int(np.floor(samplerate*winDuration))\n",
    "winOverlap=int(np.floor(samplerate*winOverlap))\n",
    "slices=sliding_window_view(data,winLength)[::winOverlap,:]\n",
    "slices_q15=sliding_window_view(dataQ15,winLength)[::winOverlap,:]\n",
    "slices_q31=sliding_window_view(dataQ31,winLength)[::winOverlap,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f49559a",
   "metadata": {},
   "source": [
    "## Energy VAD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ffe045f",
   "metadata": {},
   "source": [
    "This is a very simplified VAD algorithm based on energy threshold. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3159c218",
   "metadata": {},
   "source": [
    "### Double implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cfa5a997",
   "metadata": {},
   "outputs": [],
   "source": [
    "def clean_vad(v):\n",
    "    v = np.hstack([[0],v,[0]])\n",
    "    # Remove isolated peak\n",
    "    vmin=[np.min(l) for l in sliding_window_view(v,3)]\n",
    "    vmin = np.hstack([[0,0],vmin,[0]])\n",
    "    # Remove isolated hole\n",
    "    vmax=[np.max(l) for l in sliding_window_view(vmin,4)]\n",
    "    return(vmax)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "23a57a6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Energy of the window\n",
    "def signal_energy(window):\n",
    "    w = window - np.mean(window)\n",
    "    return(10*np.log10(np.sum(window * window)))\n",
    "    \n",
    "# Comparison with a threshold\n",
    "def signal_vad(window):\n",
    "    if signal_energy(window)>-11:\n",
    "        return(1)\n",
    "    else:\n",
    "        return(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83c9e318",
   "metadata": {},
   "source": [
    "Plotting of the signal and the VAD detection. It can be used to tune the threshold for the audio signal used in the notebook. It is unlikely that only one threshold will work with any audio signal. It is a very simplified VAD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d8cf68d2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1dbe9025c00>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16812 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_,ax=plt.subplots(1,1)\n",
    "cleaned=clean_vad([signal_vad(w) for w in slices])\n",
    "vad = np.array([[w]*winOverlap for w in cleaned]).flatten()\n",
    "#signal_energy = np.array([[signal_energy(w)]*winOverlap for w in slices]).flatten()\n",
    "ax.plot(data)\n",
    "ax.plot(vad)\n",
    "#ax.plot(signal_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3e2e0d4",
   "metadata": {},
   "source": [
    "### CMSIS-DSP Q15 implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0885101f",
   "metadata": {},
   "source": [
    "The same functions but in q15 and using CMSIS-DSP API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "37e41c05",
   "metadata": {},
   "outputs": [],
   "source": [
    "def signal_energy_q15(window):\n",
    "    mean=dsp.arm_mean_q15(window)\n",
    "    window=dsp.arm_offset_q15(window,-mean)\n",
    "    energy=dsp.arm_power_q15(window)\n",
    "    # Energy is not in q15 (refer to CMSIS-DSP documentation)\n",
    "    energy=dsp.ssat(energy>>20,16)\n",
    "    # The shift is by more than 15, so the result is not in q15\n",
    "    # (With a shift smaller there is a saturation)\n",
    "    dB=dsp.arm_vlog_q15([energy])\n",
    "    # The output of the vlog is not in q15\n",
    "    # The multiplication by 10 is missing.\n",
    "    # The result of this function is not \"equal\" to the float implementation due to the different\n",
    "    # formats used in the intermediate computations.\n",
    "    return(dB[0])\n",
    "    \n",
    "def signal_vad_q15(window):\n",
    "    # The threshold is not easily related to the float implementation because\n",
    "    # of the different intermediate formats used in the fixed poitn implementation\n",
    "    if signal_energy_q15(window)>fix.toQ15(-0.38):\n",
    "        return(1)\n",
    "    else:\n",
    "        return(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89b856d2",
   "metadata": {},
   "source": [
    "Plot of the signal and VAD which can be used to tune the threshdold. The threshold is different because we do not apply the multiplication by 10 and because the output format of the vlog is not Q15. And the output format of the power function is also not q15 (the formats are described in CMSIS-DSP documentation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "276749ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1dbeb48afb0>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16684 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_,ax=plt.subplots(1,1)\n",
    "cleaned=clean_vad([signal_vad_q15(w) for w in slices_q15])\n",
    "vad_q15 = np.array([[w]*winOverlap for w in cleaned]).flatten()\n",
    "#vad = np.array([[signal_vad(w)]*winOverlap for w in slices]).flatten()\n",
    "ax.plot(data)\n",
    "#ax.plot(vad)\n",
    "ax.plot(vad_q15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1576979d",
   "metadata": {},
   "source": [
    "### CMSIS-DSP Q31 Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "51e5db4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def signal_energy_q31(window):\n",
    "    mean=dsp.arm_mean_q31(window)\n",
    "    window=dsp.arm_offset_q31(window,-mean)\n",
    "    energy=dsp.arm_power_q31(window)\n",
    "    # Energy is not in q31 (refer to CMSIS-DSP documentation)\n",
    "    # If the scaling is not enough the int conversion will fail.\n",
    "    scaled=energy>>25 \n",
    "    if scaled > 0x7FFFFFFF:\n",
    "        energy=0x7FFFFFFF\n",
    "    elif scaled < -0x80000000:\n",
    "        energy=-0x80000000\n",
    "    else:\n",
    "        energy=int(energy>>25)\n",
    "    # The shift is by more than 17, so the result is not in q31\n",
    "    # (With a shift smaller there is a saturation and int conversion is failing)\n",
    "    dB=dsp.arm_vlog_q31([energy])\n",
    "    # The output of the vlog is not in q15\n",
    "    # The multiplication by 10 is missing.\n",
    "    # The result of this function is not \"equal\" to the float implementation due to the different\n",
    "    # formats used in the intermediate computations.\n",
    "    return(dB[0])\n",
    "    \n",
    "def signal_vad_q31(window):\n",
    "    # The threshold is not easily related to the float implementation because\n",
    "    # of the different intermediate formats used in the fixed poitn implementation\n",
    "    if signal_energy_q31(window)>fix.toQ31(-0.255):\n",
    "        return(1)\n",
    "    else:\n",
    "        return(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3e9e5692",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1dbeb50bee0>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16724 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_,ax=plt.subplots(1,1)\n",
    "cleaned=clean_vad([signal_vad_q31(w) for w in slices_q31])\n",
    "vad_q31 = np.array([[w]*winOverlap for w in cleaned]).flatten()\n",
    "ax.plot(data)\n",
    "#ax.plot(vad)\n",
    "ax.plot(vad_q31)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e278fc7",
   "metadata": {},
   "source": [
    "## Noise suppression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "964c42e4",
   "metadata": {},
   "source": [
    "Each window of samples extracted from the signal is multiplied by a Hann window in below algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d780cca7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14476 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "window=hann(winLength,sym=False)\n",
    "plt.plot(window)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35aa4c9c",
   "metadata": {},
   "source": [
    "### Double reference implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f82508ce",
   "metadata": {},
   "source": [
    "Test of the overlap and add:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4ae681e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13920 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Overlap and add\n",
    "offsets = range(0, len(data),winOverlap)\n",
    "offsets=offsets[0:len(slices)]\n",
    "res=np.zeros(len(data))\n",
    "i=0\n",
    "for n in offsets:\n",
    "    res[n:n+winLength] += slices[i]*window\n",
    "    i=i+1\n",
    "plt.plot(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e8c4ce02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480952 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audio2=Audio(data=res,rate=samplerate,autoplay=False)\n",
    "audio2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5beccc",
   "metadata": {},
   "source": [
    "#### Reference noise suppression algorithm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "cca944e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fft_length(length):\n",
    "    result=2\n",
    "    fft_shift=1\n",
    "    while result < length:\n",
    "        result = 2*result\n",
    "        fft_shift = fft_shift + 1\n",
    "    return(result,fft_shift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5706f45c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NoiseSuppression():\n",
    "    def __init__(self,slices):\n",
    "        self._windowLength=len(slices[0])\n",
    "        self._fftLen,self._fftShift=fft_length(self._windowLength)\n",
    "        \n",
    "        self._padding_left=(self._fftLen - self._windowLength)//2 \n",
    "        self._padding_right=self._fftLen- self._windowLength-self._padding_left\n",
    "         \n",
    "        self._signal=[]\n",
    "        self._slices=slices\n",
    "        self._window=None\n",
    "        \n",
    "    def window_and_pad(self,w):\n",
    "        if w.dtype==np.int32:\n",
    "            w=dsp.arm_mult_q31(w,self._window)\n",
    "        elif w.dtype==np.int16:\n",
    "            w=dsp.arm_mult_q15(w,self._window)\n",
    "        else:\n",
    "            w = w*self._window\n",
    "        sig=np.hstack([np.zeros(self._padding_left,dtype=w.dtype),w,np.zeros(self._padding_right,dtype=w.dtype)])\n",
    "        return(sig)\n",
    "    \n",
    "    def remove_padding(self,w):\n",
    "        return(w[self._padding_left:self._padding_left+self._windowLength])\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "69d451f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NoiseSuppressionReference(NoiseSuppression):\n",
    "    def __init__(self,slices):\n",
    "        # In a better version this could be computed from the signal length by taking the\n",
    "        # smaller power of two greater than the signal length.\n",
    "        NoiseSuppression.__init__(self,slices)\n",
    "        \n",
    "        # Compute the vad signal\n",
    "        self._vad=clean_vad([signal_vad(w) for w in slices])\n",
    "        self._noise=np.zeros(self._fftLen)\n",
    "        # The Hann window\n",
    "        self._window=hann(self._windowLength,sym=False)\n",
    "        \n",
    "    # Subtract the noise\n",
    "    def subnoise(self,v):\n",
    "        # This is a Wiener estimate\n",
    "        energy = v * np.conj(v) + 1e-6\n",
    "        \n",
    "        scaling = (energy - self._noise)/energy\n",
    "        scaling[scaling<0] = 0\n",
    "        \n",
    "        return(v * scaling)\n",
    "    \n",
    "    def remove_noise(self,w):\n",
    "        # We pad the signal with zero. It assumes that the padding can be divided by 2.\n",
    "        # In a better implementation we would manage also the odd case.\n",
    "        # The padding is required because the FFT has a length which is greater than the length of\n",
    "        # the window\n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        # FFT\n",
    "        fft=np.fft.fft(sig)\n",
    "        # Noise suppression\n",
    "        fft = self.subnoise(fft)\n",
    "        # IFFT\n",
    "        res=np.fft.ifft(fft)\n",
    "        # We assume the result should be real so we just ignore the imaginary part\n",
    "        res=np.real(res)\n",
    "        # We remove the padding\n",
    "        res=self.remove_padding(res)\n",
    "        return(res)\n",
    "    \n",
    "   \n",
    "    \n",
    "    def estimate_noise(self,w):\n",
    "        # Compute the padded signal\n",
    "        sig=self.window_and_pad(w)\n",
    "        fft=np.fft.fft(sig)\n",
    "        \n",
    "        # Estimate the noise energy\n",
    "        self._noise = np.abs(fft)*np.abs(fft)\n",
    "        \n",
    "        # Remove the noise\n",
    "        fft = self.subnoise(fft)\n",
    "        \n",
    "        # IFFT and we assume the result is real so we ignore imaginary part\n",
    "        res=np.fft.ifft(fft)\n",
    "        res=np.real(res)\n",
    "        res=self.remove_padding(res)\n",
    "        return(res)\n",
    "        \n",
    "    # Process all the windows using the VAD detection\n",
    "    def nr(self):\n",
    "        for (w,v) in zip(self._slices,self._vad):\n",
    "            result=None\n",
    "            if v==1:\n",
    "                # If voice detected, we only remove the noise\n",
    "                result=self.remove_noise(w)\n",
    "            else:\n",
    "                # If no voice detected, we update the noise estimate\n",
    "                result=self.estimate_noise(w)\n",
    "            self._signal.append(result)\n",
    "        \n",
    "    # Overlap and add to rebuild the signal\n",
    "    def overlap_and_add(self):\n",
    "        offsets = range(0, len(self._signal)*winOverlap,winOverlap)\n",
    "        offsets=offsets[0:len(self._signal)]\n",
    "        res=np.zeros(len(data))\n",
    "        i=0\n",
    "        for n in offsets:\n",
    "            res[n:n+winLength]+=self._signal[i]\n",
    "            i=i+1\n",
    "        return(res)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3915dc37",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12452 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=NoiseSuppressionReference(slices)\n",
    "n.nr()\n",
    "cleaned=n.overlap_and_add()\n",
    "plt.plot(cleaned)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f12a6f54",
   "metadata": {},
   "source": [
    "The audio to check the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3091d6f5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480952 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audioRef=Audio(data=cleaned,rate=samplerate,autoplay=False)\n",
    "audioRef"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ee5a76d",
   "metadata": {},
   "source": [
    "### CMSIS-DSP Implementations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73276ee1",
   "metadata": {},
   "source": [
    "Test of the overlap and add. The Hann window is converted to Q15. The CMSIS-DSP functions are not used for this conversions because it is assumed that this conversion is done outside of CMSIS-DSP and the final code (Python or C) is using an array of Q15 values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d0313891",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13760 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "offsets = range(0, len(data),winOverlap)\n",
    "offsets=offsets[0:len(slices_q15)]\n",
    "res=np.zeros(len(data))\n",
    "window_q15=fix.toQ15(window)\n",
    "i=0\n",
    "for n in offsets:\n",
    "    w = dsp.arm_mult_q15(slices_q15[i],window_q15)\n",
    "    res[n:n+winLength] = dsp.arm_add_q15(res[n:n+winLength],w)\n",
    "    i=i+1\n",
    "res_q15=fix.Q15toF32(res)\n",
    "plt.plot(res_q15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b7518435",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480952 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audio4=Audio(data=res_q15,rate=samplerate,autoplay=False)\n",
    "audio4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dafab1bf",
   "metadata": {},
   "source": [
    "Utilities functions since there is no specific datatype for complex in CMSIS-DSP. An array of complex is thus represented as an array of reals (of twice the number of samples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "3a1be0da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def imToReal1D(a):\n",
    "    ar=np.zeros(np.array(a.shape) * 2)\n",
    "    ar[0::2]=a.real\n",
    "    ar[1::2]=a.imag\n",
    "    return(ar)\n",
    "\n",
    "def realToIm1D(ar):\n",
    "    return(ar[0::2] + 1j * ar[1::2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc594498",
   "metadata": {},
   "source": [
    "#### CMSIS-DSP F32 noise suppression algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "328c0eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NoiseSuppressionF32(NoiseSuppression):\n",
    "    def __init__(self,slices):\n",
    "        NoiseSuppression.__init__(self,slices)\n",
    "\n",
    "        # VAD signal\n",
    "        self._vad= clean_vad(np.array([signal_vad(w) for w in slices]))\n",
    "        self._noise=np.zeros(self._fftLen,dtype=np.float32)\n",
    "        # The Hann window\n",
    "        self._window=hann(self._windowLength,sym=False)\n",
    "        # CFTT F32 instance\n",
    "        self._cfftF32=dsp.arm_cfft_instance_f32()\n",
    "        status=dsp.arm_cfft_init_f32(self._cfftF32,self._fftLen)\n",
    "        \n",
    "    \n",
    "    # Subtract the noise\n",
    "    def subnoise(self,v):\n",
    "        energy = dsp.arm_cmplx_mag_squared_f32(v)\n",
    "        # To avoid division by zeros\n",
    "        energy = dsp.arm_offset_f32(energy,1e-6)\n",
    "        \n",
    "        temp = dsp.arm_sub_f32(energy , self._noise)\n",
    "        # C implementation will have to loop on the array and replace negative vaues\n",
    "        temp[temp<0]=0\n",
    "            \n",
    "        scaling = np.zeros(len(temp),dtype=np.float32)\n",
    "        # C implementation will have to loop on the arrays here\n",
    "        scaling = temp / energy\n",
    "        \n",
    "        res=dsp.arm_cmplx_mult_real_f32(v,scaling)    \n",
    "                    \n",
    "        return(res)\n",
    "    \n",
    "    def remove_noise(self,w):\n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        # Convert real signal to complex        \n",
    "        signalR=np.zeros(len(sig) * 2)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "        resultR = dsp.arm_cfft_f32(self._cfftF32,signalR,0,1)\n",
    "        \n",
    "        # Remove the noise\n",
    "        resultR = self.subnoise(resultR)\n",
    "\n",
    "        # Inverse FFT and we keep only the real part\n",
    "        res = dsp.arm_cfft_f32(self._cfftF32,resultR,1,1)*self._fftLen\n",
    "        res=res[0::2]\n",
    "        # Remove the padding\n",
    "        res=self.remove_padding(res)\n",
    "        return(res)\n",
    "    \n",
    "    def estimate_noise(self,w):\n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        signalR=np.zeros(len(sig) * 2)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "        resultR = dsp.arm_cfft_f32(self._cfftF32,signalR,0,1)\n",
    "\n",
    "        self._noise = dsp.arm_cmplx_mag_squared_f32(resultR)\n",
    "        \n",
    "        # Set noise to zero. In reference code we substract the noise.\n",
    "        # It is giving the same result with the noise estimation used.\n",
    "        resultR = np.zeros(len(resultR),dtype=np.float32)\n",
    "        \n",
    "        # Not really needed since we set the values to zero. We could\n",
    "        # just set the output to zero.\n",
    "        # But in a more compelx algorithm we may still need to do the IFFT\n",
    "        # so the line is kept.\n",
    "        res = dsp.arm_cfft_f32(self._cfftF32,resultR,1,1)*self._fftLen\n",
    "        \n",
    "        res=res[0::2]\n",
    "        res=self.remove_padding(res)\n",
    "        return(res)\n",
    "        \n",
    "    def nr(self):\n",
    "        for (w,v) in zip(self._slices,self._vad):\n",
    "            result=None\n",
    "            if v==1:\n",
    "                result=self.remove_noise(w)\n",
    "            else:\n",
    "                result=self.estimate_noise(w)\n",
    "            self._signal.append(result)\n",
    "        \n",
    "    def overlap_and_add(self):\n",
    "        nbSamples = len(self._signal)*winOverlap\n",
    "        offsets = range(0, nbSamples,winOverlap)\n",
    "        offsets=offsets[0:len(self._signal)]\n",
    "        res=np.zeros(nbSamples)\n",
    "        i=0\n",
    "        for n in offsets:\n",
    "            res[n:n+winLength] = dsp.arm_add_f32(res[n:n+winLength],self._signal[i])\n",
    "            i=i+1\n",
    "        return(res)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "cca48f75",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 10436 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=NoiseSuppressionF32(slices)\n",
    "n.nr()\n",
    "cleaned_f32=n.overlap_and_add()\n",
    "plt.plot(cleaned_f32)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "fcfd59e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480808 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audioF32=Audio(data=cleaned_f32,rate=samplerate,autoplay=False)\n",
    "audioF32"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "283082da",
   "metadata": {},
   "source": [
    "#### CMSIS-DSP Q31 noise suppression algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "8ccaeae1",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NoiseSuppressionQ31(NoiseSuppression):\n",
    "    def __init__(self,slices):\n",
    "        NoiseSuppression.__init__(self,slices)\n",
    "        \n",
    "        # VAD signal\n",
    "        self._vad= clean_vad(np.array([signal_vad_q31(w) for w in slices]))\n",
    "        self._noise=np.zeros(self._fftLen,dtype=np.int32)\n",
    "        # Q31 version of the Hann window\n",
    "        self._window=fix.toQ31(hann(self._windowLength,sym=False))\n",
    "        # CFTT Q31 instance\n",
    "        self._cfftQ31=dsp.arm_cfft_instance_q31()\n",
    "        status=dsp.arm_cfft_init_q31(self._cfftQ31,self._fftLen)\n",
    "        \n",
    "    \n",
    "    # Subtract the noise\n",
    "    def subnoise(self,v):\n",
    "        \n",
    "        energy = dsp.arm_cmplx_mag_squared_q31(v)\n",
    "        \n",
    "        temp = dsp.arm_sub_q31(energy , self._noise)\n",
    "        temp[temp<0]=0\n",
    "        \n",
    "        # Used for computing the division (energy - noise) / energy\n",
    "        scalingQ31 = np.zeros(len(temp),dtype=np.int32)\n",
    "        shift = np.zeros(len(temp),dtype=np.int32)\n",
    "        \n",
    "        # The scaling factor (energy - noise) / energy is computed\n",
    "        k=0\n",
    "        # We assume that |energy-noise|<=energy\n",
    "        # otherwise we set scaling to 1\n",
    "        # If energy is 0 we also set scaling to 1\n",
    "        # When a == b shiftVal is equal to 1 because 1 is represented \n",
    "        # as 0x40000000 and shift of 1 instead of 0x7FFF for output of division\n",
    "        # We manage this case separately\n",
    "        for a,b in zip(temp,energy):\n",
    "            quotient=0x7FFFFFFF\n",
    "            shiftVal=0\n",
    "            if b!=0 and a!=b:\n",
    "                # We compute the quotient\n",
    "                status,quotient,shiftVal = dsp.arm_divide_q31(a,b)\n",
    "                if shiftVal > 0:\n",
    "                    quotient=0x7FFFFFFF\n",
    "                    shiftVal = 0\n",
    "        \n",
    "            scalingQ31[k] = quotient\n",
    "            shift[k] = shiftVal\n",
    "            \n",
    "            k = k + 1\n",
    "        \n",
    "       \n",
    "        # We scale the values\n",
    "        res=dsp.arm_cmplx_mult_real_q31(v,scalingQ31)    \n",
    "            \n",
    "        return(res)\n",
    "    \n",
    "    def remove_noise(self,w):\n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        # Convert to complex        \n",
    "        signalR=np.zeros(len(sig) * 2,dtype=np.int32)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "\n",
    "        resultR = dsp.arm_cfft_q31(self._cfftQ31,signalR,0,1)\n",
    "        \n",
    "        resultR = self.subnoise(resultR)\n",
    "\n",
    "        res = dsp.arm_cfft_q31(self._cfftQ31,resultR,1,1)\n",
    "        res = dsp.arm_shift_q31(res,self._fftShift)\n",
    "        \n",
    "        res=res[0::2]\n",
    "        res=self.remove_padding(res)\n",
    "        return(res)\n",
    "    \n",
    "    def estimate_noise(self,w):\n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        signalR=np.zeros(len(sig) * 2,dtype=np.int32)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "        resultR = dsp.arm_cfft_q31(self._cfftQ31,signalR,0,1)\n",
    "\n",
    "        self._noise = dsp.arm_cmplx_mag_squared_q31(resultR)\n",
    "        \n",
    "        resultR = np.zeros(len(resultR),dtype=np.int32)\n",
    "        \n",
    "        res = dsp.arm_cfft_q31(self._cfftQ31,resultR,1,1)\n",
    "        res = dsp.arm_shift_q31(res,self._fftShift)\n",
    "        \n",
    "        res=res[0::2]\n",
    "        res=self.remove_padding(res)\n",
    "        return(res)\n",
    "        \n",
    "    def nr(self):\n",
    "        for (w,v) in zip(self._slices,self._vad):\n",
    "            result=None\n",
    "            if v==1:\n",
    "                result=self.remove_noise(w)\n",
    "            else:\n",
    "                result=self.estimate_noise(w)\n",
    "            self._signal.append(result)\n",
    "        \n",
    "    def overlap_and_add(self):\n",
    "        nbSamples = len(self._signal)*winOverlap\n",
    "        offsets = range(0, nbSamples,winOverlap)\n",
    "        offsets=offsets[0:len(self._signal)]\n",
    "        res=np.zeros(nbSamples)\n",
    "        i=0\n",
    "        for n in offsets:\n",
    "            res[n:n+winLength] = dsp.arm_add_q31(res[n:n+winLength],self._signal[i])\n",
    "            i=i+1\n",
    "        return(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "0ccb9512",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12448 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=NoiseSuppressionQ31(slices_q31)\n",
    "n.nr()\n",
    "cleaned_q31=n.overlap_and_add()\n",
    "plt.plot(fix.Q31toF32(cleaned_q31))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "452cbc86",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480808 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audioQ31=Audio(data=fix.Q31toF32(cleaned_q31),rate=samplerate,autoplay=False)\n",
    "audioQ31"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "510a343c",
   "metadata": {},
   "source": [
    "#### CMSIS-DSP Q15 noise suppression algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ad0778f",
   "metadata": {},
   "source": [
    "When the sampling rate is higher, the FFT length is bigger and on some signals it gives some internal saturations in the FFT Q15.\n",
    "\n",
    "In those cases, the algorithm will not work.\n",
    "\n",
    "You can decrease the sampling rate or the window duration. \n",
    "\n",
    "Scaling the signal to remove the saturations will not work if the background noise is too small because there won't be enough bits left to estimate the noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "a23ff9eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NoiseSuppressionQ15(NoiseSuppression):\n",
    "    def __init__(self,slices):\n",
    "        NoiseSuppression.__init__(self,slices)\n",
    "        \n",
    "        # VAD signal\n",
    "        self._vad= clean_vad(np.array([signal_vad_q15(w) for w in slices]))\n",
    "        self._noise=np.zeros(self._fftLen,dtype=np.int32)\n",
    "        # Q15 version of the Hann window\n",
    "        self._window=fix.toQ15(hann(self._windowLength,sym=False))\n",
    "        # CFTT Q15 instance\n",
    "        self._cfftQ15=dsp.arm_cfft_instance_q15()\n",
    "        status=dsp.arm_cfft_init_q15(self._cfftQ15,self._fftLen)\n",
    "        \n",
    "        self._noise_status = -1 \n",
    "        self._noise_max = 0x7FFF\n",
    "        \n",
    "    \n",
    "    # Subtract the noise\n",
    "    def subnoise(self,v,status,the_max):\n",
    "        \n",
    "        # We cannot compute the energy in Q15 because otherwise lots of values are 0.\n",
    "        # The noise signal is too small for its energy to be representable in Q15.\n",
    "        # So, we convert into Q31 and do noise subtraction in Q31\n",
    "        vq31 = dsp.arm_q15_to_q31(v)\n",
    "        energy = dsp.arm_cmplx_mag_squared_q31(vq31)\n",
    "        \n",
    "        # The energy for the signal and noise have been computed on rescaled signal.\n",
    "        # So we remove the scaling from the values before computing the ratio (energy-noise)/energy\n",
    "        if status==0:\n",
    "            the_max_q31=dsp.arm_q15_to_q31([the_max])[0]\n",
    "            energy=dsp.arm_scale_q31(energy,the_max_q31,0)\n",
    "            energy=dsp.arm_scale_q31(energy,the_max_q31,0)\n",
    "            \n",
    "        noise = self._noise \n",
    "        if self._noise_status==0:\n",
    "            the_max_q31=dsp.arm_q15_to_q31([self._noise_max])[0]\n",
    "            noise=dsp.arm_scale_q31(noise,the_max_q31,0)\n",
    "            noise=dsp.arm_scale_q31(noise,the_max_q31,0)\n",
    "            \n",
    "        \n",
    "        temp = dsp.arm_sub_q31(energy , noise)\n",
    "        temp[temp<0]=0\n",
    "        \n",
    "        scalingQ31 = np.zeros(len(temp),dtype=np.int32)\n",
    "        shift = np.zeros(len(temp),dtype=np.int32)\n",
    "        \n",
    "        # The scaling factor (energy - noise) / energy is computed\n",
    "        k=0\n",
    "        # We assume that |energy-noise|<=energy\n",
    "        # otherwise we set scaling to 1\n",
    "        # If energy is 0 we also set scaling to 1\n",
    "        # When a == b shiftVal is equal to 1 because 1 is represented \n",
    "        # as 0x40000000 and shift of 1 instead of 0x7FFF for output of division\n",
    "        # We manage this case separately\n",
    "        for a,b in zip(temp,energy):\n",
    "            quotient=0x7FFFFFFF\n",
    "            shiftVal=0\n",
    "            if b!=0 and a!=b:\n",
    "                # We compute the quotient\n",
    "                status,quotient,shiftVal = dsp.arm_divide_q31(a,b)\n",
    "                if shiftVal > 0:\n",
    "                    quotient=0x7FFFFFFF\n",
    "                    shiftVal = 0\n",
    "        \n",
    "            scalingQ31[k] = quotient\n",
    "            shift[k] = shiftVal\n",
    "            \n",
    "            k = k + 1\n",
    "        \n",
    "       \n",
    "        res=dsp.arm_cmplx_mult_real_q31(vq31,scalingQ31)     \n",
    "        resQ15 = dsp.arm_q31_to_q15(res)\n",
    "            \n",
    "        return(resQ15)\n",
    "    \n",
    "    # To have maximum accuracy with the FFT Q15, the signal is rescaled before computing the FFT\n",
    "    def rescale(self,w):\n",
    "        the_max,index=dsp.arm_absmax_q15(w)\n",
    "        \n",
    "        quotient=0x7FFF \n",
    "        the_shift=0\n",
    "        status = -1\n",
    "        if the_max != 0:\n",
    "            status,quotient,the_shift = dsp.arm_divide_q15(0x7FFF,the_max)\n",
    "            if status == 0:\n",
    "                w=dsp.arm_scale_q15(w,quotient,the_shift)\n",
    "        return(w,status,the_max)\n",
    "    \n",
    "    # The scaling is removed after computing the IFFT\n",
    "    def undo_scale(self,w,the_max):\n",
    "        w=dsp.arm_scale_q15(w,the_max,0)\n",
    "        return(w)\n",
    "    \n",
    "    \n",
    "    def remove_noise(self,w):\n",
    "        w,status,the_max = self.rescale(w)\n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        # Convert to complex        \n",
    "        signalR=np.zeros(len(sig) * 2,dtype=np.int16)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "\n",
    "        resultR = dsp.arm_cfft_q15(self._cfftQ15,signalR,0,1)\n",
    "        \n",
    "        resultR = self.subnoise(resultR,status,the_max)\n",
    "\n",
    "        res = dsp.arm_cfft_q15(self._cfftQ15,resultR,1,1)\n",
    "        res = dsp.arm_shift_q15(res,self._fftShift)\n",
    "        \n",
    "        res=res[0::2]\n",
    "        res=self.remove_padding(res)\n",
    "        \n",
    "        if status == 0:\n",
    "            res=self.undo_scale(res,the_max)\n",
    "        return(res)\n",
    "    \n",
    "    def estimate_noise(self,w):\n",
    "        w,status,the_max = self.rescale(w)\n",
    "        self._noise_status = status \n",
    "        self._noise_max = the_max\n",
    "        \n",
    "        sig=self.window_and_pad(w)\n",
    "        \n",
    "        signalR=np.zeros(len(sig) * 2)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "        resultR = dsp.arm_cfft_q15(self._cfftQ15,signalR,0,1)\n",
    "\n",
    "        resultRQ31 = dsp.arm_q15_to_q31(resultR)\n",
    "        \n",
    "        \n",
    "        self._noise = dsp.arm_cmplx_mag_squared_q31(resultRQ31)\n",
    "        \n",
    "        \n",
    "        resultR = np.zeros(len(resultR),dtype=np.int16)\n",
    "        \n",
    "        res = dsp.arm_cfft_q15(self._cfftQ15,resultR,1,1)\n",
    "        res = dsp.arm_shift_q15(res,self._fftShift)\n",
    "        \n",
    "        res=res[0::2]\n",
    "        res=self.remove_padding(res)\n",
    "        \n",
    "        if status == 0:\n",
    "            res=self.undo_scale(res,the_max)\n",
    "            \n",
    "        return(res)\n",
    "    \n",
    "    def do_nothing(self,w):\n",
    "        w,status,the_max = self.rescale(w)\n",
    "        sig=self.window_and_pad(w)\n",
    "       \n",
    "        \n",
    "        # Convert to complex        \n",
    "        signalR=np.zeros(len(sig) * 2,dtype=np.int16)\n",
    "        signalR[0::2]=sig\n",
    "        \n",
    "\n",
    "        resultR = dsp.arm_cfft_q15(self._cfftQ15,signalR,0,1)\n",
    "        res = dsp.arm_cfft_q15(self._cfftQ15,resultR,1,1)\n",
    "        res = dsp.arm_shift_q15(res,self._fftShift)\n",
    "        \n",
    "        res=res[0::2]\n",
    "        \n",
    "        res=self.remove_padding(res)\n",
    "        \n",
    "        if status == 0:\n",
    "            res=self.undo_scale(res,the_max)\n",
    "            \n",
    "        return(res)\n",
    "        \n",
    "        \n",
    "    def nr(self,nonr=False):\n",
    "        for (w,v) in zip(self._slices,self._vad):\n",
    "            result=None\n",
    "            if nonr:\n",
    "                result = self.do_nothing(w)\n",
    "            else:\n",
    "                if v==1:\n",
    "                    result=self.remove_noise(w)\n",
    "                else:\n",
    "                    result=self.estimate_noise(w)\n",
    "            self._signal.append(result)\n",
    "        \n",
    "    def overlap_and_add(self):\n",
    "        nbSamples = len(self._signal)*winOverlap\n",
    "        offsets = range(0, nbSamples,winOverlap)\n",
    "        offsets=offsets[0:len(self._signal)]\n",
    "        res=np.zeros(nbSamples,dtype=np.int16)\n",
    "        i=0\n",
    "        for n in offsets:\n",
    "            res[n:n+winLength] = dsp.arm_add_q15(res[n:n+winLength],self._signal[i])\n",
    "            i=i+1\n",
    "        return(res)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cf1b685",
   "metadata": {},
   "source": [
    "Processing with no noise suppression and only a FFT/IFFT. It is useful to check if the Q15 FFT/IFFT is working without any saturation. If the reconstructed signal is bad, then you probably need to decrease the FFT length by decreasing the window duration or the sampling frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "48a8041b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13696 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=NoiseSuppressionQ15(slices_q15)\n",
    "n.nr(nonr=True)\n",
    "cleaned_q15=n.overlap_and_add()\n",
    "plt.plot(fix.Q15toF32(cleaned_q15))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "a8cac01a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480808 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audioQ15_nonr=Audio(data=fix.Q15toF32(cleaned_q15),rate=samplerate,autoplay=False)\n",
    "audioQ15_nonr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9543f723",
   "metadata": {},
   "source": [
    "Same but with noise suppression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "f99e212c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12444 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=NoiseSuppressionQ15(slices_q15)\n",
    "n.nr()\n",
    "cleaned_q15=n.overlap_and_add()\n",
    "plt.plot(fix.Q15toF32(cleaned_q15))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abd36813",
   "metadata": {},
   "source": [
    "The audio to check the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "b31616ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 480808 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "audioQ15=Audio(data=fix.Q15toF32(cleaned_q15),rate=samplerate,autoplay=False)\n",
    "audioQ15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7204cf08",
   "metadata": {},
   "source": [
    "## C Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d83edfd5",
   "metadata": {},
   "source": [
    "Once the Python code is working, the C code should be easy to write since the API of the CMSIS-DSP Python wrapper is as close as possible to the C one.\n",
    "\n",
    "For instance, let's look at the F32 noise suppression:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e99015e9",
   "metadata": {},
   "source": [
    "```python\n",
    "energy = dsp.arm_cmplx_mag_squared_f32(v)\n",
    "# To avoid division by zeros\n",
    "energy = dsp.arm_offset_f32(energy,1e-6)\n",
    "        \n",
    "temp = dsp.arm_sub_f32(energy , self._noise)\n",
    "# C implementation will have to loop on the array and replace negative vaues\n",
    "temp[temp<0]=0\n",
    "            \n",
    "scaling = np.zeros(len(temp),dtype=np.float32)\n",
    "# C implementation will have to loop on the arrays here\n",
    "scaling = temp / energy\n",
    "        \n",
    "res=dsp.arm_cmplx_mult_real_f32(v,scaling) \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd93e60a",
   "metadata": {},
   "source": [
    "The C code could be:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "540d8fb5",
   "metadata": {},
   "source": [
    "```C\n",
    "float32_t *v,*energy,*noise,*temp,*scaling,*res;\n",
    "\n",
    "arm_cmplx_mag_squared_f32(v,energy,windowLength);\n",
    "// To avoid division by zeros\n",
    "arm_offset_f32(energy,1e-6,energy,windowLength);\n",
    "        \n",
    "arm_sub_f32(energy , noise,temp,windowLength);\n",
    "\n",
    "for(int i=0;i < windowLength; i++)\n",
    "{\n",
    "    if (temp[i]<0)\n",
    "    {\n",
    "        temp[i]=0;\n",
    "    }\n",
    "}\n",
    "            \n",
    "for(int i=0;i < windowLength; i++)\n",
    "{\n",
    "  scaling[i] = temp[i] / energy[i];\n",
    "}\n",
    "        \n",
    "arm_cmplx_mult_real_f32(v,scaling,res,windowLength) ;\n",
    "```"
   ]
  }
 ],
 "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
