add inverse DFT

Hadryan · Sep 14, 2019 · 6e3109d · 6e3109d
1 parent 85aeda4
commit 6e3109d
Show file tree

Hide file tree

Showing 6 changed files with 1,262 additions and 74 deletions.
diff --git a/CQT.ipynb b/CQT.ipynb
@@ -2,12 +2,13 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [],
    "source": [
     "import Spectrogram\n",
     "from librosa.core import note_to_hz, cqt\n",
+    "from librosa import filters\n",
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
@@ -16,7 +17,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -37,6 +38,145 @@
     "    + np.sin(2*np.pi*800*s, dtype=np.float32) + np.sin(2*np.pi*1600*s, dtype=np.float32)  "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 31,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def cqt_filter_fft(sr, fmin, n_bins, bins_per_octave, tuning,\n",
+    "                     filter_scale, norm, sparsity, hop_length=None,\n",
+    "                     window='hann'):\n",
+    "    '''Generate the frequency domain constant-Q filter basis.'''\n",
+    "\n",
+    "    basis, lengths = filters.constant_q(sr,\n",
+    "                                        fmin=fmin,\n",
+    "                                        n_bins=n_bins,\n",
+    "                                        bins_per_octave=bins_per_octave,\n",
+    "                                        tuning=tuning,\n",
+    "                                        filter_scale=filter_scale,\n",
+    "                                        norm=norm,\n",
+    "                                        pad_fft=True,\n",
+    "                                        window=window)\n",
+    "\n",
+    "    # Filters are padded up to the nearest integral power of 2\n",
+    "    n_fft = basis.shape[1]\n",
+    "\n",
+    "    if (hop_length is not None and\n",
+    "            n_fft < 2.0**(1 + np.ceil(np.log2(hop_length)))):\n",
+    "\n",
+    "        n_fft = int(2.0 ** (1 + np.ceil(np.log2(hop_length))))\n",
+    "\n",
+    "    # re-normalize bases with respect to the FFT window length\n",
+    "    basis *= lengths[:, np.newaxis] / float(n_fft)\n",
+    "\n",
+    "    # FFT and retain only the non-negative frequencies\n",
+    "    fft = get_fftlib()\n",
+    "    fft_basis = fft.fft(basis, n=n_fft, axis=1)[:, :(n_fft // 2)+1]\n",
+    "\n",
+    "    # sparsify the basis\n",
+    "    fft_basis = util.sparsify_rows(fft_basis, quantile=sparsity)\n",
+    "\n",
+    "    return fft_basis, n_fft, lengths, basis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 32,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sr = fs\n",
+    "fmin_t = fmin\n",
+    "n_filters = 12\n",
+    "bins_per_octave = 12\n",
+    "window= 'hann'\n",
+    "from librosa import get_fftlib\n",
+    "from librosa import util\n",
+    "fft_basis, n_fft, _, basis = cqt_filter_fft(sr, fmin_t,\n",
+    "                                    n_filters,\n",
+    "                                    bins_per_octave,\n",
+    "                                    tuning=0,\n",
+    "                                    filter_scale=1,\n",
+    "                                    norm=None,\n",
+    "                                    sparsity=0,\n",
+    "                                    window=window)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([-0.02515722+0.j        , -0.01870174-0.01804653j,\n",
+       "       -0.00143447-0.03432906j, ...,  0.02076852+0.04632163j,\n",
+       "       -0.00143447+0.03432906j, -0.01870174+0.01804653j])"
+      ]
+     },
+     "execution_count": 41,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.fft.fft(basis[0].real) + np.fft.fft(basis[0].imag)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ 0.00081404+0.j        ,  0.00063303-0.01746888j,\n",
+       "        0.00015239-0.033281j  , ..., -0.00045847+0.04500784j,\n",
+       "        0.00015239+0.033281j  ,  0.00063303+0.01746888j])"
+      ]
+     },
+     "execution_count": 40,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x7fd8a0fa9ac8>]"
+      ]
+     },
+     "execution_count": 30,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEtpJREFUeJzt3X+wndVd7/H3hwSobaeElJRiftyk06gTq73gkVLrvWOltoCdpnemKoxjY2Umc3upt2rHGuSPjjrOWO8dsZ1b0UyJ0ju9UsQqGS4WgbY6/gElsS0lQORIxSTyI7WU2qJAyNc/9jrnbGJCyNk7ZydnvV8ze87zrGed/ax1VmZ/8qznx05VIUnq1ymTboAkabIMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnlk66AS/GWWedVWvXrp10MyTppLJz586vVdWKo9U7KYJg7dq17NixY9LNkKSTSpKHX0w9p4YkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSepct0HwZ1/cy7efPjDpZkjSxHUZBDsf/jq/+Kkv86HtuybdFEmauC6D4FtPPwfAY9/8twm3RJImr8sgkCTNGUsQJFmW5MYkDyS5P8kbkyxPcluSB9vPM1vdJPlokukk9yQ5bxxtOBZVtdC7lKQT1riOCD4CfKaqvgd4PXA/sAW4o6rWA3e0dYCLgfXttRm4ZkxtOGZJJrVrSTphjBwESc4A/itwLUBVPVNV3wA2Ate1atcB72zLG4FP1MCdwLIk54zaDknS/IzjiGAdsB/4wyRfTPLxJC8Dzq6qR1qdR4Gz2/JKYM/Q7+9tZc+TZHOSHUl27N+/fwzNlCQdzjiCYClwHnBNVZ0LfJu5aSAAajApf0wT81W1taqmqmpqxYqjfq/CMfEMgSTNGUcQ7AX2VtVdbf1GBsHw2MyUT/v5eNu+D1g99PurWpkkaQJGDoKqehTYk+S7W9GFwH3AdmBTK9sE3NSWtwPvblcPXQA8OTSFJElaYOP6qsqfBz6Z5DTgIeA9DELmhiSXAw8DP9nq3gJcAkwDT7W6kqQJGUsQVNWXgKnDbLrwMHULuGIc+5Ukja7rO4u9i0CSOg8CSZJBIEndMwgkqXMGgSR1rs8g8NZiSZrVZxA0PnxUkjoPAklSp0FQzg1J0qwug0CSNMcgkKTOGQSS1DmDQJI613UQePWoJHUeBJIkg0CSumcQSFLnDAJJ6lyXQVDeWCxJs7oMAknSHINAkjrXZRDMTA3F51BLUp9BIEmaM7YgSLIkyReT3NzW1yW5K8l0kk8lOa2Vn97Wp9v2teNqgyTp2I3ziOD9wP1D6x8Grq6q1wJPAJe38suBJ1r51a2eJGlCxhIESVYBPw58vK0H+FHgxlblOuCdbXljW6dtvzBO1kvSxIzriOB3gQ8CB9v6K4FvVNWBtr4XWNmWVwJ7ANr2J1v950myOcmOJDv2798/pmYeso/j8q6SdHIZOQiSvB14vKp2jqE9s6pqa1VNVdXUihUrxvnWkqQhS8fwHm8C3pHkEuAlwCuAjwDLkixt/+tfBexr9fcBq4G9SZYCZwD/PIZ2SJLmYeQjgqq6sqpWVdVa4FLgs1X108DngHe1apuAm9ry9rZO2/7ZKh/6IEmTcjzvI/gV4JeSTDM4B3BtK78WeGUr/yVgy3FsgyTpKMYxNTSrqj4PfL4tPwScf5g6/wb8xDj3e6w8/JCkOd5ZLEmdMwgkqXNdB4G3sUlS50EgSTIIJKl7BoEkda7LIPD+NUma02UQzPFssSR1HgSSJINAkjpnEEhS5wwCSepcl0HgNUOSNKfLIJAkzTEIJKlzXQeBD52TpM6DQJJkEEhS9wwCSeqcQSBJnesyCHz4qCTN6TIIZnjRkCSNIQiSrE7yuST3JdmV5P2tfHmS25I82H6e2cqT5KNJppPck+S8UdsgSZq/cRwRHAA+UFUbgAuAK5JsALYAd1TVeuCOtg5wMbC+vTYD14yhDZKkeRo5CKrqkar627b8L8D9wEpgI3Bdq3Yd8M62vBH4RA3cCSxLcs6o7ZiPf3ryXyexW0k6oYz1HEGStcC5wF3A2VX1SNv0KHB2W14J7Bn6tb2tbAENzhbfu++bC7tbSToBjS0Ikrwc+FPgF6rqeZ+wNfiS4GO6VifJ5iQ7kuzYv3//uJopSTrEWIIgyakMQuCTVfXpVvzYzJRP+/l4K98HrB769VWt7HmqamtVTVXV1IoVK8bRTEnSYYzjqqEA1wL3V9XvDG3aDmxqy5uAm4bK392uHroAeHJoCkmStMCWjuE93gT8DPCVJF9qZb8K/BZwQ5LLgYeBn2zbbgEuAaaBp4D3jKENkqR5GjkIqupvOPK9WRcepn4BV4y633GpKuLzqCV1rOs7iwGeO+jzJiT1rfsgMAck9a7LIBh+6NxBn0AnqXNdBsEwg0BS7wwCc0BS57oPAk8WS+pd90FQTg1J6lyXQTD80e8RgaTedRkEw8wBSb0zCJwaktQ5g8AgkNS57oPAcwSSetd9EHhAIKl33QeBRwSSetd9EHiOQFLvugwCHzonSXO6DIJhzgxJ6l33QeA5Akm96z4InBqS1Lsug6CGnjZkDkjqXZ9BMPTh79SQpN51GQTDnBqS1Lsug2D4o98gkNS7iQVBkouS7E4ynWTLpNrhzJCk3k0kCJIsAT4GXAxsAC5LsmGh9j/8rWSeI5DUu0kdEZwPTFfVQ1X1DHA9sHEhdvytpw/w25/ZPbv+1DMHFmK3knTCWjqh/a4E9gyt7wXeMO6dPPnUs/yP/7eTA88Vzx0sDhws9j7xFF/71jOzdT78F7vZ9jf/wDPPHWRJQjLuVkjS/L32VS/n1ze+7rjuY1JBcFRJNgObAdasWTPPN4Gnnz3IklPC6aeewktPOYWpVyznM7sena2y+7F/AWDZS0/lwMGDs+VVGAqSJu7Z547/9PWkgmAfsHpofVUrm1VVW4GtAFNTU/P6S5zxHady43t/6D+Ur93y/2eX//qX38yaV750Pm8vSYvCpM4R3A2sT7IuyWnApcD2STTEEJDUu4kcEVTVgSTvA24FlgDbqmrXJNoiSb2b2DmCqroFuGVS+5ckDXR5Z7EkaY5BIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkda7LIPALZyRpTp9BMOkGSNIJpMsgOMVDAkma1WUQmAOSNKfPIHBySJJmdRkE5oAkzekyCMwBSZrTZRB4sliS5nQZBOaAJM3pMgg8IpCkOV0GgTEgSXNGCoIk/yvJA0nuSfJnSZYNbbsyyXSS3UneNlR+USubTrJllP3Pv+ET2asknZBGPSK4DXhdVX0/8HfAlQBJNgCXAt8LXAT8XpIlSZYAHwMuBjYAl7W6C8ockKQ5IwVBVf1lVR1oq3cCq9ryRuD6qnq6qr4KTAPnt9d0VT1UVc8A17e6C+qUU4wCSZoxznMEPwf8RVteCewZ2ra3lR2pfEEZA5I0Z+nRKiS5HXj1YTZdVVU3tTpXAQeAT46rYUk2A5sB1qxZM663BbxqSJKGHTUIquotL7Q9yc8CbwcurKpqxfuA1UPVVrUyXqD80P1uBbYCTE1N1eHqzJc5IElzRr1q6CLgg8A7quqpoU3bgUuTnJ5kHbAe+AJwN7A+ybokpzE4obx9lDbMs90LvUtJOmEd9YjgKP4PcDpwW/twvbOq/ntV7UpyA3AfgymjK6rqOYAk7wNuBZYA26pq14htOGbGgCTNGSkIquq1L7DtN4HfPEz5LcAto+x3VB4QSNKcTu8sNgkkaUafQWAOSNKsLoPAy0claU6XQbBy2XdMugmSdMLoMgh+/2d+YNJNkKQTRpdBsPxlp026CZJ0wugyCCRJcwwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnxhIEST6QpJKc1daT5KNJppPck+S8obqbkjzYXpvGsX9J0vwtHfUNkqwG3gr841DxxcD69noDcA3whiTLgQ8BU0ABO5Nsr6onRm2HJGl+xnFEcDXwQQYf7DM2Ap+ogTuBZUnOAd4G3FZVX28f/rcBF42hDZKkeRopCJJsBPZV1ZcP2bQS2DO0vreVHalckjQhR50aSnI78OrDbLoK+FUG00Jjl2QzsBlgzZo1x2MXkiReRBBU1VsOV57k+4B1wJeTAKwC/jbJ+cA+YPVQ9VWtbB/wI4eUf/4I+90KbAWYmpqqw9WRJI1u3lNDVfWVqnpVVa2tqrUMpnnOq6pHge3Au9vVQxcAT1bVI8CtwFuTnJnkTAZHE7eO3g1J0nyNfNXQEdwCXAJMA08B7wGoqq8n+Q3g7lbv16vq68epDZKkF2FsQdCOCmaWC7jiCPW2AdvGtV9J0mi8s1iSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzo0cBEl+PskDSXYl+e2h8iuTTCfZneRtQ+UXtbLpJFtG3b8kaTRLR/nlJG8GNgKvr6qnk7yqlW8ALgW+F/hO4PYk39V+7WPAjwF7gbuTbK+q+0ZphyRp/kYKAuC9wG9V1dMAVfV4K98IXN/Kv5pkGji/bZuuqocAklzf6hoEkjQho04NfRfwX5LcleSvkvxgK18J7Bmqt7eVHalckjQhRz0iSHI78OrDbLqq/f5y4ALgB4EbkrxmHA1LshnYDLBmzZpxvOXzXP1Tr+fsV7xk7O8rSSebowZBVb3lSNuSvBf4dFUV8IUkB4GzgH3A6qGqq1oZL1B+6H63AlsBpqam6mjtPFb/7dxV435LSTopjTo19OfAmwHayeDTgK8B24FLk5yeZB2wHvgCcDewPsm6JKcxOKG8fcQ2SJJGMOrJ4m3AtiT3As8Am9rRwa4kNzA4CXwAuKKqngNI8j7gVmAJsK2qdo3YBknSCDL43D6xTU1N1Y4dOybdDEk6qSTZWVVTR6vnncWS1DmDQJI6ZxBIUucMAknqnEEgSZ07Ka4aSrIfeHiEtziLwf0NPeipr9BXf+3r4nW8+vufqmrF0SqdFEEwqiQ7XswlVItBT32FvvprXxevSffXqSFJ6pxBIEmd6yUItk66AQuop75CX/21r4vXRPvbxTkCSdKR9XJEIEk6gkUdBEkuSrI7yXSSLZNuz6iSrE7yuST3JdmV5P2tfHmS25I82H6e2cqT5KOt//ckOW+yPZifJEuSfDHJzW19XftWvOkkn2qPNKc99vxTrfyuJGsn2e5jlWRZkhuTPJDk/iRvXMxjm+QX27/je5P8cZKXLJaxTbItyePtycwzZcc8lkk2tfoPJtl0vNq7aIMgyRLgY8DFwAbgsiQbJtuqkR0APlBVGxh8K9wVrU9bgDuqaj1wR1uHQd/Xt9dm4JqFb/JYvB+4f2j9w8DVVfVa4Ang8lZ+OfBEK7+61TuZfAT4TFV9D/B6Bn1elGObZCXwP4Gpqnodg8fSX8riGds/Ai46pOyYxjLJcuBDwBsYfOf7h2bCY+yqalG+gDcCtw6tXwlcOel2jbmPNwE/BuwGzmll5wC72/IfAJcN1Z+td7K8GHyL3R3AjwI3A2Fw483SQ8eZwfdcvLEtL231Muk+vMh+ngF89dD2LtaxZe77y5e3sboZeNtiGltgLXDvfMcSuAz4g6Hy59Ub52vRHhEw9w9txt5Wtii0Q+NzgbuAs6vqkbbpUeDstrwY/ga/C3wQONjWXwl8o6oOtPXhPs32t21/stU/GawD9gN/2KbBPp7kZSzSsa2qfcD/Bv4ReITBWO1kcY7tjGMdywUb48UcBItWkpcDfwr8QlV9c3hbDf7rsCguBUvyduDxqto56bYsgKXAecA1VXUu8G3mpg6ARTe2ZwIbGQTgdwIv4z9OpSxaJ9pYLuYg2AesHlpf1cpOaklOZRACn6yqT7fix5Kc07afAzzeyk/2v8GbgHck+QfgegbTQx8BliWZ+ZrV4T7N9rdtPwP454Vs8Aj2Anur6q62fiODYFisY/sW4KtVtb+qngU+zWC8F+PYzjjWsVywMV7MQXA3sL5dhXAagxNR2yfcppEkCXAtcH9V/c7Qpu3AzBUFmxicO5gpf3e7KuEC4MmhQ9MTXlVdWVWrqmotg/H7bFX9NPA54F2t2qH9nfk7vKvVP2H+1/VCqupRYE+S725FFzL4zu9FObYMpoQuSPLS9u96pr+LbmyHHOtY3gq8NcmZ7Qjqra1s/CZ9QuU4n6y5BPg74O+BqybdnjH054cZHE7eA3ypvS5hMFd6B/AgcDuwvNUPgyun/h74CoMrNCbej3n2/UeAm9vya4AvANPAnwCnt/KXtPXptv01k273MfbxPwM72vj+OXDmYh5b4NeAB4B7gf8LnL5Yxhb4YwbnPp5lcLR3+XzGEvi51udp4D3Hq73eWSxJnVvMU0OSpBfBIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6pxBIEmdMwgkqXP/DuFepEEaYgwfAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(fft_basis.toarray()[0].real)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},